2 Computing the self-force: a 2010 literature survey

Much progress has been achieved in the development of practical methods for computing the self-force. We briefly summarize these efforts in this section, with the goal of introducing the main ideas and some key issues. A more detailed coverage of the various implementations can be found in Barack’s excellent review [9Jump To The Next Citation Point]. The 2005 collection of reviews published in Classical and Quantum Gravity [118] is also recommended for an introduction to the various aspects of self-force theory and numerics. Among our favourites in this collection are the reviews by Detweiler [49Jump To The Next Citation Point] and Whiting [183].

An important point to bear in mind is that all the methods covered here mainly compute the self-force on a particle moving on a fixed world line of the background spacetime. A few numerical codes based on the radiative approximation have allowed orbits to evolve according to energy and angular-momentum balance. As will be emphasized below, however, these calculations miss out on important conservative effects that are only accounted for by the full self-force. Work is currently underway to develop methods to let the self-force alter the motion of the particle in a self-consistent manner.

2.1 Early work: DeWitt and DeWitt; Smith and Will

The first evaluation of the electromagnetic self-force in curved spacetime was carried out by DeWitt and DeWitt [132Jump To The Next Citation Point] for a charge moving freely in a weakly curved spacetime characterized by a Newtonian potential Φ ≪ 1. In this context the right-hand side of Eq. (1.33View Equation) reduces to the tail integral, because the particle moves in a vacuum region of the spacetime, and there is no external force acting on the charge. They found that the spatial components of the self-force are given by

2M-- 2- 2dg- fem = e r3 ˆr + 3 e dt, (2.1 )
where M is the total mass contained in the spacetime, r = |x| is the distance from the centre of mass, ˆr = x∕r, and g = − ∇ Φ is the Newtonian gravitational field. (In these expressions the bold-faced symbols represent vectors in three-dimensional flat space.) The first term on the right-hand side of Eq. (2.1View Equation) is a conservative correction to the Newtonian force mg. The second term is the standard radiation-reaction force; although it comes from the tail integral, this is the same result that would be obtained in flat spacetime if an external force mg were acting on the particle. This agreement is necessary, but remarkable!

A similar expression was obtained by Pfenning and Poisson [141Jump To The Next Citation Point] for the case of a scalar charge. Here

2M-- 1- 2dg- fscalar = 2ξq r3 ˆr + 3 q dt , (2.2 )
where ξ is the coupling of the scalar field to the spacetime curvature; the conservative term disappears when the field is minimally coupled. Pfenning and Poisson also computed the gravitational self-force acting on a point mass moving in a weakly curved spacetime. The expression they obtained is in complete agreement (within its domain of validity) with the standard post-Newtonian equations of motion.

The force required to hold an electric charge in place in a Schwarzschild spacetime was computed, without approximations, by Smith and Will [163]. As we reviewed previously in Section 1.10, the self-force contribution to the total force is given by

r 2M-- 1∕2 fself = e r3 f , (2.3 )
where M is the black-hole mass, r the position of the charge (in Schwarzschild coordinates), and f := 1 − 2M ∕r. When r ≫ M, this expression agrees with the conservative term in Eq. (2.1View Equation). This result was generalized to Reissner–Nordström spacetime by Zel’nikov and Frolov [186]. Wiseman [185] calculated the self-force acting on a static scalar charge in Schwarzschild spacetime. He found that in this case the self-force vanishes. This result is not incompatible with Eq. (2.2View Equation), even for nonminimal coupling, because the computation of the weak-field self-force requires the presence of matter, while Wiseman’s scalar charge lives in a purely vacuum spacetime.

2.2 Mode-sum method

Self-force calculations involving a sum over modes were pioneered by Barack and Ori [16, 7], and the method was further developed by Barack, Ori, Mino, Nakano, and Sasaki [15Jump To The Next Citation Point, 8, 18Jump To The Next Citation Point, 20Jump To The Next Citation Point, 19Jump To The Next Citation Point, 127Jump To The Next Citation Point]; a somewhat related approach was also considered by Lousto [117]. It has now emerged as the method of choice for self-force calculations in spacetimes such as Schwarzschild and Kerr. Our understanding of the method was greatly improved by the Detweiler–Whiting decomposition [53Jump To The Next Citation Point] of the retarded field into singular and regular pieces, as outlined in Sections 1.4 and 1.8, and subsequent work by Detweiler, Whiting, and their collaborators [51Jump To The Next Citation Point].

Detweiler–Whiting decomposition; mode decomposition; regularization parameters

For simplicity we consider the problem of computing the self-force acting on a particle with a scalar charge q moving on a world line γ. (The electromagnetic and gravitational problems are conceptually similar, and they will be discussed below.) The potential Φ produced by the particle satisfies Eq. (1.34View Equation), which we rewrite schematically as

□Φ = qδ(x,z ), (2.4 )
where □ is the wave operator in curved spacetime, and δ(x,z) represents a distributional source that depends on the world line γ through its coordinate representation z (τ). From the perspective of the Detweiler–Whiting decomposition, the scalar self-force is given by
F = q∇ Φ := q(∇ Φ − ∇ Φ ), (2.5 ) α α R α α S
where Φ, ΦS, and ΦR are the retarded, singular, and regular potentials, respectively. To evaluate the self-force, then, is to compute the gradient of the regular potential.

From the point of view of Eq. (2.5View Equation), the task of computing the self-force appears conceptually straightforward: Either (i) compute the retarded and singular potentials, subtract them, and take a gradient of the difference; or (ii) compute the gradients of the retarded and singular potentials, and then subtract the gradients. Indeed, this is the basic idea for most methods of self-force computations. However, the apparent simplicity of this sequence of steps is complicated by the following facts: (i) except for a very limited number of cases, the retarded potential of a point particle cannot be computed analytically and must therefore be obtained by numerical means; and (ii) both the retarded and singular potential diverge at the particle’s position. Thus, any sort of subtraction will generally have to be performed numerically, and for this to be possible, one requires representations of the retarded and singular potentials (and/or their gradients) in terms of finite quantities.

In a mode-sum method, these difficulties are overcome with a decomposition of the potential in spherical-harmonic functions:

