{'fraction_non_alphanumeric': 0.03279660338059761, 'fraction_numerical': 0.011343427060802692, 'mean_word_length': 4.955725190839694, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Self-force: Computational Strategies

Building on substantial foundational progress in understanding the effect of a small body's self-field on its own motion, the past 15 years has seen the emergence of several strategies for explicitly computing self-field corrections to the equations of motion of a small, point-like charge. These approaches broadly fall into three categories: (i) mode-sum regularization, (ii) effective source approaches and (iii) worldline convolution methods. This paper reviews the various approaches and gives details of how each one is implemented in practice, highlighting some of the key features in each case.


Introduction
Compact-object binaries are amongst the most compelling sources of gravitational waves. In particular, the ubiquity of supermassive black holes residing in galactic centres [1] has made the extreme mass ratio regime a prime target for the eLISA mission [2][3][4][5][6]. Meanwhile, the comparable and intermediate mass ratio regimes are an intriguing target for study by the imminent Advanced LIGO detector [7]. In order to maximise the scientific gain realised from gravitational-wave observations, highly accurate models of gravitational-wave sources are essential. For the case of extreme mass ratio inspirals (EMRIs) -binary systems in which a compact, solar mass object inspirals into an approximately million solar mass black hole -the demands of gravitational wave astronomy are particularly stringent; the promise of groundbreaking scientific advances -including precision tests of general relativity in the strong-field regime [8][9][10] and a better census of black hole populations -hinges on our ability to track the phase of their gravitational waveforms throughout the long inspiral, with an accuracy of better than 1 part in 10,000 [2]. This, in turn requires highly accurate, long time models of the orbital motion.
For the past two decades, these demands have stimulated an intense period of EMRI research among the gravitational physics community. Despite of the impressive progress made by numerical relativists towards tackling the two-body problem in general relativity (for reviews see Refs. [11][12][13]), the disparity of length scales characterising the EMRI regime is a significant roadblock for existing numerical relativity techniques. Indeed, to this day EMRIs remain intractable by current numerical relativity methods, and successful approaches have instead tackled the problem perturbatively or through post-Newtonian approximations (see Ref. [14] for a review). This article will focus on the first of these; by treating the smaller object as a perturbation to the larger mass, the so-called "self-force approach" reviewed here has been a resoundingly successful tool for EMRI research.
Within the self-force approach, the smaller mass, µ << M , is assumed to be sufficiently small that it may be used as a perturbative expansion parameter in the background of the larger mass, M . Expanding the Einstein equation in µ, we see the smaller object as an effective point particle generating a perturbation about the background of the larger mass. At zeroth order in µ, the smaller object merely follows a geodesic of the background. At first order in µ, it deviates from this geodesic due to its interaction with its self-field. Viewing this deviation as a force acting on the smaller object, the calculation of this self-force is critical to the accurate modelling of the evolution of the system. For the purposes of producing an accurate waveform for space-based detectors, and for producing accurate intermediate mass ratio inspiral (IMRI) models, it will be necessary to include further effects up to second perturbative order [15,16]. Indeed, recent compelling work [17,18] suggests that IMRIs -and even comparable mass binariesmay be modelled using self-force techniques.
A naïve calculation of the first order perturbation due to a point particle leads to a retarded field which diverges at the location of the particle. The self-force, being the derivative of the field, also diverges at the location of the particle and one obtains equations of motion which are not well-defined and must be regularized. A series of formal derivations of the regularized first order equations of motion (now commonly referred to as the MiSa-TaQuWa equations, named after Mino, Sasaki, Tanaka [19] and Quinn and Wald [20] who first derived them) for a point particle in curved spacetime have been developed [19][20][21][22][23][24][25][26][27][28][29][30], culminating in a rigorous work by Gralla and Wald [31] and Pound [32] in the gravitational case and by Gralla, et al. [33] in the electromagnetic case. This was subsequently extended to second perturbative order by Rosenthal [34][35][36][37], Pound [38][39][40], Gralla [41] and Detweiler [42]. These derivations eliminate the ambiguities associated with the divergent self-field of a point particle and provide a well-defined, finite equation of motion. Building upon this foundational progress, several practical computational strategies have emerged from these formal derivations: • Dissipative self-force approaches: While the full first-order self-force is divergent, it turns out that the dissipative component is finite and requires no regularization. This fact has prompted the development of methods for computing the dissipative component alone, sidestepping the issue of regularization altogether. These dissipative approaches fall into two categories: 1. Flux methods: By measuring the orbit-averaged flux of gravitational waves onto the horizon of the larger black hole and out to infinity, the fact that the field is evaluated far away from the worldline means that no divergent quantities are even encountered. This approach yields the time-averaged 1 dissipative component of the self-force [43][44][45][46][47][48][49][50].
2. Local/instantaneous dissipative self-force: The time averaging element of flux methods can be eliminated by instead computing the local, instantaneous dissipative component of the self-force from the half-advanced-minus-half-retarded field [43,[51][52][53].
Both methods, however, fundamentally rely on neglecting potentially important conservative effects which can significantly alter the orbital phase of the system.
• The mode-sum approach: Introduced in Refs. [54,55], and having since been successfully used in many applications, the approach relies on the decomposition of the retarded field into spherical harmonic modes (which are finite, but not differentiable at the particle), numerically solving for each mode independently and subtracting analytically-derived "regularization parameters", then summing over modes.
• The effective source approach: Proposed in [56] and [57], the approach implements the regularization before solving the wave equation. This has the advantage that all quantities are finite throughout the calculation and one can directly solve a wave equation for the regularized field.
• The worldline convolution approach: First suggested in [58,59], one computes the regularized retarded field as a convolution of the retarded Green function along the past worldline of the particle. Although the approach is the most closely related to the early formal derivations, it is only recently that it has been successfully applied to calculations in black hole spacetimes.
For a comprehensive review of the self-force problem, see Refs. [60][61][62][63]. In this paper, I will review the various approaches and give details of how each one is implemented in practice, highlighting the advantages and disadvantages in each case. This paper follows the conventions of Misner, Thorne and Wheeler [64]; a "mostly positive" metric signature, (−, +, +, +), is used for the spacetime metric, the connection coefficients are defined by Γ λ µν = 1 2 g λσ (g σµ,ν + g σν,µ − g µν,σ ), the Riemann tensor is R α λµν = Γ α λν,µ − Γ α λµ,ν + Γ α σµ Γ σ λν − Γ α σν Γ σ λµ , the Ricci tensor and scalar are R αβ = R µ αµβ and R = R α α , and the Einstein equations are G αβ = R αβ − 1 2 g αβ R = 8πT αβ . Standard geometrized units are used, with c = G = 1. Greek indices are used for four-dimensional spacetime components, symmetrisation of indices is denoted using parenthesis [e.g. (αβ)], anti-symmetrisation is denoted using square brackets (e.g. [αβ]) and indices are excluded from symmetrisation by surrounding them by vertical bars [e.g. (α|β|γ)]. Latin letters starting from i are used for indices summed only over spatial dimensions and capital letters are used to denote the spinorial/tensorial indices appropriate to the field being considered. Either x or x µ are used when referring to a spacetime field point and z(τ ) or z µ (τ ) are used when referring to a point on a worldline parametrised by proper time τ . Finally, a retarded (or source) point is denoted using a prime, i.e. z ′ .

Equations of Motion
The formal equations of motion of a compact object moving in a curved spacetime are now well established up to second perturbative order. Writing the perturbed spacetime in terms of a background plus perturbation, αβ + h αβ , the equations of motion essentially amount to those of an accelerated worldline in the background spacetime, with the acceleration given by a well-defined regular field which is sourced by the worldline. To order µ, this coupled system of equations for the worldline and its self-field are commonly referred to as the MiSaTaQuWa equations and are given (in Lorenz gauge, assuming a Ricci-flat background spacetime) by Here, µ is mass of the object, g α ′ α is the bivector of parallel transport, C αβγδ is the Weyl tensor of the background spacetime, and we use the trace- αβ h, h = h γ γ . One can also consider compact objects possessing other types of charge. For example, a particularly simple case is that of a scalar charge q with mass m and scalar field Φ, in which case the equations of motion are given by 2,3 Here, R is the Ricci scalar of the background spacetime and ξ is the coupling to scalar curvature. Similarly, for an electric charge, e, one obtains equations of motion which are given in Lorenz gauge by 2 In the scalar case, it is important to distinguish between the self-force F α = ∇ α Φ and the self-acceleration, which is given by projecting to self-force orthogonal to the worldline, a α = (g αβ (0) + u α u β )F β . In the electromagnetic and gravitational cases the self-force has no component along the worldline and the two may be used interchangeably. 3 We assume that the mass m is small and ignore its effect on the equations of motion.
where R αβ is the Ricci tensor of the background spacetime and A µ is the vector potential.
The key component in all instances is the identification of the appropriate regularized field on the worldline. Detweiler and Whiting identified a particularly elegant choice for the regularized field, written in terms of the difference between the retarded field and a locally-defined singular field, In addition to giving the physically-correct self-force, the Detweiler-Whiting regular field has the appealing feature of being a solution of the homogeneous field equations in the vicinity of the worldline. Most computational strategies essentially amount to differing ways of representing this singular field 4 and obtaining the regularized field on the worldline.

Numerical regularization strategies
In a numerical implementation, it is essential to avoid the evaluation of divergent quantities. In the case of self-force calculations, both the retarded and singular fields diverge on the worldline so one must avoid evaluating them there. Several strategies for doing so have emerged over the years (see Table 1 for a summary).
One option is to only ever evaluate finite, dissipative quantities (e.g. the retarded field far from the worldline or the half-advanced-minus-halfretarded field on the worldline). This is the basis of the dissipative methods mentioned in the introduction. Since these methods effectively avoid the problem of regularization, they will not be discussed further here; we will return to them in Sec. 4. It is worth noting, however, that the methods typically used by dissipative calculations are essentially the same as those used by mode-sum regularization for computing the retarded field, but without the additional regularization step.
This leaves three regularization strategies which allow the regularized field to be computed on the worldline without encountering numerical divergences: worldline convolution, mode-sum regularization, and the effective source approach.

Worldline convolution
The worldline convolution method relies on a split of the regularized selfforce into an "instantaneous" piece and a history-dependent term. The instantaneous piece is easily calculated from local quantities evaluated at the particle's position, The history-dependent term is much more difficult to calculate, as it is given in terms of a convolution of the derivative of the retarded Green function along the worldline's entire past-history (see Fig. 1), The covariant derivatives here are taken with respect to the first argument of the retarded Green function. Likewise, the regularized self-field can be obtained from a worldline convolution of the retarded Green function itself; similar formulae can also be derived for higher derivatives. The retarded Green function appearing in these equations is a solution of the wave equation with a delta-function source, with boundary conditions such that the solutions correspond to purely outgoing radiation at infinity and no radiation emerging from the horizon. The regularization (i.e. subtraction of the Detweiler-Whiting singular field) is formally achieved by the limiting procedure in the upper limit of integration, i.e. by cutting off the integration at τ − , slightly before the Timelike worldline of the particle Current location of the particle -z(!)

Integral in quasilocal region
Integral outside quasilocal region Boundary of causal domain where Hadamard form is valid Figure 1: Schematic representation of the worldline convolution method. The equations shown are for the case of the self-force on a scalar charge, but are representative of similar equations in the electromagnetic and gravitational cases. Reproduced from Ref. [69]. coincidence point z = x. In practice, this is done by taking the integration all the way up to coincidence, but excluding the direct contribution to the Green function at coincidence (i.e. the term proportional to δ(σ) in Hadamard form for the Green function). Although the retarded Green function also diverges at certain other points along the past worldline (in particular at null-geodesic intersections where the particle sees null rays emitted from its past), it turns out that these are all integrable singularities (of the form 1/σ and δ(σ)) and the integral may be accurately evaluated to give a finite and accurate value for the regularized self-force.
Despite having been proposed in the early days of EMRI self-force calculations [58], the worldline convolution method was largely ignored for a long time by numerical implementations (the notable exception being investigative studies [59,127] which did not complete a full calculation of the self-force). A likely reason is that the method relies on knowledge of the retarded Green function for all points on the particle's past worldline. While methods have existed for decades for computing portions of the re-tarded Green function, it turns out to be very difficult to obtain it accurately everywhere that it is needed for the worldline convolution.
Thankfully, recent progress has led to two practical computational strategies, both of which have been successfully applied to compute the self-force in black hole spacetimes. The first of these is based on a frequency-domain decomposition of the Green function, and builds on the rich history of black hole perturbation theory developed over the past several decades. The second, a time-domain approach, has been a very recent development and shows a great deal of promise for the future.
Independently of whether a frequency-domain or a time-domain scheme is used, a common problem with both is that they fail at early times when source and field points are close together; in the frequency-domain case the convergence is poor at early times, while in the time-domain case features from the numerical approximation pollute the data at early times. A relatively straightforward solution which has been successfully applied in both scenarios is to only rely on their results at late times and to supplement them at early times with a quasilocal Taylor series expansion. Provided a sufficiently early time can be chosen where both the distant past and quasilocal calculations converge, this yields a global approximation for the Green function which is sufficient for producing an accurate result for the self-force. To date this has been shown to be possible in the case of a scalar charge in Nariai [128], Schwarzschild [69,70] and Kerr [129] spacetimes.
Given an approximation for the Green function valid throughout the past worldline, it is then trivial to numerically integrate Eq. (6) to obtain the selfforce (see Fig. 1). As mentioned previously, the divergences in the retarded Green function at null-geodesic intersections on the past worldline are all integrable singularities and do not pose a significant obstacle to accurate numerical evaluation of the integrals.

Quasilocal expansion
In order to obtain an approximation to the retarded Green function which is valid at early times, it is convenient to start with the Hadamard form for the Green function (in Lorenz gauge) where Θ − is analogous to the Heaviside step-function, being 1 when x ′ is in the causal past of x, and 0 otherwise, δ(σ) is the covariant form of the Dirac delta function, U AB ′ and V AB ′ are symmetric bi-spinors/tensors and are regular for x ′ → x. The bi-scalar σ (x, x ′ ) is the Synge world function, which is equal to one half of the squared geodesic distance between x and x ′ . In particular, σ(x, x) = 0 and σ(x, z) < 0 when x and z are timelike separated. Because of the limiting procedure in the history integral, Eq. (6), only the term involving V AB ′ is non-zero, and at all required points along the worldline Θ − = 1 = Θ(−σ). The problem of determining the retarded Green function at early times therefore reduces to finding an approximation for V AB ′ (x, x ′ ) which is valid for x and x ′ close together. Several methods have been developed for computing approximations to V AB ′ . Fundamentally, they rely on either the use of a series expansion, or on the use of numerically evolved transport equations (ordinary differential equations defined along a worldline). The series expansion approach has been the most fruitful to date with results including: leading-order coordinate expansions in Schwarzschild and Kerr spacetimes for scalar [68,73] and gravitational cases [101], high-order coordinate expansions in spherically symmetric spacetimes (including Schwarzschild) [67], formal covariant expansions in generic spacetimes [130][131][132][133], and moderately high-order coordinate expansions in Schwarzschild [82] and Kerr [90] spacetimes. The only numerical calculation I am aware of was done in [133] for generic spacetimes (with an example application in Schwarzschild spacetime).
The series expansion method produces an expression for V (x, x ′ ) as a power series in the coordinate distance between x and x ′ . For example, for the scalar case in Schwarzschild spacetime it takes the form where γ is the angular separation of the points and the v ijk are analytic functions of r and M . It is straightforward to take partial derivatives of these expressions at either spacetime point to obtain the derivative of the Green function. Although this series on its own may be sufficient for use in the quasilocal component of a worldline convolution, it turns out that some simple tricks allow for a vast improvement in accuracy. It turns out that, because V (x, x ′ ) diverges at the edge of the normal neighbourhood, the series approximation benefits significantly from Padé resummation which incorporates information about the form of the divergence [67]. One minor caveat is that since Padé re-summation is only well defined for series expansions in a single variable, it is necessary to first expand r ′ and γ in a Taylor series in t − t ′ , using the equations of motion to determine the higher derivatives appearing in the series coefficients. Then, with V (x, x ′ ) written as a power series in t − t ′ alone, a standard diagonal Padé approximant provides an accurate representation of the Green function in the quasi-local region.

Frequency domain methods
Frequency domain methods for computing the retarded Green function rely crucially on the separability of the wave equation. In the scalar, Schwarzschild case that represents the current state-of-the art 5 [69] (also see [128] for a related calculation in Nariai spacetime), this can be achieved by writing the Green function as a sum of spherical harmonic and Fourier modes Here, ǫ > 0 is a formal parameter to ensure the correct boundary conditions are satisfied for a retarded Green function. Substituting this into the wave equation, Eq. (7a), one obtains an independent set of ordinary differential equations forĝ ℓ (r, r ′ ; ω), one equation for each ℓ and ω, with Here, r * ≡ r + 2M ln r 2M − 1 is the radial tortoise coordinate. Given two linearly-independent solutions, p(r; ω) and q(r; ω), of the homogeneous version of (11),ĝ ℓ (r, r ′ ; ω) is given bŷ where r > ≡ max(r, r ′ ), r < ≡ min(r, r ′ ) and W (p, q) is the Wronskian. One is then faced with the integral over frequencies in the inverse Fourier transform appearing in (10). This may be achieved through straightforward integration  along the real-ω axis (see Fig. 2); the only caveat is that the formally-infinite integral over frequencies should be cut off at some finite maximum frequency using a smooth window function, for example Alternatively, as was done in [69], the integration contour can be deformed into the complex frequency plane following a proposal by Leaver [144,145]. In Schwarzschild spacetime,ĝ ℓ (r, r ′ ; ω) has simple poles in the lower semi-plane at the quasinormal mode frequencies, along with a branch cut starting at ω = 0 and continuing along the negative imaginary axis (see Fig. 3). The deformed frequency integral is therefore given by an integral over a high-frequency arc, an integral around a branch cut, and an integral over the residues at the poles 6 .
Re Ω Im Ω HF BC QNMs Figure 3: Contour deformation in the complex-frequency plane. The residue theorem of complex analysis allows one to re-express the integral over (just above) the real line of the Fourier modesĝ ℓ (r, r ′ ; ω) as an integral over a high-frequency arc plus an integral around a branch cut of the Fourier modes plus an integral over the residues at the poles of the Fourier modes. Reproduced from Ref. [69].
The high-frequency arcs can be disregarded as they are likely to only contribute at very early times, where the quasilocal approximation can be used. There are well-established methods for accurately computing the location of the poles (quasinormal mode frequencies) and the residues at the poles. The biggest difficulty in the frequency-domain approach is in evaluating the branch cut integral, about which very little was known until recently (other than asymptotic approximations for e.g. large radius or late times). Substantial recent progress has established methods for calculating this branch cut contribution [146]; these methods were used in [69] to compute the self-force in Schwarzschild spacetime.
The final step in the frequency domain approach is to sum over spherical harmonic modes to produce the full Green function. Here, the distributional parts of the full Green function can cause poor convergence in the modesum. Fortunately, there is a straightforward solution to this problem: by smoothly cutting off the mode sum at large ℓ, one obtains a mollified retarded Green which is appropriate for use in a self-force calculation, and whose sum over ℓ converges. Empirically, it has been found that choosing ℓ cut ≈ ℓ max /5 gives good results.

Time domain approaches
The frequency domain approach to computing the retarded Green function has several shortcomings. It relies on relatively difficult technical calculations and has poor convergence properties at early times and in scenarios where the worldline is not well-represented by a discrete spectrum of frequencies. Recent work has shown that the Green function may also be accurately computed (at least for the purposes of self-force calculations) in the time domain using straightforward numerical evolutions. A time domain calculation sidesteps issues related to a wide frequency spectrum and appears to exhibit much better convergence properties at early times. There are two closely-related proposals for computing the Green function in the time domain. In [71], Eq. (7a) was numerically solved as an initial value problem, with the delta-function source approximated by a narrow Gaussian. A reformulation of this approach as a homogeneous problem with Gaussian initial data was subsequently given in [70]. It was found that these "numerical Gaussian" approaches are able to approximate the retarded Green function remarkably well in a large region of the space required by worldline convolutions. The size of the Gaussian limits the scale of the smallest features which can be resolved, but it turns out that this is not significantly detrimental to a self-force calculation.
The only regimes where the numerical Gaussian approach is not well suited to computing the retarded Green function are at very early and very late times. At early times the direct δ(σ(x, x)) term (which must be excluded from the worldline convolution of retarded Green function) is smeared out and is difficult to isolate from the rest of the Green function. There is no conceptual difficulty at late times, but the reality of a numerical evolution is that it can only be run for a finite time. Fortunately, both of these issues are easily overcome; the former by using the quasilocal expansion at early times, and the latter by using a late-time expansion of the branch cut integral at late times (see Fig. 4).  With an eccentric orbit of period ≈ 317M , we see that a good approximation to the self-field is obtained by including the contributions from less than one full orbit. Note that the formal divergences of the Green function are integrable and do not cause significant numerical difficulty in computing the self-force.
In the time domain approach, each numerical time domain evolution gives the Green function G(x 0 , x ′ ) for a single base point x 0 , and for all source points, x. As a result, a single numerical calculation can only be used to compute the self-force at a single spacetime point, but it can be computed for any past-worldline ending up at that point. This is in contrast to other self-force methods, where a single orbit is considered at a time, but a single calculation gives the self-force at all points on the orbit. The problem of efficiently spanning the parameter space of base points x 0 is an ideal application of reduced order methods [129,147].

Mode-sum regularization
The mode-sum regularization scheme, first proposed by Barack and Ori [54], has proven highly successful as a computational self-force strategy, having been applied to the computation of the the self force for a variety of configurations in Schwarzschild [74-81, 83, 102, 107-110, 113, 114, 148, 149] and Kerr [89,91,92,124,150] spacetimes. The success of the method hinges on the fact that the first order retarded field diverges in a way which can be effectively smoothed out by a spectral decomposition in the angular directions. More specifically, the divergence appears as an odd power of 1/s, where s 2 = (g αβ + u α u β )σ ;α σ ;β is an appropriate measure of distance from the worldline. The result is that the 1/s divergence of the field near the worldline turns into an infinite sum of modes, each of which are individually finite (but possibly discontinuous) on the worldline. The divergence of the field then manifests itself through the failure of the (infinite) sum over modes to converge. Conveniently, this odd-in-1/s property of the retarded field also holds for its derivatives, both the first derivatives required for the self-force and higher derivatives which are useful for computing higher-multipole gauge invariants [104,105]. As a result, an arbitrary number of derivatives of the first order retarded field is represented by an infinite sum of modes, each of which are individually finite (but possibly discontinuous) on the worldline. The trade-off is that as the number of derivatives increases the divergence of the sum over modes becomes increasingly strong.
Since the individual modes are finite, numerical calculations of the retarded field and its derivatives can be done in the reduced (t, r) space without encountering any numerical divergences. The remaining piece of the problem is a method for rendering the divergent sum over modes finite. An analytic decomposition of the angular dependence of the Detweiler-Whiting singular field yields "regularization parameters" which may be subtracted mode-by-mode from the numerical retarded field values. Provided all parts of the Detweiler-Whiting singular field which don't vanish on the worldline are included in the calculation, the regularized sum over modes is convergent on the worldline and one never encounters any numerical divergences.
It is important to note that the use of the Detweiler-Whiting singular field is not merely a convenience; it also provides well-motivated physical grounds to justify the mode-sum regularization approach. Other ad-hoc approaches based on identifying the asymptotically divergent contributions to the mode-sum often produce the correct regularization parameters, but their use is difficult to justify rigorously and can easily lead to incorrect results. Crucially, there is no way of knowing whether such an ad-hoc regularization procedure is producing a physically correct result, or if important contributions are being overlooked. They should therefore not be relied on without extreme care. Unfortunately, despite its success in first order calculations, mode-sum regularization alone is not an effective tool for computing the second order self-force. The reason for this is straightforward: the second order retarded field contains divergent terms which appears in the form of even powers of 1/s. Intuitively, this arises from the fact that the second order field contains terms which depend quadratically on the first order field. Unlike the oddin-1/s case, the angular decomposition of 1/s 2 leads to individual modes which diverge logarithmically as the worldline is approached. Fortunately, all is not lost for the mode-sum method as a computational tool at second order. Provided the leading order logarithmic divergence is subtracted by some other means (for example, using the effective source approach), regularization parameters may be used to accelerate the rate of convergence of the mode sum. Such a hybrid scheme complements the generality of the effective source approach with the computational efficiency of mode-sum regularization.

Regularization parameters
The computation of regularization parameters has been addressed by a series of calculations stemming from the original Barack-Ori derivation. Barack and Ori's original work gave the first two self-force regularization parameters for scalar, electromagnetic and Lorenz-gauge gravitational cases in both Schwarzschild and Kerr spacetimes. These are sufficient for computing the regularized self-force, but they yield mode-sums which have relatively poor quadratic convergence with the number of modes included. Since the Detweiler-Whiting regularized field is a smooth, homogeneous function in the vicinity of the worldline, its mode-sum representation converges expo-nentially. This spectral convergence is spoiled by the fact that only the portion of the Detweiler-Whiting singular field which does not vanish on the worldline is subtracted by the leading order regularization parameters. The mode decomposition effectively contains information about the extension of the field off the worldline to the entire two-sphere, so the regularized field contains residual pieces of the Detweiler-Whiting singular field off the worldline, but on a two-sphere of the same radius. By deriving higher-order regularization parameters one can subtract these residual pieces order-byorder, leaving a mode-sum which is more and more rapidly convergent (see Fig. 5). The derivation of these higher-order parameters is closely related to the computation of the quasilocal expansion of the Green function and has been addressed in a series of papers: the first higher order parameter was given in [76] for the case of a scalar charge on a circular orbit in Schwarzschild spacetime, and for eccentric geodesic orbits in [80]. This was subsequently extended by several further orders for equatorial geodesic motion in the scalar, electromagnetic and gravitational cases in both Schwarzschild [82] and Kerr [90] spacetimes. Recent work has also produced parameters for accelerated worldlines [72,93]. Figure 5: Spherical harmonic modes of the self-force for a point scalar charge on a circular orbit of radius r 0 = 6M in Schwarzschild spacetime. By subtracting analytically determined high-order "regularization parameters", the sum over modes is rendered more and more rapidly convergent. Each regularization parameter, F r , behaves asymptotically for large-ℓ as 1/ℓ n .
Within these derivations of high-order regularization parameters, a subtle ambiguity appears through the elevation of the four-velocity from a quantity defined on a worldline to a quantity defined everywhere on the twosphere. A natural, covariant choice is to define this off-worldline extension through parallel transport, u α ′ = g α α ′ u α . However, in practical calculations it is often convenient to make a coordinate choice. For example, a common choice is to define the extension in terms of "constant coordinate components", i.e. to define u α ′ such that its components in some coordinate system have a constant value everywhere on the two-sphere. This is a perfectly valid choice, as it any other choice where u α ′ agrees with the actual four-velocity when evaluated on the worldline. The only caveat is that the regularization parameters beyond the leading two orders change depending on the particular choice of extension. As a result, in order to use higherorder regularization parameters it is essential that a compatible choice of off-worldline extension is used in the retarded field calculation.

Choice of basis and gauge
There are two other factors which must be considered in the mode-sum scheme: the choice of angular basis functions for the spectral decomposition and the choice of gauge in the electromagnetic and gravitational cases.
The choice of basis functions is typically motivated by their ability to produce separability of the retarded field equations; an appropriate choice of basis functions yields an independent set equations in (t, r) space for each individual mode. For the case of a scalar charge in Schwarzschild spacetime, the appropriate basis functions are the standard spherical harmonics. In the Kerr case the spheroidal harmonics are required for separability. For electromagnetic and gravitational cases, there are several choices. For Schwarzschild spacetime one can choose between a tensor harmonic basis and a basis of spin-weighted spherical harmonics. In the Kerr case only the spinweighted spheroidal harmonics are known to yield separability [134,135].
This choice of angular basis also affects the regularization parameters. The parameters for scalar spherical harmonic modes, tensor harmonic modes, spin-weighted spherical harmonic modes and spheroidal harmonic modes are all potentially different. However, it is always possible project the tensor, spin-weighted, or spheroidal harmonic modes onto the scalar spherical harmonic basis. In fact, since the regularization parameters are most easily obtained for the scalar spherical harmonic basis, calculations of the retarded field are typically done using the computationally most convenient basis and the result is then projected onto the scalar spherical harmonics before the regularization and sum-over-modes steps are done [113].
A more difficult issue in the electromagnetic and gravitational cases is the choice of gauge. The Detweiler-Whiting singular field is defined in Lorenz gauge (since it is derived from a Lorenz gauge Green function), but numerical calculations of the retarded field are more easily done in either radiation or Regge-Wheeler gauge. The existence of tensor spherical harmonics makes a Lorenz-gauge calculation possible, if somewhat cumbersome in the Schwarzschild case [111,118,120]. One obtains a coupled set of 10 equations for the metric perturbation, but there is no coupling between different different tensor-harmonic modes. Unfortunately, this does not hold for the Kerr case as there are no known tensor spheroidal harmonics. Rather than trying to derive tensor spheroidal harmonics, a better approach is to work with the relatively straightforward Teukolsky equation in radiation gauge [124]. The difficulty then is in identifying the appropriate regularization parameters, particularly since the gauge transformation from Lorenz gauge is itself often singular. This remained an open problem until recently, when an understanding of how to apply the mode-sum regularization scheme in radiation gauge was finally established [119] and applied in a self-force calculation [108].

Mode-sum regularization in the frequency domain
The mode-sum scheme provides a method for computing the regularized self-force using numerical data for the modes of the retarded field in the reduced (t, r) space. There are, however, several possibilities for computing this numerical data. One option relies on a further decomposition of the time dependence into Fourier-frequency modes, in an analogous way to the frequency-domain Green function described in Sec. 3.1.2 above. Using this spin-weighted spheroidal harmonic basis, this leads to a radial equation for the Teukolsky function, with the potential given by Here, ∆ ≡ r 2 − 2M r + a 2 , K ≡ (r 2 + a 2 )ω − am, λ is an eigenvalue of the spheroidal equation, and a and M are the Kerr spin and mass parameters. For a = 0, λ = (ℓ − s)(ℓ + s + 1) and this reduces to the Schwarzschild wave equation decomposed into spin-weighted spherical harmonics (in radiation gauge for the gravitational case), while for s = 0 it reduces to the scalar wave equation.
A decomposition of the point-particle source term yields a source involving δ(r − r p ) (and in some cases its derivative), where r p is the radial location of the particle. In practice, the frequency domain mode-sum approach proceeds in the same way as for the Green function: two independent solutions of the homogeneous radial equation are obtained and are matched at the particle's location 7 . Then, the distributional sources do not introduce any numerical difficulty as they simply appear as jumps when matching the homogeneous inner and outer solutions. The solutions of the Teukolsky equation for a given (ℓ, m, ω) can be obtained either through straightforward numerical integration of the radial ordinary differential equation (typically with some modifications to improve numerical accuracy [151,152]) or as an approximation in the form of an infinite convergent series of hypergeometric functions. The latter method is based an idea originally developed by [153] and now commonly referred to as the "MST" method, after Mano, Suzuki, and Takasugi [154] who reformulated it into its current form. It provides an efficient and highly accurate method for computing solutions of the radial equation. For example, recent results have used it to compute solutions accurate to several hundred decimal places, allowing the solutions to be used to determine previously unknown high-order post-Newtonian parameters [155][156][157][158][159]. For a comprehensive review of the MST approach, see the living review by Sasaki and Tagoshi [160] and references therein.
The frequency-domain approach is particularly appropriate in scenarios where the worldline is well represented by a discrete spectrum involving a small number of frequencies. In such cases, mode-sum regularization is by far the optimal choice and is unparalleled in its accuracy and computational efficiency [161]. The prototypical example is a circular orbit, in which only a single frequency is present. In the case of eccentric equatorial orbits (and inclined circular orbits in the Kerr case), there are two fundamental frequencies, and also an infinite number of higher harmonics produced from combinations of the fundamental frequencies. For mildly eccentric orbits this does not cause a great deal of difficulty. However, for more eccentric orbits (with eccentricities e 0.5) an increasingly large number of frequencies must be included and the competitive advantage of the frequency domain approach is lost [161]. Even worse, generic geodesic orbits in Kerr spacetime have three fundamental frequencies and the computational difficulty is so high that a calculation has yet to be attempted. Similarly, unbound orbits and other cases such as radial infall are not well suited to frequency domain methods. Apart from these deficiencies, the frequency domain mode-sum approach has been highly successful for producing results.

Mode-sum regularization in the time domain
The mode-sum scheme may also be applied in the time domain by skipping the Fourier decomposition step and instead solving a set of 1 + 1D partial differential equations in (t, r) space. The main difficulty then is in appropriately handling the distributional source term which has the form δ(r−r p (t)). One solution, used in [49,79,98,113,149,162], is to use a discretised representation of a delta function and to construct the computational grid such that the worldline only ever passes between grid points.
An alternative approach is to reformulate the problem in an analogous way to the frequency-domain method. By splitting the computational domain into two domains -one either side of the particle -the delta-function source can be reformulated in terms of a jump in the fields and their derivatives at the interface of the two domains. This method is well suited to highly-accurate spectral numerical methods as all of the numerically evolved fields are smooth functions. It was implemented using discontinuous-Galerkin methods in [106,115] and using Chebyshev pseudospectral methods in [78,81]. Eccentric orbits present a small additional complexity in this case as the particle must always lie on the domain interface. This is easily achieved by introducing a time-dependent mapping between the computational and physical coordinates of the system [81,106].

Limitations of the mode-sum regularization scheme
Despite its resounding success to date, the mode-sum regularization scheme has some unfortunate disadvantages which make it ill-suited as a generalpurpose method for self-force calculations: 1. Its application in the Kerr case is only straightforward in the frequency domain, since the field equations are not separable in the time domain 8 2. It relies on the use of Lorenz gauge for regularization in the gravitational case. This is not a major issue in the Schwarzschild case since the tensor spherical harmonics may be used to decouple the Lorenz gauge field equations in the angular directions (leaving 10 coupled 1 + 1D equations for each ℓ, m mode). However, there are no know tensor spheroidal harmonics which would be required for the Kerr case, and even if there were they would likely not be applicable in the time domain (again, it is conceivable that a coupled system of equations involving tensor spherical harmonics could be used, but the coupling would result in considerable complexity).
3. It is not applicable beyond first perturbative order, since the modes of the second order perturbation diverge logarithmically near the worldline.
The first two issues are not necessarily showstoppers. There have been several attempts at 1 + 1D time-domain implementations using coupled spherical-harmonic modes in the Kerr case [163,164], and recent progress on reformulating mode-sum regularization for radiation gauge has clarified the regularization issue [119]. The third point, however, appears insurmountable. For these reasons, among others, the effective source method, described in the next section, was developed.

Effective source approach
Proposed in 2007 as a solution to the shortcomings of mode-sum regularization, the effective source approach 9 provides an alternative method for handling the divergence of the retarded field. Rather than first computing the 8 It may still be possible to use mode-sum regularization for the Kerr case in the time domain by decomposing the field equations into spherical harmonics and evolving the resulting infinitely coupled set of 1 + 1D partial differential equations in a similar manner to the Schwarzschild case. 9 Note that the effective source proposed by Lousto and Nakano [65] is similar in spirit, but differs in that it is not derived from the Detweiler-Whiting singular field. retarded field and then subtracting the singular piece as a post-processing step, one can instead work directly with an equation for the regular field. This idea -independently proposed by Barack and Golbourn [56] and by Vega and Detweiler [57] -has the distinct advantage of involving only regular quantities from the outset, making it applicable in a wider variety of scenarios than the mode-sum scheme. In particular, since it does not rely on a mode decomposition of the retarded field, it can be used by any numerical prescription for solving the retarded field equations, whether in the frequency domain (where the method is really just a generalisation of modesum regularization) or in the time domain as a 1 + 1D, 2 + 1D or even 3 + 1D problem.
The basic idea is to use the split of the retarded field into regular and singular pieces, Eq. (4), to rewrite the field equations, Eqs. (1a), (2a) and (3a) as equations for the Detweiler-Whiting regular field, If the singular field used in the subtraction is exactly the Detweiler-Whiting singular field, then the two terms on the right hand side of this equation cancel and the regularized field would be a homogeneous solution of the wave equation. Unfortunately, one typically does not have access to an exact expression for the singular field. Indeed, the Detweiler-Whiting singular field is only defined through a Hadamard parametrix which is not even defined globally. Instead, the best one can typically do is a local expansion which is valid only in the vicinity of the worldline. Borrowing the language of Barack and Golbourn, we refer to an approximation to the singular field as a "puncture" field, Φ S ≈ Φ P , A S α ≈ A P α ,h S αβ ≈h P αβ . Then, the corresponding approximate regular field -referred to as the "residual field" -is no longer a solution of the homogeneous equation, but instead is a solution of the sourced equation with an effective source (see Fig. 6) which is defined to be the right-hand side of Eq. (19) with the puncture field substituted for the singular field Note that the presence of a distributional component of the source on Figure 6: The effective source for a scalar particle on a circular orbit in Schwarzschild spacetime. The source is generically non-smooth in the vicinity of the worldline (left), but the smoothness can be improved by incorporating higher-order parts of the Detweiler-Whiting singular field into the source. the worldline is merely a formal prescription; in practice the puncture field is chosen so that it exactly cancels this distributional component on the worldline. This effective source is then finite everywhere, but has limited differentiability on the worldline. This makes it well-suited to numerical implementations since no divergent quantities are ever encountered. The only numerical difficulty arises from the non-smoothness of the source in the vicinity of the worldline (see Fig. 6), which leads to numerical noise in the computed self force. The noise can be reduced by making the source smoother using higher-order parts of the Detweiler-Whiting singular field. As shown in Fig. 7, at the same numerical resolution a higher order source (C 2 in this case) eliminates the vast majority of the numerical noise that is present when using a lower order source (C 0 in the case in the figure). The cost of this improved accuracy is a source which is considerably more complicated, and costly to compute in a numerical code. Figure 7: Radial self-force for a scalar particle on an eccentric orbit in Schwarzschild spacetime, computed with a 3 + 1D implementation of the effective source scheme. Here, the independent variable χ is a "relativistic anomaly parameter" defined through r = pM/(1 + e cos χ), with p the semilatus rectum and e the eccentricity. The high-frequency errors using a continuous but non-differentiable source (blue) are dramatically decreased by using a twice differentiable source (orange) obtained from a higher-order approximation to the Detweiler-Whiting singular field. For reference, a highly accurate value computed using frequency-domain mode-sum regularization is also included (dashed, black). Figure based on version presented in Ref. [165].
An additional level of complexity arises from the fact that the puncture field is defined only in the vicinity of the worldline. To avoid ambiguities in its definition far from the worldline, one must ensure that the puncture field goes to zero there. This is most easily achieved by multiplying it by a window function, W, with properties such that it only modifies the puncture field in a way that its local expansion about the worldline is preserved to some chosen order. In a first-order calculation of the self-force, it suffices to choose W such that W(x p ) = 1, W ′ (x p ) = 0, W ′′ (x p ) = 0 and W = 0 far away from the worldline. The residual field (see Fig. 8) then obeys The ℓ = 1, m = 1 spherical-harmonic mode of the residual field, Φ res , for a scalar point particle on a circular orbit in Schwarzschild spacetime. This was produced in [86] using a frequency-domain effective source approach.
and has the properties As the residual field coincides with the retarded field far from the particle we can use the usual retarded field boundary conditions when solving Eqs. (21a), (21b) and (21c). The details of a numerical implementation of the effective source approach are then much the same as for the Green function and mode-sum regularization schemes. One can use either a frequency domain or a time domain method for solving the field equations, the key differences now being that there is no restriction to 1 + 1D, and that the effective source must be included as a source term. A more thorough review of the effective source approach -including a detailed description of methods for obtaining explicit expressions for the effective source -can be found in Refs. [126,166].
The effective source approach has been successfully applied in the frequency domain [86], and in the time domain in 1+1D [57], 2+1D [84,96,122] and 3 + 1D [85,87,88] contexts. It has also been formulated -but not yet implemented -in the gravitational case to second order in the mass ratio [40,41]; for the second order problem it is currently the only viable computational strategy. Since the effects of the second-order metric perturbation will be very small, being suppressed by two orders of the mass ratio relative to test body effects it is likely that a highly accurate numerical scheme will be required, suggesting a frequency domain treatment of the problem where one encounters ordinary differential equations (ODEs) which are relatively easy to compute numerically to high accuracy.

Evolution schemes
The calculation of the self-force is only the first stage in the production of an inspiral model. Another critical component is a scheme for evolving the orbit using this self-force information. Various approaches have been proposed, each of which brings with it its own advantages and disadvantages. Approaches which make approximations in the self-force used to drive the inspiral can give substantial decreases in the computational cost of an inspiral simulation, but come at the cost of ignoring potentially relevant physical effects.

Dissipation driven inspirals
A straightforward model for the inspiral can be obtained from energy balance considerations. Using a flux calculation -in which fluxes of energy and angular momentum are obtained by evaluating the point-particle retarded field near the horizon and out at infinity -the entire issue of regularization is avoided and one obtains an approximation to the contributions to the inspiral coming from dissipative self-force effects. However, this is inadequate for capturing all relevant effects from the first-order self-force. Being a dissipative approximation, it completely misses all conservative corrections to the motion. Furthermore, there is no well-defined way of associating fluxes far from the worldline with an instantaneous local self-force on the worldline; a flux model can only be used to drive an inspiral in a time-averaged sense, ignoring some potentially important dissipative contributions. Nevertheless, the relative straightforwardness and computational efficiency of their implementation have made flux models a compelling approach to assessing qualitative features of self-force driven inspirals. These "kludge" models have been used to produce kludge waveforms for EMRI systems which capture at least some of the qualitative physical effects [44,45,48,49,162].  To improve on the flux-based dissipative model, one can instead make use of a half-retarded minus half-advanced scheme [43,[51][52][53] to compute a local, instantaneous dissipative self-force (see Fig. 9). This approach captures all dissipative effects responsible for driving the inspiral, but still neglects small but potentially important conservative corrections to the orbital phase of the system [121].

Osculating self-forced geodesics
In order to account for both conservative and dissipative effects from the first order self-force, it is essential to build an inspiral model around a full calculation of the local self-force. This, however, still leaves some flexibility in the choice of method for computing this local self-force. One option, based on the osculating geodesics framework [167] is particularly compelling as it allows for dramatic improvements in computational efficiency by separat-ing the coupled problem of simultaneously solving for the orbit and for the regularized retarded field into two independent, largely uncoupled computational problems. The basic idea is to perturbatively expand the worldline about a geodesic of the background spacetime, z(τ ) = z 0 (τ )+· · · . Then, the self-force at first order is a functional not of the evolving orbit z(τ ), but of the geodesic orbit, z 0 (τ ), which is instantaneously tangent to the worldline. This expansion is valid in the adiabatic regime, where the orbit is evolving slowly and the difference between the geodesic and evolving orbits is small; the error introduced by the approximation appears at the same order in the equations of motion as a second-order perturbative correction to the field equations.
The improvements in computational efficiency brought about by the osculating geodesics approximation are dramatic. The problem of selfconsistently computing the regularized self-force coupled to an arbitrarily evolving orbit is reduced to the much simpler problem of determining the self-force for geodesic worldlines. Since orbits of black hole spacetimes are parametrised by at most three conserved quantities (energy and angular momentum in the Schwarzschild case, supplemented by the Carter constant in the Kerr case), it is computationally tractable to span the entire parameter space of moderately-eccentric geodesic orbits. Even better, the highly-accurate frequency domain mode-sum method can be used because bound geodesic orbits are efficiently represented by a small frequency spectrum. In Ref. [121] this approach was explored over an entire radiationreaction timescale for moderately-eccentric inspirals in the Schwarzschild spacetime 10 . This was achieved by tabulating the values for the self-force in a relevant portion of the energy-angular momentum phase space of geodesic orbits, and using an interpolated model of the tabulated results to drive an orbital inspiral.

Self-consistent evolution
Unfortunately, the adiabatic approximation responsible for the dramatic improvements in computational efficiency brought by the osculating geodesics framework also introduces errors in the equation of motion at second perturbative order, making the method inadequate for the purposes of precision EMRI astrophysics. One possible solution, implemented in [88] for the scalar field case, is to avoid expanding out the worldline and instead evolve the self-consistent coupled system of equations, Eqs. (1a) and (1b); (2a), (2b) and (2c); or (3a) and (3b). This is a computationally much more difficult problem as the field equations must be solved in the time domain and one cannot rely on an efficient off-line tabulation of self-force values. This selfconsistent evolution can in principal be implemented using a time-domain mode-sum scheme, but in practice implementations have instead used the effective source approach since that will be necessary at second perturbative order (see Sec. 3.2 for an explanation of this issue).
While a self-consistent evolution incorporates all effects contributing to the first order self-force, and only neglects second order contributions from the second order field equations, this comes at the cost of computational efficiency. Whereas an osculating geodesics framework can evolve a large number (∼ 10, 000) of orbits with ease once the initial off-line tabulation phase is complete [121], the self-consistent scheme requires a long calculation of the solutions of the first-order field equations for each new orbit. In reality, this limits the practicality of the scheme to tens of orbits using current methods. The self-consistent evolution scheme is therefore most valuable as a benchmark against which other, less accurate but more efficient methods can be validated. In fact, comparisons between the self-consistent and osculating geodesics scheme for the scalar case indicate that the osculating geodesics scheme performs remarkably well [169].

Two-timescale expansions
The osculating geodesics scheme is fast, but inaccurate, while the selfconsistent evolution scheme is accurate, but slow. This begs the question of whether there is a middle ground which incorporates most of the accuracy of a self-consistent evolution while maintaining the computational efficiency of the osculating geodesics scheme. One promising possibility is that the use of a two-timescale expansion could give just such a scheme. In the two-timescale scheme, rather than using an expansion about a background geodesic (which is valid over short timescales characterised by the orbital motion), one introduces an additional radiation reaction timescale into the problem and incorporates the relevant effects over this radiation reaction timescale.
The relevant two-timescale expansion of the worldline equations of motion was completed in Ref. [52], and follow-up work has improved our understanding of important resonance effects for inspiral orbits [53,[170][171][172][173][174][175][176][177][178][179]. In order to consistently incorporate a two-timescale expansion into an orbital evolution scheme, it will also be necessary to have a two-time expansion of the field equations. Promising progress towards such an expansion was recently reported in [180], indicating that it is likely that a self-consistent orbital evolution using a two-timescale expansion is indeed feasible.

Discussion
The techniques described in this review represent the current state of the art of self-force calculations. The three primary approaches: worldline convolutions, mode-sum regularization and the effective source approach can be considered complimentary, with each having regimes where they are most appropriate: • The mode-sum scheme gives unparalleled accuracy (particularly in the frequency domain) for orbits which have a small frequency spectrum.
• The worldline convolution method provides valuable physical insight and can be easily applied to arbitrary orbital configurations, including those inaccessible by other means.
• The effective source approach can be used in arbitrary spacetimes without relying on any symmetries, and also stands out as the most applicable to a second order calculation.
There still remain important developments to be made in each case. For example, the mode sum and effective source schemes are still under development for the Kerr gravitational case, and worldline convolution approaches have yet to be fully applied to the gravitational case for any spacetime. Finally, it should be pointed out that this review is not an exhaustive exposition of all self-force computation strategies. Many other calculations have not been described, including alternative regularization strategies [66,181,182], near-horizon waveform calculations [183][184][185][186], methods based on effective field theory [187][188][189][190], analytic calculations in particularly simple cases [95,191,192], black holes in higher dimensions [193], and calculations in non-black hole spacetimes such as wormholes [194,195] and cosmological models [196,197]. More details can be found in the reviews [60,61], and references therein.