Φ = ∑ Φlm (t,r)Ylm(𝜃,ϕ ). (2.6 ) lm
When the background spacetime is spherically symmetric, Eq. (2.4View Equation) gives rise to a fully decoupled set of reduced wave equations for the mode coefficients Φlm(t,r), and these are easily integrated with simple numerical methods. The dimensional reduction of the wave equation implies that each Φlm (t,r ) is finite and continuous (though nondifferentiable) at the position of the particle. There is, therefore, no obstacle to evaluating each l-mode of the field, defined by
∑ l lm lm (∇ αΦ )l := lim ∇α [Φ (t,r)Y (𝜃,ϕ )]. (2.7 ) x→z m=− l
The sum over modes, however, must reproduce the singular field evaluated at the particle’s position, and this is infinite; the mode sum, therefore, does not converge.

Fortunately, there is a piece of each l-mode that does not contribute to the self-force, and that can be subtracted out; this piece is the corresponding l-mode of the singular field ∇α ΦS. Because the retarded and singular fields share the same singularity structure near the particle’s world line (as described in Section 1.6), the subtraction produces a mode decomposition of the regular field ∇ αΦR. And because this field is regular at the particle’s position, the sum over all modes q(∇ αΦR )l is guaranteed to converge to the correct value for the self-force. The key to the mode-sum method, therefore, is the ability to express the singular field as a mode decomposition.

This can be done because the singular field, unlike the retarded field, can always be expressed as a local expansion in powers of the distance to the particle; such an expansion was displayed in Eqs. (1.28View Equation) and (1.29View Equation). (In a few special cases the singular field is actually known exactly [43, 114, 33Jump To The Next Citation Point, 86Jump To The Next Citation Point, 162].) This local expansion can then be turned into a multipole decomposition. Barack and Ori [18, 15, 20, 19Jump To The Next Citation Point, 9], and then Mino, Nakano, and Sasaki [127], were the first to show that this produces the following generic structure:

C D E (∇ αΦS )l = (l + 12)A α + B α + --α-1+ ----1--α---3- + -----3-----1-α----3-----5- + ⋅⋅⋅, (2.8 ) l + 2 (l − 2)(l + 2) (l − 2)(l − 2)(l + 2)(l + 2)
where A α, B α, C α, and so on are l-independent functions that depend on the choice of field (i.e., scalar, electromagnetic, or gravitational), the choice of spacetime, and the particle’s state of motion. These so-called regularization parameters are now ubiquitous in the self-force literature, and they can all be determined from the local expansion for the singular field. The number of regularization parameters that can be obtained depends on the accuracy of the expansion. For example, expansions accurate through order 0 r such as Eqs. (1.28View Equation) and (1.29View Equation) permit the determination of Aα, B α, and Cα; to obtain D α one requires the terms of order r, and to get Eα the expansion must be carried out through order r2. The particular polynomials in l that accompany the regularization parameters were first identified by Detweiler and his collaborators [51Jump To The Next Citation Point]. Because the D α term is generated by terms of order r in the local expansion of the singular field, the sum of 1 3 − 1 [(l − 2 )(l + 2)] from l = 0 to l = ∞ evaluates to zero. The sum of the polynomial in front of E α also evaluates to zero, and this property is shared by all remaining terms in Eq. (2.8View Equation).

Mode sum

With these elements in place, the self-force is finally computed by implementing the mode-sum formula

∑ L [ F α = q (∇ αΦ)l − (l + 1)A α − B α −-Cα--− -----D-α----- 2 l + 12 (l − 12)(l + 32) l=0 ] − ------------Eα------------ − ⋅⋅⋅ + remainder, (2.9 ) (l − 32)(l − 12)(l + 32)(l + 52)
where the infinite sum over l is truncated to a maximum mode number L. (This truncation is necessary in practice, because in general the modes must be determined numerically.) The remainder consists of the remaining terms in the sum, from l = L + 1 to l = ∞; it is easy to see that since the next regularization term would scale as l−6 for large l, the remainder scales as L−5, and can be made negligible by summing to a suitably large value of l. This observation motivates the inclusion of the D α and Eα terms within the mode sum, even though their complete sums evaluate to zero. These terms are useful because the sum must necessarily be truncated, and they permit a more rapid convergence of the mode sum. For example, exclusion of the D α and Eα terms in Eq. (2.9View Equation) would produce a remainder that scales as L−1 instead of L−5; while this is sufficient for convergence, the rate of convergence is too slow to permit high-accuracy computations. Rapid convergence therefore relies on a knowledge of as many regularization parameters as possible, but unfortunately these parameters are not easy to calculate. To date, only A α, B α, C α, and D α have been calculated for general orbits in Schwarzschild spacetime [51Jump To The Next Citation Point, 87Jump To The Next Citation Point], and only Aα, B α, Cα have been calculated for orbits in Kerr spacetime [19]. It is possible, however, to estimate a few additional regularization parameters by fitting numerical results to the structure of Eq. (2.8View Equation); this clever trick was first exploited by Detweiler and his collaborators [51Jump To The Next Citation Point] to achieve extremely high numerical accuracies. This trick is now applied routinely in mode-sum computations of the self-force.

Case study: static electric charge in extreme Reissner–Nordström spacetime

The practical use of the mode-sum method can be illustrated with the help of a specific example that can be worked out fully and exactly. We consider, as in Section 1.10, an electric charge e held in place at position r = r0 in the spacetime of an extreme Reissner–Nordström black hole of mass M and charge Q = M. The reason for selecting this spacetime resides in the resulting simplicity of the spherical-harmonic modes for the electromagnetic field.

The metric of the extreme Reissner–Nordström spacetime is given by

ds2 = − f dt2 + f −1dr2 + r2dΩ2, (2.10 )
where 2 f = (1 − M ∕r). The only nonzero component of the electromagnetic field tensor is Ftr = − Er, and this is decomposed as
∑ Ftr = Fltmr (r)Y lm(𝜃,ϕ). (2.11 ) lm
This field diverges at r = r0, but the modes lm Ftr (r) are finite, though discontinuous. The multipole coefficients of the field are defined to be
∑l (Ftr)l = lim Ftlmr Y lm, (2.12 ) m= −l
where the limit is taken in the direction of the particle’s position. The charge can be placed on the axis 𝜃 = 0, and this choice produces an axisymmetric field with contributions from m = 0 only. Because Y l0 = [(2l + 1)∕4π]1∕2Pl(cos 𝜃) and Pl(1) = 1, we have
∘ ------- 2l + 1 (Ftr)l = ------ lim Flt0r(r0 + Δ ). (2.13 ) 4π Δ →0
The sign of Δ is arbitrary, and (Ftr)l depends on the direction in which r0 is approached.

The charge density of a static particle can also be decomposed in spherical harmonics, and the mode coefficients are given by

∘ ------- r2jl0 = e 2l +-1-f δ(r − r), (2.14 ) t 4π 0 0
where f = (1 − M ∕r )2 0 0. If we let
l 2 l0 Φ := − r Ftr, (2.15 )
then Gauss’s law in the extreme Reissner–Nordström spacetime can be shown to reduce to
∘ ------- l(l + 1 ) 2l + 1 (f Φ′)′ −----2---Φ = 4πe ------f0δ′(r − r0), (2.16 ) r 4π
in which a prime indicates differentiation with respect to r, and the index l on Φ is omitted to simplify the expressions. The solution to Eq. (2.16View Equation) can be expressed as Φ(r) = Φ >(r)Θ (r − r0) + Φ <(r)Θ(r0 − r), where Φ> and Φ < are each required to satisfy the homogeneous equation ′ ′ 2 (f Φ ) − l(l + 1)Φ ∕r = 0, as well as the junction conditions
∘ ------- 2l + 1 [Φ] = 4πe ------, [Φ ′] = 0, (2.17 ) 4π
with [Φ ] := Φ> (r0) − Φ <(r0) denoting the jump across r = r0.

For l = 0 the general solution to the homogeneous equation is c r∗ + c 1 2, where c 1 and c 2 are constants and ∗ ∫ −1 r = f dr. The solution for r < r0 must be regular at r = M, and we select Φ < = constant. The solution for r > r0 must produce a field that decays as −2 r at large r, and we again select Φ > = constant. Since each constant is proportional to the total charge enclosed within a sphere of radius r, we arrive at

√ --- Φ< = 0, Φ > = 4πe, (l = 0). (2.18 )
For l ⁄= 0 the solutions to the homogeneous equation are
( ) r − M l( ) Φ < = c1e r--−-M-- lr + M (2.19 ) 0
( ) r0 − M l+1[ ] Φ > = c2e -------- (l + 1)r − M . (2.20 ) r − M
The constants c1 and c2 are determined by the junction conditions, and we get
∘ ------- ∘ ------- 4π 1 4π 1 c1 = − --------, c2 = --------. (2.21 ) 2l + 1 r0 2l + 1 r0
The modes of the electromagnetic field are now completely determined.

According to the foregoing results, and recalling the definition of Eq. (2.13View Equation), the multipole coefficients of the electromagnetic field at + r = r0 + 0 are given by

( ) e ( ) ( ) ( 1) e F >tr 0 = − -2, Ft>r l = e l + 12 − -2 − --3(r0 − 2M ). (2.22 ) r0 r0 2r0
For r = r0 + 0− we have instead
( ) ( < ) ( <) ( 1) 1 e Ftr 0 = 0, F tr l = e l + 2 + -2 − --3(r0 − 2M ). (2.23 ) r0 2r0
We observe that the multipole coefficients lead to a diverging mode sum. We also observe, however, that the multipole structure is identical to the decomposition of the singular field displayed in Eq. (2.8View Equation). Comparison of the two expressions allows us to determine the regularization parameters for the given situation, and we obtain
e e A = ∓ -2, B = − --3-(r0 − 2M ), C = D = E = ⋅⋅⋅ = 0. (2.24 ) r0 2r0
Regularization of the mode sum via Eq. (2.9View Equation) reveals that the modes l ⁄= 0 give rise to the singular field, while the regular field comes entirely from the mode l = 0. In this case, therefore, we can state that the exact expression for the regular field evaluated at the position of the particle is F R = (Ftr)0 − 1A − B tr 2, or FR (r0) = − eM ∕r3 tr 0. Because the regular field must be a solution to the vacuum Maxwell equations, its monopole structure guarantees that its value at any position is given by
eM ∕r FtRr(r) = − ---2-0. (2.25 ) r
This is the field of an image charge e′ = +eM ∕r0 situated at the centre of the black hole.

The self-force acting on the static charge is then

r ∘ --- R e2M--∘ --- e2M-- f = − e f0F tr(r0) = r3 f0 = r3 (1 − M ∕r0). (2.26 ) 0 0
This expression agrees with the Smith-Will force of Eq. (1.50View Equation). The interpretation of the result in terms of an interaction between e and the image charge e′ was elaborated in Sec. 1.10.

Computations in Schwarzschild spacetime

The mode-sum method was successfully implemented in Schwarzschild spacetime to compute the scalar and electromagnetic self-forces on a static particle [31, 36] . It was used to calculate the scalar self-force on a particle moving on a radial trajectory [10], circular orbit [30, 51, 87, 37Jump To The Next Citation Point], and a generic bound orbit [84Jump To The Next Citation Point]. It was also developed to compute the electromagnetic self-force on a particle moving on a generic bound orbit [85Jump To The Next Citation Point], as well as the gravitational self-force on a point mass moving on circular [21Jump To The Next Citation Point, 1] and eccentric orbits [22Jump To The Next Citation Point]. The mode-sum method was also used to compute unambiguous physical effects associated with the gravitational self-force [50Jump To The Next Citation Point, 157, 11Jump To The Next Citation Point], and these results were involved in detailed comparisons with post-Newtonian theory [50Jump To The Next Citation Point, 29Jump To The Next Citation Point, 28Jump To The Next Citation Point, 44Jump To The Next Citation Point, 11Jump To The Next Citation Point]. These achievements will be described in more detail in Section 2.6.

An issue that arises in computations of the electromagnetic and gravitational self-forces is the choice of gauge. While the self-force formalism is solidly grounded in the Lorenz gauge (which allows the formulation of a wave equation for the potentials, the decomposition of the retarded field into singular and regular pieces, and the computation of regularization parameters), it is often convenient to carry out the numerical computations in other gauges, such as the popular Regge–Wheeler gauge and the Chrzanowski radiation gauge described below. Compatibility of calculations carried out in different gauges has been debated in the literature. It is clear that the singular field is gauge invariant when the transformation between the Lorenz gauge and the adopted gauge is smooth on the particle’s world line; in such cases the regularization parameters also are gauge invariant [17Jump To The Next Citation Point], the transformation affects the regular field only, and the self-force changes according to Eq. (1.49View Equation). The transformations between the Lorenz gauge and the Regge–Wheeler and radiation gauges are not regular on the world line, however, and in such cases the regularization of the retarded field must be handled with extreme care.

Computations in Kerr spacetime; metric reconstruction

The reliance of the mode-sum method on a spherical-harmonic decomposition makes it generally impractical to apply to self-force computations in Kerr spacetime. Wave equations in this spacetime are better analyzed in terms of a spheroidal-harmonic decomposition, which simultaneously requires a Fourier decomposition of the field’s time dependence. (The eigenvalue equation for the angular functions depends on the mode’s frequency.) For a static particle, however, the situation simplifies, and Burko and Liu [35] were able to apply the method to calculate the self-force on a static scalar charge in Kerr spacetime. More recently, Warburton and Barack [181] carried out a mode-sum calculations of the scalar self-force on a particle moving on equatorial orbits of a Kerr black hole. They first solve for the spheroidal multipoles of the retarded potential, and then re-express them in terms of spherical-harmonic multipoles. Fortunately, they find that a spheroidal multipole is well represented by summing over a limited number of spherical multipoles. The Warburton–Barack work represents the first successful computations of the self-force in Kerr spacetime, and it reveals the interesting effect of the black hole’s spin on the behaviour of the self-force.

The analysis of the scalar wave equation in terms of spheroidal functions and a Fourier decomposition permits a complete separation of the variables. For decoupling and separation to occur in the case of a gravitational perturbation, it is necessary to formulate the perturbation equations in terms of Newman–Penrose (NP) quantities [172], and to work with the Teukolsky equation that governs their behaviour. Several computer codes are now available that are capable of integrating the Teukolsky equation when the source is a point mass moving on an arbitrary geodesic of the Kerr spacetime. (A survey of these codes is given below.) Once a solution to the Teukolsky equation is at hand, however, there still remains the additional task of recovering the metric perturbation from this solution, a problem referred to as metric reconstruction.

Reconstruction of the metric perturbation from solutions to the Teukolsky equation was tackled in the past in the pioneering efforts of Chrzanowski [41], Cohen and Kegeles [42, 105], Stewart [166], and Wald [179]. These works have established a procedure, typically attributed to Chrzanowski, that returns the metric perturbation in a so-called radiation gauge. An important limitation of this method, however, is that it applies only to vacuum solutions to the Teukolsky equation. This makes the standard Chrzanowski procedure inapplicable in the self-force context, because a point particle must necessarily act as a source of the perturbation. Some methods were devised to extend the Chrzanowski procedure to accommodate point sources in specific circumstances [121, 134], but these were not developed sufficiently to permit the computation of a self-force. See Ref. [184] for a review of metric reconstruction from the perspective of self-force calculations.

A remarkable breakthrough in the application of metric-reconstruction methods in self-force calculations was achieved by Keidl, Wiseman, and Friedman [107Jump To The Next Citation Point, 106, 108Jump To The Next Citation Point], who were able to compute a self-force starting from a Teukolsky equation sourced by a point particle. They did it first for the case of an electric charge and a point mass held at a fixed position in a Schwarzschild spacetime [107Jump To The Next Citation Point], and then for the case of a point mass moving on a circular orbit around a Schwarzschild black hole [108Jump To The Next Citation Point]. The key conceptual advance is the realization that, according to the Detweiler–Whiting perspective, the self-force is produced by a regularized field that satisfies vacuum field equations in a neighbourhood of the particle. The regular field can therefore be submitted to the Chrzanowski procedure and reconstructed from a source-free solution to the Teukolsky equation.

More concretely, suppose that we have access to the Weyl scalar ψ 0 produced by a point mass moving on a geodesic of a Kerr spacetime. To compute the self-force from this, one first calculates the singular Weyl scalar ψS0 from the Detweiler–Whiting singular field hSαβ, and subtracts it from ψ0. The result is a regularized Weyl scalar ψR0, which is a solution to the homogeneous Teukolsky equation. This sets the stage for the metric-reconstruction procedure, which returns (a piece of) the regular field R h αβ in the radiation gauge selected by Chrzanowski. The computation must be completed by adding the pieces of the metric perturbation that are not contained in ψ0; these are referred to either as the nonradiative degrees of freedom (since ψ0 is purely radiative), or as the l = 0 and l = 1 field multipoles (because the sum over multipoles that make up ψ0 begins at l = 2). A method to complete the Chrzanowski reconstruction of hR αβ was devised by Keidl et al. [107, 108], and the end result leads directly to the gravitational self-force. The relevance of the l = 0 and l = 1 modes to the gravitational self-force was emphasized by Detweiler and Poisson [52].

Time-domain versus frequency-domain methods

When calculating the spherical-harmonic components lm Φ (t,r) of the retarded potential Φ – refer back to Eq. (2.6View Equation) – one can choose to work either directly in the time domain, or perform a Fourier decomposition of the time dependence and work instead in the frequency domain. While the time-domain method requires the integration of a partial differential equation in t and r, the frequency-domain method gives rise to set of ordinary differential equations in r, one for each frequency ω. For particles moving on circular or slightly eccentric orbits in Schwarzschild spacetime, the frequency spectrum is limited to a small number of discrete frequencies, and a frequency-domain method is easy to implement and yields highly accurate results. As the orbital eccentricity increases, however, the frequency spectrum broadens, and the computational burden of summing over all frequency components becomes more significant. Frequency-domain methods are less efficient for large eccentricities, the case of most relevance for extreme-mass-ratio inspirals, and it becomes advantageous to replace them with time-domain methods. (See Ref. [25] for a quantitative study of this claim.) This observation has motivated the development of accurate evolution codes for wave equations in 1+1 dimensions.

Such codes must be able to accommodate point-particle sources, and various strategies have been pursued to represent a Dirac distribution on a numerical grid, including the use of very narrow Gaussian pulses [116, 110, 34] and of “finite impulse representations” [168]. These methods do a good job with waveform and radiative flux calculations far away from the particle, but are of very limited accuracy when computing the potential in a neighborhood of the particle. A numerical method designed to provide an exact representation of a Dirac distribution in a time-domain computation was devised by Lousto and Price [120] (see also Ref. [123]). It was implemented by Haas [84, 85] for the specific purpose of evaluating Φlm (t,r) at the position of particle and computing the self-force. Similar codes were developed by other workers for scalar [176Jump To The Next Citation Point] and gravitational [21, 22Jump To The Next Citation Point] self-force calculations.

Most extant time-domain codes are based on finite-difference techniques, but codes based on pseudo-spectral methods have also been developed [67, 68, 37, 38]. Spectral codes are a powerful alternative to finite-difference codes, especially when dealing with smooth functions, because they produce much faster convergence. The fact that self-force calculations deal with point sources and field modes that are not differentiable might suggest that spectral convergence should not be expected in this case. This objection can be countered, however, by placing the particle at the boundary between two spectral domains. Functions are then smooth in each domain, and discontinuities are handled by formulating appropriate boundary conditions; spectral convergence is thereby achieved.

2.3 Effective-source method

The mode-sum methods reviewed in the preceding subsection have been developed and applied extensively, but they do not exhaust the range of approaches that may be exploited to compute a self-force. Another set of methods, devised by Barack and his collaborators [12Jump To The Next Citation Point, 13Jump To The Next Citation Point, 60Jump To The Next Citation Point] as well as Vega and his collaborators [176Jump To The Next Citation Point, 177Jump To The Next Citation Point, 175], begin by recognizing that an approximation to the exact singular potential can be used to regularize the delta-function source term of the original field equation. We shall explain this idea in the simple context of a scalar potential Φ.

We continue to write the wave equation for the retarded potential Φ in the schematic form

□Φ = qδ(x,z ), (2.27 )
where □ is the wave operator in curved spacetime, and δ(x,z) is a distributional source term that depends on the particle’s world line γ through its coordinate representation z(τ). By construction, the exact singular potential ΦS satisfies the same equation, and an approximation to the singular potential, denoted &tidle;ΦS, will generally satisfy an equation of the form
&tidle; n □ ΦS = qδ(x, x0) + O (r ) (2.28 )
for some integer n > 0, where r is a measure of distance to the world line. A “better” approximation to the singular potential is one with a higher value of n. From the approximated singular potential we form an approximation to the regular potential by writing
&tidle; &tidle; ΦR := Φ − W ΦS, (2.29 )
where W is a window function whose properties will be specified below. The approximated regular potential is governed by the wave equation
( ) □ &tidle;ΦR = qδ(x,z) − □ W &tidle;ΦS := S (x,z), (2.30 )
and the right-hand side of this equation defines the effective source term S(x,z). This equation is much less singular than Eq. (2.27View Equation), and it can be integrated using numerical methods designed to handle smooth functions.

To see this, we write the effective source more specifically as

S (x,z) = − &tidle;Φ □W − 2∇ W ∇ α&tidle;Φ − W □ &tidle;Φ + qδ(x,z). (2.31 ) S α S S
With the window function W designed to approach unity as x → z, we find that the delta function that arises from the third term on the right-hand side precisely cancels out the fourth term. To keep the other terms in S well behaved on the world line, we further restrict the window function to satisfy p ∇ αW = O(r ) with p ≥ 2; this ensures that multiplication by &tidle; −2 ∇ αΦS = O (r ) leaves behind a bounded quantity. In addition, we demand that □W = O(rq) with q ≥ 1, so that multiplication by &tidle;Φ = O (r−1) S again produces a bounded quantity. It is also useful to require that W (x) have compact (spatial) support, to ensure that the effective source term S(x,z ) does not extend beyond a reasonably small neighbourhood of the world line; this property also has the virtue of making &tidle;ΦR precisely equal to the retarded potential Φ outside the support of the window function. This implies, in particular, that Φ&tidle;R can be used directly to compute radiative fluxes at infinity. Another considerable virtue of these specifications for the window function is that they guarantee that the gradient of &tidle; ΦR is directly tied to the self-force. We indeed see that
( ) lim ∇ αΦ&tidle;R = lim ∇α Φ − W ∇ α &tidle;ΦS − lim &tidle;ΦS∇ αW x→z x→z ( ) x→z = lim ∇α Φ − ∇ α&tidle;ΦS x→z = q− 1F α, (2.32 )
with the second line following by virtue of the imposed conditions on W, and the third line from the properties of the approximated singular field.

The effective-source method therefore consists of integrating the wave equation

□ &tidle;Φ = S (x, z), (2.33 ) R
for the approximated regular potential &tidle;ΦR, with a source term S (x,z) that has become a regular function (of limited differentiability) of the spacetime coordinates x. The method is also known as a “puncture approach,” in reference to a similar regularization strategy employed in numerical relativity. It is well suited to a 3+1 integration of the wave equation, which can be implemented on mature codes already in circulation within the numerical-relativity community. An important advantage of a 3+1 implementation is that it is largely indifferent to the choice of background spacetime, and largely insensitive to the symmetries possessed by this spacetime; a self-force in Kerr spacetime is in principle just as easy to obtain as a self-force in Schwarzschild spacetime.

The method is also well suited to a self-consistent implementation of the self-force, in which the motion of the particle is not fixed in advance, but determined by the action of the computed self-force. This amounts to combining Eq. (2.33View Equation) with the self-force equation

Du μ ( μν μ ν) m -----= q g + u u ∇ ν&tidle;ΦR, (2.34 ) dτ
in which the field is evaluated on the dynamically determined world line. The system of equations is integrated jointly, and self-consistently. The 3+1 version of the effective-source approach presents a unique opportunity for the numerical-relativity community to get involved in self-force computations, with only a minimal amount of infrastructure development. This was advocated by Vega and Detweiler [176], who first demonstrated the viability of the approach with a 1+1 time-domain code for a scalar charge on a circular orbit around a Schwarzschild black hole. An implementation with two separate 3+1 codes imported from numerical relativity was also accomplished [177].

The work of Barack and collaborators [12, 13] is a particular implementation of the effective-source approach in a 2+1 numerical calculation of the scalar self-force in Kerr spacetime. (See also the independent implementation by Lousto and Nakano [119].) Instead of starting with Eq. (2.27View Equation), they first decompose Φ according to

Φ(x) = ∑ Φm (t,r,𝜃)exp (im ϕ) (2.35 ) m
and formulate reduced wave equations for the Fourier coefficients Φm. Each coefficient is then regularized with an appropriate singular field &tidle;ΦmS, which eliminates the delta-function from Eq. (2.27View Equation). This gives rise to regularized source terms for the reduced wave equations, which can then be integrated with a 2+1 evolution code. In the final stage of the computation, the self-force is recovered by summing over the regularized Fourier coefficients. This strategy, known as the m-mode regularization scheme, is currently under active development. Recently it was successfully applied by Dolan and Barack [60] to compute the self-force on a scalar charge in circular orbit around a Schwarzschild black hole.

2.4 Quasilocal approach with “matched expansions”

As was seen in Eqs. (1.33View Equation), (1.40View Equation), and (1.47View Equation), the self-force can be expressed as an integral over the past world line of the particle, the integrand involving the Green’s function for the appropriate wave equation. Attempts have been made to compute the Green’s function directly [132Jump To The Next Citation Point, 141Jump To The Next Citation Point, 33Jump To The Next Citation Point, 86Jump To The Next Citation Point], and to evaluate the world-line integral. The quasilocal approach, first introduced by Anderson and his collaborators [4, 3, 6, 5], is based on the expectation that the world-line integral might be dominated by the particle’s recent past, so that the Green’s function can be represented by its Hadamard expansion, which is restricted to the normal convex neighbourhood of the particle’s current position. To help with this enterprise, Ottewill and his collaborators [136, 182, 137, 39] have pushed the Hadamard expansion to a very high order of accuracy, building on earlier work by Décanini and Folacci [48].

The weak-field calculations performed by DeWitt and DeWitt [132] and Pfenning and Poisson [141] suggest that the world-line integral is not, in fact, dominated by the recent past. Instead, most of the self-force is produced by signals that leave the particle at some time in the past, scatter off the central mass, and reconnect with the particle at the current time; such signals mark the boundary of the normal convex neighbourhood. The quasilocal evaluation of the world-line integral must therefore be supplemented with contributions from the distant past, and this requires a representation of the Green’s function that is not limited to the normal convex neighbourhood. In some spacetimes it is possible to express the Green’s function as an expansion in quasi-normal modes, as was demonstrated by Casals and his collaborators for a static scalar charge in the Nariai spacetime [40]. Their study provided significant insights into the geometrical structure of Green’s functions in curved spacetime, and increased our understanding of the non-local character of the self-force.

2.5 Adiabatic approximations

The accurate computation of long-term waveforms from extreme-mass-ratio inspirals (EMRIs) involves a lengthy sequence of calculations that include the calculation of the self-force. One can already imagine the difficulty of numerically integrating the coupled linearized Einstein equation for durations much longer than has ever been attempted by existing numerical codes. While doing so, the code would also have to evaluate the self-force reasonably often (if not at each time step) in order to remain close to the true dynamics of the point mass. Moreover, gravitational-wave data analysis via matched filtering require full evolutions of the sort just described for a large sample of systems parameters. All these considerations underlie the desire for simplified approximations to fully self-consistent self-force EMRI models.

The adiabatic approximation refers to one such class of potentially useful approximations. The basic assumption is that the secular effects of the self-force occur on a timescale that is much longer than the orbital period. In an extreme-mass-ratio binary, this assumption is valid during the early stage of inspiral; it breaks down in the final moments, when the orbit transitions to a quasi-radial infall called the plunge. From the adiabaticity assumption, numerous approximations have been formulated: For example, (i) since the particle’s orbit deviates only slowly from geodesic motion, the self-force can be calculated from a field sourced by a geodesic; (ii) since the radiation-reaction timescale trr, over which a significant shrinking of the orbit occurs due to the self-force, is much longer than the orbital period, periodic effects of the self-force can be neglected; and (iii) conservative effects of the self-force can be neglected (the radiative approximation).

A seminal example of an adiabatic approximation is the Peters–Mathews formalism [140, 139], which determines the long-term evolution of a binary orbit by equating the time-averaged rate of change of the orbital energy E and angular momentum L to, respectively, the flux of gravitational-wave energy and angular momentum at infinity. This formalism was used to successfully predict the decreasing orbital period of the Hulse–Taylor pulsar, before more sophisticated methods, based on post-Newtonian equations of motion expanded to 2.5pn order, were incorporated in times-of-arrival formulae.

In the hope of achieving similar success in the context of the self-force, considerable work has been done to formulate a similar approximation for the case of an extreme-mass-ratio inspiral [124Jump To The Next Citation Point, 125, 126, 98, 61Jump To The Next Citation Point, 62Jump To The Next Citation Point, 159Jump To The Next Citation Point, 158Jump To The Next Citation Point, 78Jump To The Next Citation Point, 128, 94Jump To The Next Citation Point]. Bound geodesics in Kerr spacetime are specified by three constants of motion – the energy E, angular momentum L, and Carter constant C. If one could easily calculate the rates of change of these quantities, using a method analogous to the Peters–Mathews formalism, then one could determine an approximation to the long-term orbital evolution of the small body in an EMRI, avoiding the lengthy process of regularization involved in the direct integration of the self-forced equation of motion. In the early 1980s, Gal’tsov [77] showed that the average rates of change of E and L, as calculated from balance equations that assume geodesic source motion, agree with the averaged rates of change induced by a self-force constructed from a radiative Green’s function defined as Grad := 1(G − − G+ ) 2. As discussed in Section 1.4, this is equal to the regular two-point function GR in flat spacetime, but Grad ⁄= GR in curved spacetime; because of its time-asymmetry, it is purely dissipative. Mino [124], building on the work of Gal’tsov, was able to show that the true self-force and the dissipative force constructed from Grad cause the same averaged rates of change of all three constants of motion, lending credence to the radiative approximation. Since then, the radiative Green’s function was used to derive explicit expressions for the rates of change of E, L, and C in terms of the particle’s orbit and wave amplitudes at infinity [159, 158, 78], and radiative approximations based on such expressions have been concretely implemented by Drasco, Hughes and their collaborators [99, 61, 62].

The relevance of the conservative part of the self-force – the part left out when using Grad – was analyzed in numerous recent publications [32, 148Jump To The Next Citation Point, 146Jump To The Next Citation Point, 147Jump To The Next Citation Point, 94Jump To The Next Citation Point, 97Jump To The Next Citation Point]. As was shown by Pound et al. [148, 146, 147], neglect of the conservative effects of the self-force generically leads to long-term errors in the phase of an orbit and the gravitational wave it produces. These phasing errors are due to orbital precession and a direct shift in orbital frequency. This shift can be understood by considering a conservative force acting on a circular orbit: the force is radial, it alters the centripetal acceleration, and the frequency associated with a given orbital radius is affected. Despite these errors, a radiative approximation may still suffice for gravitational-wave detection [94Jump To The Next Citation Point]; for circular orbits, which have minimal conservative effects, radiative approximations may suffice even for parameter-estimation [97]. However, at this point in time, these analyses remain inconclusive because they all rely on extrapolations from post-Newtonian results for the conservative part of the self-force. For a more comprehensive discussion of these issues, the reader is referred to Ref. [94Jump To The Next Citation Point, 143Jump To The Next Citation Point].

Hinderer and Flanagan performed the most comprehensive study of these issues [69Jump To The Next Citation Point], utilizing a two-timescale expansion [109Jump To The Next Citation Point, 145Jump To The Next Citation Point] of the field equations and self-forced equations of motion in an EMRI. In this method, all dynamical variables are written in terms of two time coordinates: a fast time t and a slow time &tidle; t := (m ∕M )t. In the case of an EMRI, the dynamical variables are the metric and the phase-space variables of the world line. The fast-time dependence captures evolution on the orbital timescale ∼ M, while the slow-time dependence captures evolution on the radiation-reaction timescale ∼ M 2∕m. One obtains a sequence of fast-time and slow-time equations by expanding the full equations in the limit of small m while treating the two time coordinates as independent. Solving the leading-order fast-time equation, in which &tidle; t is held fixed, yields a metric perturbation sourced by a geodesic, as one would expect from the linearized field equations for a point particle. The leading-order effects of the self-force are incorporated by solving the slow-time equation and letting &tidle;t vary. (Solving the next-higher-order slow-time equation determines similar effects, but also the backreaction that causes the parameters of the large black hole to change slowly.)

Using this method, Hinderer and Flanagan identified the effects of the various pieces of the self-force. To describe this we write the self-force as

( ) 2( ) f μ = m-- fμ + fμ + -m-- f μ + fμ + ⋅⋅⋅, (2.36 ) M (1)rr (1)c M 2 (2)rr (2)c
where ‘rr’ denotes a radiation-reaction, or dissipative, piece of the force, and ‘c’ denotes a conservative piece. Hinderer and Flanagan’s principal result is a formula for the orbital phase (which directly determines the phase of the emitted gravitational waves) in terms of these quantities:
2( ) ϕ = M--- ϕ(0)(&tidle;t) + m-ϕ (1)(&tidle;t) + ⋅⋅⋅ , (2.37 ) m M
where ϕ(0) depends on an averaged piece of f μ (1)rr, while ϕ(1) depends on f μ (1)c, the oscillatory piece of μ f(1)rr, and the averaged piece of μ f (2)rr. From this result, we see that the radiative approximation yields the leading-order phase, but fails to determine the first subleading correction. We also see that the approximations (i) – (iii) mentioned above are consistent (so long as the parameters of the ‘geodesic’ source are allowed to vary slowly) at leading order in the two-timescale expansion, but diverge from one another beyond that order. Hence, we see that an adiabatic approximation is generically insufficient to extract parameters from a waveform, since doing so requires a description of the inspiral accurate up to small (i.e., smaller than order-1) errors. But we also see that an adiabatic approximation based on the radiative Green’s function may be an excellent approximation for other purposes, such as detection.

To understand this result, consider the following naive analysis of a quasicircular EMRI — that is, an orbit that would be circular were it not for the action of the self-force, and which is slowly spiraling into the large central body. We write the orbital frequency as (1) ω (0)(E) + (m ∕M )ω 1 (E) + ⋅⋅⋅, where ω (0)(E ) is the frequency as a function of energy on a circular geodesic, and (m ∕M )ω (1)(E) 1 is the correction to this due to the conservative part of the first-order self-force (part of the correction also arises due to oscillatory zeroth-order effects combining with oscillatory first-order effects, but for simplicity we ignore this contribution). Neglecting oscillatory effects, we write the energy in terms only of its slow-time dependence: E = E (0)(&tidle;t) + (m ∕M )E (1)(&tidle;t) + ⋅⋅⋅. The leading-order term E (0) is determined by the dissipative part of first-order self-force, while E (1) is determined by both the dissipative part of the second-order force and a combination of conservative and dissipative parts of the first-order force. Substituting this into the frequency, we arrive at

m [ ] ω = ω(0)(E (0)) + --- ω (11)(E (0)) + ω (21)(E (0),E (1)) + ⋅⋅⋅, (2.38 ) M
where (1) ω2 = E (1)∂ ω(0)∕∂E, in which the partial derivative is evaluated at E = E (0). Integrating this over a radiation-reaction time, we arrive at the orbital phase of Eq. (2.37View Equation). (In a complete description, E (t) will have oscillatory pieces, which are functions of t rather than t&tidle;, and one must know these in order to correctly determine (1) ϕ.) Such a result remains valid even for generic orbits, where, for example, orbital precession due to the conservative force contributes to the analogue of ω (11).

2.6 Physical consequences of the self-force

To be of relevance to gravitational-wave astronomy, the paramount goal of the self-force community remains the computation of waveforms that properly encode the long-term dynamical evolution of an extreme-mass-ratio binary. This requires a fully consistent orbital evolution fed to a wave-generation formalism, and to this day the completion of this program remains as a future challenge. In the meantime, a somewhat less ambitious, though no less compelling, undertaking is that of probing the physical consequences of the self-force on the motion of point particles.

Scalar charge in cosmological spacetimes

The intriguing phenomenon of a scalar charge changing its rest mass because of an interaction with its self-field was studied by Burko, Harte, and Poisson [33Jump To The Next Citation Point] and Haas and Poisson [86Jump To The Next Citation Point] in the simple context of a particle at rest in an expanding universe. The scalar Green’s function could be computed explicitly for a wide class of cosmological spacetimes, and the action of the field on the particle determined without approximations. It is found that for certain cosmological models, the mass decreases and then increases back to its original value. For other models, the mass is restored only to a fraction of its original value. For de Sitter spacetime, the particle radiates all of its rest mass into monopole scalar waves.

Physical consequences of the gravitational self-force

The earliest calculation of a gravitational self-force was performed by Barack and Lousto for the case of a point mass plunging radially into a Schwarzschild black hole [14]. The calculation, however, depended on a specific choice of gauge and did not identify unambiguous physical consequences of the self-force. To obtain such consequences, it is necessary to combine the self-force (computed in whatever gauge) with the metric perturbation (computed in the same gauge) in the calculation of a well-defined observable that could in principle be measured. For example, the conservative pieces of the self-force and metric perturbation can be combined to calculate the shifts in orbital frequencies that originate from the gravitational effects of the small body; an application of such a calculation would be to determine the shift (as measured by frequency) in the innermost stable circular orbit of an extreme-mass-ratio binary, or the shift in the rate of periastron advance for eccentric orbits. Such calculations, however, must exclude all dissipative aspects of the self-force, because these introduce an inherent ambiguity in the determination of orbital frequencies.

A calculation of this kind was recently achieved by Barack and Sago [22, 23], who computed the shift in the innermost stable circular orbit of a Schwarzschild black hole caused by the conservative piece of the gravitational self-force. The shift in orbital radius is gauge dependent (and was reported in the Lorenz gauge by Barack and Sago), but the shift in orbital frequency is measurable and therefore gauge invariant. Their main result – a genuine milestone in self-force computations – is that the fractional shift in frequency is equal to 0.4870m ∕M; the frequency is driven upward by the gravitational self-force. Barack and Sago compare this shift to the ambiguity created by the dissipative piece of the self-force, which was previously investigated by Ori and Thorne [135] and Sundararajan [167]; they find that the conservative shift is very small compared with the dissipative ambiguity. In a follow-up analysis, Barack, Damour, and Sago [11Jump To The Next Citation Point] computed the conservative shift in the rate of periastron advance of slightly eccentric orbits in Schwarzschild spacetime.

Conservative shifts in the innermost stable circular orbit of a Schwarzschild black hole were first obtained in the context of the scalar self-force by Diaz-Rivera et al. [55]; in this case they obtain a fractional shift of 0.0291657q2∕ (mM ), and here also the frequency is driven upward.

Detweiler’s redshift factor

In another effort to extract physical consequences from the gravitational self-force on a particle in circular motion in Schwarzschild spacetime, Detweiler discovered [50Jump To The Next Citation Point] that ut, the time component of the velocity vector in Schwarzschild coordinates, is invariant with respect to a class of gauge transformations that preserve the helical symmetry of the perturbed spacetime. Detweiler further showed that 1 ∕ut is an observable: it is the redshift that a photon suffers when it propagates from the orbiting body to an observer situated at a large distance on the orbital axis. This gauge-invariant quantity can be calculated together with the orbital frequency Ω, which is a second gauge-invariant quantity that can be constructed for circular orbits in Schwarzschild spacetime. Both ut and Ω acquire corrections of fractional order m ∕M from the self-force and the metric perturbation. While the functions ut(r) and Ω (r) are still gauge dependent, because of the dependence on the radial coordinate r, elimination of r from these relations permits the construction of t u (Ω), which is gauge invariant. A plot of t u as a function of Ω therefore contains physically unambiguous information about the gravitational self-force.

The computation of the gauge-invariant relation ut(Ω) opened the door to a detailed comparison between the predictions of the self-force formalism to those of post-Newtonian theory. This was first pursued by Detweiler [50], who compared t u(Ω ) as determined accurately through second post-Newtonian order, to self-force results obtained numerically; he reported full consistency at the expected level of accuracy. This comparison was pushed to the third post-Newtonian order [29, 28, 44, 11]. Agreement is remarkable, and it conveys a rather deep point about the methods of calculation. The computation of ut(Ω ), in the context of both the self-force and post-Newtonian theory, requires regularization of the metric perturbation created by the point mass. In the self-force calculation, removal of the singular field is achieved with the Detweiler–Whiting prescription, while in post-Newtonian theory it is performed with a very different prescription based on dimensional regularization. Each prescription could have returned a different regularized field, and therefore a different expression for ut(Ω). This, remarkably, does not happen; the singular fields are “physically the same” in the self-force and post-Newtonian calculations.

A generalization of Detweiler’s redshift invariant to eccentric orbits was recently proposed and computed by Barack and Sago [24], who report consistency with corresponding post-Newtonian results in the weak-field regime. They also computed the influence of the conservative gravitational self-force on the periastron advance of slightly eccentric orbits, and compared their results with full numerical relativity simulations for modest mass-ratio binaries. Thus, in spite of the unavailability of self-consistent waveforms, it is becoming clear that self-force calculations are already proving to be of value: they inform post-Newtonian calculations and serve as benchmarks for numerical relativity.

  Go to previous page Go up Go to next page