{'fraction_non_alphanumeric': 0.03377471706711902, 'fraction_numerical': 0.008355064941641137, 'mean_word_length': 4.415878239407651, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Characteristic formulation of the Regge-Wheeler and Zerilli Green functions

We present a characteristic initial value approach to calculating the Green function of the Regge-Wheeler and Zerilli equations. We combine well-known numerical methods with newly derived initial data to obtain a scheme which can in principle be generalised to any desired order of convergence. We demonstrate the approach with implementations up to sixth-order in the grid spacing. By combining the results of our numerical code with late-time tail expansions and methods of subtracting the direct part of the Green function, we show that the scalar self-force in Schwarzschild spacetime can be computed to better accuracy than previous Green-function based approaches. We also demonstrate agreement with frequency-domain methods for computing the Green function in the gravitational case. Finally, we apply the Regge-Wheeler and Zerilli Green functions to the computation of the gravitational energy flux.


I. INTRODUCTION
While still in its infancy relative to electromagnetic astronomy, gravitational wave astronomy has already seen enormous success in recent years. The LIGO-VIRGO collaboration has reached the stage of regularly detecting the merger of stellar mass black holes and neutron stars [1]. Looking towards future detectors, the planned European Space Agency mission, LISA (Laser Interferometer Space Antenna), will provide access to an entirely new frequency range for the gravitational wave astronomy community [2]. Operating in the 0.1 -100 mHz range, LISA will not only be capable of detecting stellar mass binaries long before merging, but will also detect sources involving much larger masses. Among the new sources which will be detected for the first time are Extreme Mass Ratio Inspirals (EMRIs). These systems, comprising a stellar mass compact object such as a black hole, or even a neutron star, orbiting a massive or supermassive black hole (10 6 -10 9 M ), are expected to be found in the centres of galaxies and lie in the most sensitive part of the LISA band. There are a number of expected formation channels for EMRIs, with the most likely being dynamical friction in the dense cluster about the supermassive black hole at a galaxy's centre [3]. Stars and compact objects within this cluster undergo gravitational interactions which can lead to objects falling into a close orbit around the central black hole. Main sequence stars are unlikely to survive in this environment due to tidal disruption, but compact objects such as black holes and neutron stars can withstand these tidal forces to produce a long-lived, slowly inspiralling binary. The expected event rates for detections of EMRIs by LISA remain quite uncertain, being anywhere from 1-10 3 yr −1 [4]. However, they present some of the most exciting opportunities to study black holes, and also some of the most difficult challenges.
As with LIGO-VIRGO, LISA will identify many of the sources in its frequency band via matched-filtering, comparing detected signals to large banks of precomputed waveforms to pluck out those signals likely to be astrophysical in origin. This is because most, if not all EMRI signals will actually be so weak as to lie below the noiselevel of the detector. Thus, one of the primary goals of the next decade for the general relativity community is to generate high-accuracy waveforms for all of the expected sources. Unlike typical LIGO sources, or even other LISA sources such as massive/supermassive black hole binaries, EMRIs spend a large amount of their lifetimes in a region of the parameter space where many techniques for modelling black hole binaries are not applicable, or not currently practical. For instance, the orbit of the smaller companion object is typically not circularised, and can be highly inclined. The system is typically long lived, spending weeks or months within the LISA band. The smaller object can also reach large velocities, comparable to c, during its orbit. These properties combine to make EMRIs a unique and extremely powerful laboratory in which to test general relativity. However, they also rule out post-Newtonian theory and numerical relativity as suitable methods to model the system, except in certain limits (eg. at large radii where the companion's velocity is small, where post-Newtonian theory is applicable [5]). Instead, the standard approach to modelling EMRIs is to compute the self-force (SF) [6][7][8].
In the SF approach, the effect of the companion (mass m) is treated as a small perturbation to the background spacetime of the central black hole (mass M ). This leads to an expansion of the metric, order-by-order, in the mass ratio ≡ m M , where g is the exact metric of the full two-body spacetime, g is the metric of the background, and h (n) is the n-th order contribution to the metric perturbation.
Schematically, Fig. 1 shows the behaviour of the system. At zeroth order, the companion would follow a geodesic of the background spacetime. When higher-order effects are included, the emission of gravitational waves causes the orbit of the companion to slowly evolve and inspiral into the central black hole. The goal of the self-force approach is to compute the metric perturbation and use it to derive the inspiralling In the SF picture, the compact object orbiting the central black hole moves along a worldline (red) which deviates from a geodesic of the background spacetime (dashed, black line) due to the emission of gravitational waves. This effect can be described in terms of a (fictitious) force, Fα which we call the self-force, as it results from the particle's own bending of spacetime. The effect is exaggerated here, as in reality a background geodesic would be a good approximation to the worldline of the compact object for multiple orbits. However, the magnitude of the SF increases at points on eccentric or hyperbolic orbits which are close to periastron. motion of the companion along with the associated gravitational waveform. There are a number of technical difficulties that arise in the process. One of the most fundamental is that the metric perturbation is formally singular along the worldline of the companion. Much of the work in the SF community over the last three decades has been aimed at tackling this problem. An approach which has seen significant development in recent years is that of worldline convolution [9][10][11][12]. This approach makes use of Green functions to construct the metric perturbation by performing convolution integrals over the past worldline of the particle. Once the Green function is available, it is straightforward to use it to compute the self-force, but this merely moves the challenge to that of computing the Green function.
In this work, we develop a characteristic formulation of the Green function for the Regge-Wheeler and Zerilli equations, which govern perturbations of Schwarzschild spacetime. Similar approaches have previously applied a characteristic method to compute the Green function for the scalar wave equation [13,14]. Our approach more fully develops the characteristic formalism, extending these earlier works to be applicable to the Regge-Wheeler and Zerilli equations, and also to arbitrary convergence order in the numerical scheme. The layout of the paper is as follows: in Sec. II we review the theory of Green functions for the Regge-Wheeler equation; in Sec. III we describe our numerical method for solving a characteristic initial value problem for the Green function; in Sec. IV we derive initial conditions that enable the method to be extended to arbitrary convergence order; in Sec. V we demonstrate our method by applying it to compute the scalar self-force and the gravitational energy flux. Finally, we provide some concluding remarks in Sec. VI.
Throughout this work we use geometrised units such that the speed of light and the gravitational constant are set to unity (G = c = 1).

II. GREEN FUNCTIONS FOR THE REGGE-WHEELER FORMALISM
The Regge-Wheeler formalism is based on constructing solutions to the linearised Einstein equations from solutions to the modified scalar wave equation, where s is the spin of the field (s = 0 for scalar fields, s = 1 for electromagnetic fields and s = 2 for gravitational fields). Working in Schwarzschild coordinates, {t, r, θ, φ}, this equation is separable using the ansatz where Y m (θ, φ) are the spherical harmonics. Transforming to double null coordinates, {u, v} = {t − r * , t + r * } where r * = r + 2M ln( r 2M − 1) is the Schwarzschild tortoise coordinate, the Regge-Wheeler master function Ψ s m satisfies the Regge-Wheeler equation, Here, the modes of the source are defined in the same way as the master function, We are interested in constructing the retarded Green function, which satisfies where δ (4) (x, x ) is the invariant Dirac-delta distribution. Primed coordinates denote the position of the particle, while unprimed coordinates denote a point in the particle's past history. The retarded Green function can be written in terms of its decomposition into spherical harmonic modes, Note that time translation invariance means the retarded Green function only depends on the time difference, ∆t ≡ t − t. It is convenient to exploit spherical symmetry by using the addition theorem for the spherical harmonics to rewrite the mode decomposition in terms of the angle γ between x and x , (2 + 1)P (cos γ)G ret s (r, r ; ∆t).
Substituting into Eq. (6) and using the completeness relation for the spherical harmonics, we find that the modes of the retarded Green function satisfy the Regge-Wheeler equation with a distributional source, The modes of the retarded Green function for the Regge-Wheeler equation are thus Green functions for the flatspace 2D wave equation with a potential. Writing these as we find that g s (u, v) satisfies the homogeneous Regge-Wheeler equation, with characteristic initial conditions

III. NUMERICAL SOLUTION OF THE CHARACTERISTIC INITIAL VALUE PROBLEM
The characteristic initial value problem represented by Eqs. (11) and Eqs. (12) is well-suited to a numerical treatment. All of the numerical methods, as well as the method of deriving initial conditions, are applicable to any equation of the form of a flat-space wave equation with a potential. As such, we formulate the scheme for a generic potential, P (u, v), but remind the reader that in the case of the Regge-Wheeler equation this potential is given by where r is a function of v − u.
In addition, in the gravitational case (s = 2) we will consider the Zerilli potential, The Zerilli equation, Eq. (11) with the potential replaced by P Zer (u, v), governs the even parity perturbations to the Schwarzschild spacetime, while the odd parity perturbations obey the Regge-Wheeler equation with s = 2.
In self-force applications this is the location at which the selffield and self-force components for a given worldline are computed. The domain contains a large family of past worldlines which pass through this point, several of which are shown [projected onto the 2D (u, v) submanifold]. These include a circular orbit (blue), and a number of geodesics with increasing eccentricity. Thus for a single run of the numerical code, we can compute the self-force at this point for a large number of past worldlines, as opposed to at all points on a given worldline as would be done with other methods. Note, the horizon is to the left in this plot, while radial infinity is to the right.

A. Numerical integration stencils
The retarded Green function only has support inside the past light cone, so we need only consider the value of g s and its derivatives on and inside the lightcone. This suggests a natural manner in which to subdivide our domain into a grid with spacing h, using lines of constant u or v, see Fig. 2. Doing so, we can derive a numerical scheme which will converge as a desired power of h. We begin by considering a single cell, see Fig. 3. Integrating Eq. (11) over this cell yields (omitting g's dependence on where we have introduced the notation to denote the value of the function g s at the point Note that the first term in the integral is exact. The second term, however, cannot be integrated exactly and must be evaluated numerically with an approximation valid to some order in h. Using methods from [15,16] there is a clear, systematic way to calculate the integral to any desired order, O(h n+2 ), where n is even. With this error in a single cell, we can ensure the global error is O(h n ) after integrating over O(h −2 ) intervals. The fundamental idea behind the method is to Taylor expand the integrand, which we denote by H(u, v) = P (u, v)g s (u, v) for conciseness, about the centre point of a given cell, We obtain a factor of h 2 from the double integral over the cell, and so we need only take this expansion to O(h n ), where n is the desired order of global accuracy. This will result in 1 2 n(n + 1) unknown coefficients H i,j , which can be solved for by using the value of the function at known points.
Following this approach, a second order (n = 2) algorithm is given by This scheme has a simple stencil, shown in the top panel of Fig. 3. Assuming we have already solved for, or provided a priori, the values of g s at the points (u − h, v), The overall numerical implementation of this algorithm is then straightforward. We provide initial data  18) is performed. In the higher order schemes, we show only the coordinate which changes along each ray, for simplicity. Note that in each of these diagrams, the lightcone is toward the top, so the lower points are further in the past relative to the base point at which we will ultimately compute the self-field and self-force components.
along the lightcone (described in more detail in Sec. IV), and evolve through the domain. The most straightforward evolution method is to evolve down each ray of fixed u or v, then move on to the next, u + h or v + h. We have implemented this algorithm in a C code, with two rays stored in memory at any one time, and specific rays saved for later output. This allows for the code to use comparatively little memory, even for large domains and fine resolutions. Higher order schemes are straightforward to obtain us-ing the same method. We have taken this to sixth order as a demonstration, with the algorithms derived using the computer algebra software Mathematica. The only obstacle to even higher order schemes is the computa-tional cost of the derivation, though it is not expected that anything beyond eighth order would prove necessary for most practical applications. A fourth order algorithm is given by (2) s + 10 P 11 g 11 s + P 01 g 01 s + P 10 g 10 s − 4 P 20 g 20 s + P 02 g 02 s + P 30 g 30 s + P 03 g 03 s − P 12 g 12 while a sixth order algorithm is given by The stencils for these schemes are shown in the lower panels of Fig. 3.
Note that in the fourth order scheme, the second order approximation appears. This is because the point at which we wish to calculate the value of the Green function is being used in the derivation of the algorithm. As we do not yet know the value at this point, we must provide it to a suitable approximation. Given that it is multiplied by an overall factor of h 2 , we can use the second order scheme, and the overall error will remain O(h 6 ). We use a similar approach for the sixth order scheme, with the fourth order approximation used in the algorithm.
Our algorithms are by no means unique. A notable alternative is the predictor-corrector method used in [15]. However, the Lorenz gauge equations solved there involve first derivatives, which reduce the order by the factor of h. We can thus avoid employing such a scheme, which would impact numerical efficiency due to being an iterative method, as we do not have any such first order derivatives in the Regge-Wheeler or Zerilli equations.
As is clear from Fig. 3, as we go to higher order we require a greater amount of past information in order to calculate the value of the Green function at a given point. This has significant implications for the initial data which must be provided, as will be discussed in detail in Sec. IV. Aside from initial data issues, the evolution through the numerical domain can be implemented in the same manner as for the second order scheme, though requiring more rays to be stored in memory during the computation. However, this scales only as n, whereas the scheme is accurate to O(h n ). Thus moving to higher order schemes, while computationally slower due to the increased number of operations, is not a significant additional burden on memory resources.

B. Calculation of derivatives of the Green function
In addition to solving for the Green function, we can also use the method outlined above to solve directly for the derivatives of the Green function, as opposed to computing such derivatives in post-processing by means of finite difference. For the first derivative with respect to the base point (the vertex of the lightcone),which is required for computing the self force, we can simply apply the above algorithms to equations for the derivatives, Higher order derivatives with respect to the base point are equally straightforward. However, we can also compute derivatives with respect to the field point (a point inside the past worldline), as well as mixed derivatives with respect to both the base point and field point. Applying the desired differential operator to the generalised equation, for example, we have equations to which our schemes may be applied. g s (u, v) and lower order derivatives that appear can be computed alongside the desired function. Further combinations of base point and field point derivatives, as well as higher order derivatives can be obtained in a similar fashion.

C. Convergence
The schemes outlined here have been implemented in a C code, up to sixth order for the Regge-Wheeler Green function, and fourth order for its derivatives. The Zerilli Green function and its derivatives have been taken to second order, though extending these to higher orders is straightforward. The convergence order of these schemes (once sufficiently accurate initial data is provided) can be verified by comparing numerical solutions using three different resolutions: h L , h M and h H , which we denote by g The value of n for a single mode of the Regge-Wheeler Green function computed using each of the three numerical algorithms is shown in Fig. 4, and verifies that they do indeed converge at the expected order. The black, dashed lines indicate the expected convergence orders for the schemes. The spikes in the convergence order are associated with zero-crossings, and are not unusual in numerical methods. Note, in addition, the noise in the sixth order scheme. This is due to roundoff error as we go to higher resolutions.

D. Code optimisation considerations
Two additional steps have been taken to improve the efficiency and computation time of the code.
1. OpenMP parallelization has been implemented, with independent -modes computed in parallel.
2. Given that the potential, P and f are effectively functions of v − u, not general functions of (u, v), they need only be computed at points on the lightcone, as a vertical line through the domain shown in Fig. 2 marks a line of constant r. This also holds for the computation of r from r * , which is performed numerically and is the most expensive individual pointwise calculation performed in the code. Thus, significant gains in computational cost can be obtained by precomputing these before evolution through the domain.
The current implementation of our code can compute 101 modes, using the fourth order scheme with h = 10 −2 , up to ∆t = 240 in the past of the base point, and parallelized over 5 threads in 217 s, on a desktop computer with a Ryzen 7 processor. Scripts to automate the calculation over a range of base points make it possible to compute the Green function in a large region outside a Schwarzschild black hole, though this does lead to a substantial data storage problem. We expect the use of reduced order surrogate models to provide a solution to this in the future.

IV. INITIAL DATA
With the numerical scheme outlined, it remains to provide suitable initial data to compute the retarded Green function. Previous implementations of similar numerical schemes have not been concerned with initial data, as they chose to allow junk radiation sourced by inconsistent initial data to radiate away, and only consider the solution far from the lightcone where the initial data is imposed [15,16]. However, for the purposes of computing the Green function using a characteristic initial value formulation, the early time behaviour near the the light cone is crucial to obtaining the correct Green function. As discussed in Sec. III, the higher the order of the scheme, the greater the number of values of u and v at which the value of the Green function must be known a priori. As a result, for higher orders, we must provide initial data on rays inside the lightcone, not just on the lightcone itself.
Providing values on the lightcone alone would only be sufficient for the second order scheme. The value of g s (u, v; u , v ) along null rays connecting (u, v) and (u , v ) is known to be 1 2 , and this alone provides sufficient initial data for the second order scheme.
For the higher order schemes, rather than solve for the value of the Green function on points inside the lightcone exactly, we use an expansion of g s near the lightcone to the same order in h as our numerical scheme. Writing where ∆u = u − u scales as h for rays near the ∆u = 0 side of the lightcone and ∆v scales as h for rays near the ∆v = 0 side of the lightcone, this approximation is accu-rate to O(h n ) 1 . Since the initial data is only required on a small, h-independent number of rays near the lightcone, this would then be sufficient to ensure that the global error remains O(h n ).
We can derive analytic expressions for the V k 's for the both Regge-Wheeler and Zerilli functions by substituting the above expansion into the relevant equation and solving the transport equations for each coefficient, order by order. Considering the equation with a general potential and demanding that it be satisfied at each power of ∆u∆v 2 we obtain the transport equations (k + 1)(∆u∂ u + ∆v∂ v + (k + 1))V k+1 = −2 ∂ 2 uv − P V k (25) along with the initial condition V 0 = 1 2 for the value of the Green function on the lightcone.
For higher order terms, we employ an ansatz for each term, and solve the resulting equations for the coefficients p km 1 Near the the vertex of the lightcone, both ∆u and ∆v scale as h so the approximation is even better and is accurate to O(h 2n ).
using Mathematica. We applied this procedure to compute up to V RW 8 and V Zer 6 . We do not give the higher order coefficients here due to the length of their expressions, but have provided them electronically as supplemental material. Although deriving further V k 's is trivial, we chose to stop at V RW 8 as that is already more than is required for our sixth-order numerical algorithm.
Note that there is no t dependence in V RW/Zer 1 or V RW/Zer 2 . In fact, it is trivial to show, by induction, that V RW/Zer k depends only on r for all k; the t dependence in the initial conditions is contained entirely within the powers of − ∆u∆v 2 . This is what allowed for the transport equations to be simplified to the point that they could be solved analytically.

A. Initial data near the vertex of the lightcone
The expressions for the V RW/Zer k 's appear to be singular at coincidence, where the denominator ∆r * vanishes. However, non-singularity is explicitly enforced when solving for them so we know that the apparent singularity will be cancelled by the numerators also vanishing. This cancellation is not explicit due to the non-trivial relationship between r and r * . In order to avoid this becoming an issue in the numerical implementation, where cancellation is only every guaranteed to the level of round-off, we employ a series expansion about coincidence for points very near the lightcone. Figure 5 shows the effect of this numerical cancellation and its resolution via our use of a series expansion. versus radius, for the case of the vertex of the lightcone located at r = 6. The behaviour at r = r can be avoided by use of the series expansion.

B. Initial data for the derivatives of the Green function
By differentiating the functions V k , we can obtain initial data for the derivatives of the Green function. This requires deriving the coefficients for one order higher than for the Green function, but there is no significant challenge to this. For instance, for the first u derivative, say, we have the expression Having lost an order in h in the first term in the sum, we must take this sum to at least k = 6 for fourth-order convergence, k = 8 for sixth-order, etc. A similar equation holds for the v derivative. Computation of the derivatives of the functions V k , as well as series expansions close to coincidence, is straightforward. It is also possible to obtain initial data for the higher order and mixed derivatives of the Green function, by simply applying the relevant differential operator to Eq. (24).

V. NUMERICAL RESULTS
We have implemented the schemes outlined in Sec. III along with the initial data derived in Sec. IV as both Mathematica and C codes. In Fig. 6 we show some representative Regge-Wheeler results from the C code. These were obtained by computing the modes g s using our numerical code and following Ref. [11] in performing a smoothed sum over , (2 +1)P (cos γ)G ret s (r, r ; ∆t)e − cut 2 , (31) where we chose max = 100 and cut = 20. This smoothing is required to obtain a convergent mode-sum, but in turn introduces another problem.
To understand this, consider the Hadamard decomposition of the retarded Green function, which is valid in a causal domain of the point x (i.e. where x and x are connected by a unique geodesic). Here, the term involving U (x, x ) is known as the direct part, and the term involving V (x, x ) is the tail part. As we are constructing the Green function using a smoothed sum, the δ-distribution in the direct part is effectively smeared into a Gaussian of finite width centred on the base point. This can clearly be seen in the early-time behaviour of the red dashed curve in Fig. 6. To circumvent this, we use the method developed in [17], where it was shown that an -mode decomposition of the direct part can be subtracted mode-by-mode from the modes of the full retarded Green function, G ret s . Then, summing over modes we get the full retarded Green function minus its direct part. This regularised Green function is no longer contaminated by the smeared direct part, as illustrated by the solid black curve in Fig. 6.
Our results are consistent with other approaches to computing the Green function, as demonstrated in Figs. 7 and 8, where we compare with two other approaches to computing it: an expansion valid for late times (Fig. 7) and a numerical calculation based on an inverse Fourier transform of the frequency-domain Green function (Fig. 8). Indeed at sufficiently late times we can save computational cost in Green function calculations by using the late-time expansions in place of our numerically computed solutions.
From Fig. 6, it may appear that there is some nonnegligible disagreement between the retarded solution and the late-time tail, particularly in the r derivative of the Green function. However, we can verify that this   6. The scalar Green function, G ret 0 (x, x ) and its derivatives for points separated along an eccentric orbit with p = 7.2 and e = 0.5, with base point r = 6. The effect of subtracting the direct modes using the method of [17] to obtain only the tail contribution can be clearly seen in each case. At late times, the numerical solutions also agree with a late-time approximation of the Green function obtained using the method of [18].
is consistent with expectations. The leading order term of the late-time tail is t −3 and the expansion we use is accurate up to t −5 log(t), taken from [18]. Thus the latetime tail for the derivatives should be accurate up to t −6 log(t). We can see in Fig. 7, that as we go to later times, the agreement between the retarded solution and the late-time tail scales as the next term in the tail expansion. Thus, if greater agreement with the late-time tail is required, it is a simple matter to define a larger numerical domain and compute the retarded modes to later times. In the self-force calculation described below, however, we found that the relative change from altering the matching time between the retarded solution and late-time expansion is below our numerical error in the resulting self-field and self-force values.

A. Scalar self-force
As a demonstration of the utility of our method, we now show how it can be used to compute the scalar selfforce. We obtain results for the self-field and self-forces along both a circular and an eccentric orbit. Both orbits considered lie close to the separatrix between bound and unbound orbits, with p = 7.2, e = 0.5 in the eccentric case and p = 6 in the circular case.
Despite only seeing much of its progress recently [10], the worldline convolution method we use was one of the first methods proposed to tackle the self-force problem.
In two independent works [19,20] it was shown that the self-force is given by the so-called MiSaTaQuWa equation, which in the scalar field case and for geodesic motion is given by The integral here is over the past worldline of the compact object, truncated at a proper time just before reaching its current location. It is clear that due to the cutoff on the upper limit of integration, only the tail part of the Green function contributes to the self-force. It is also clear that the Hadamard parametrix will be inadequate: the Hadamard decomposition is restricted to a causal domain whereas the integral in Eq. (33) extends over the entire past worldline, which is not confined to a causal domain. Our numerically calculated regularised Green function does not suffer from this limitation and is therefore ideally suited to computing the self-force using the worldline convolution approach.
With the regularised Green function in hand (and a late-time tail expansion matched on at very late times), we can then compute the self-field and self-force components by straightforward numerical integration. The results of doing so are shown in Table I 10 t'-t t'-t that a high degree of accuracy has been achieved. These results were obtained with max = 100, cut = 20, matching the tail on at ∆t = 200, and integrating the tail to ∆t = 2000. We also used a numerical resolution of h = 0.01 and sampled the numerical solution with resolution of 0.1 (i.e. every 10 steps).
After performing an analysis on the primary contributions to the numerical error in our results, the dominant sources were found to be the values of max and cut . It is worth noting, however, that we already obtain results of comparable accuracy (particularly in the circular case) to [12], with only half as many -modes used. We attribute this to use of a characteristic initial value formulation, rather than the time-domain Gaussian approach used in [12]. This can be verified by replacing our characteristic initial data with a narrow Gaussian of width 0.1, centred at the base point, and matching on a quasilocal expansion as described in [21] at ∆t = 14, whereby we obtain relative errors of of O(1)% for the circular orbit self-field when 100 -modes are considered.
The derivatives of the direct modes defined in [17] are another significant contribution to the error. For convenience, these were computed using finite-differencing (and differentiation of series expansions at early times). We found that finite differencing error introduced a significant amount of noise if too small a step in t or r was chosen or if too high an order of finite differencing was used. Using larger values of ∆t or ∆r , or a lower order finite difference stencil reduced the noise in the derivatives and lead to changes in the results at O(10 −2 )%. Ideally, we would compute these directly rather than via finite difference, but this is beyond the scope of this work. Another source of error which contributed at a similar order is the sampling resolution. Increasing the sampling resolution improves the results by O(10 −2 )%.
The next most significant contributions considered include the time at which the tail is matched on to the numerical results, the maximum time to which the tail is integrated, and the numerical resolution. These all contribute at O(10 −3 )% only.
Finally, the orders of the Bessel series expansions (used for the retarded modes at very early times, as described in [17]) and the series expansion for the direct modes (again, described in detail in [17]) contribute negligibly to the error, at O(10 −5 ) and O(10 −8 )% respectively.

B. Gravitational Green function
In addition to the scalar calculations demonstrated thus far, we can also solve for modes of the Green function and its derivatives in the s = 2 gravitational case. While no formal procedure yet exists for subtracting the direct modes or matching on a quasilocal piece at early times, we can still draw comparisons between the results here and frequency domain solutions for the vacuum Regge-Wheeler equation. Some sample results are show in Fig. 8. These results, for = 2 along a circular orbit with r = 6, show that there is excellent agreement between both the time domain and frequency domain approaches. While there looks to be larger disagreement at early times, it can be seen in the accompanying plot that increasing the maximum frequency integrated to in the Fourier transformation from frequency domain to time domain improves the agreement at these early times. It has also been observed that the magnitude of the resid-ual increases with , but again improves with larger maximum frequency. Thus, the primary source of the residual is the accuracy of the frequency domain results, not the time domain results.
With the accuracy of the gravitational results verified, we next apply the Regge-Wheeler and Zerilli Green functions to the computation of the local force F t on a gravitational perturbation. This is directly related to the gravitational energy flux radiated to null infinity and into the horizon, and we use this fact to provide robust reference values against which to check our results. As demonstration of the method, we choose a circular orbit, r = 10, and compute the total flux for specific ( , m) modes,Ė 21 , E 22 . This requires no regularisation of the individual modes, and can be constructed easily from the master functions, Ψ RW/Zer (t , r * ), and their derivatives. In this work, we compute the Moncrief versions of the master functions, detailed in [22].
To obtain these, we must integrate the Green function against the source. The point particle source for the Regge-Wheeler master function, in general, takes the form S(t , r * ) = s 1 (t, r * )δ(r * −r * )+s 2 (t, r * )∂ r * δ(r * −r * ). (34) Using integration by parts and Eq. (10), we obtain an expression for, Ψ at the base point, g 2, (r * , r * ; t) ∂ r * s 2 (t, r * )| r * =r * − s 1 (t, r * ) + ∂ r * g 2, (r * , r * ; t )| r * =r * s 2 (t, r * ) dt This expression holds for both the Regge-Wheeler-Moncrief and Zerilli-Moncrief master functions, with g 2, being the mode of the s = 2 case of the corresponding Green function. The ± comes from the Heaviside step functions in Eq. (10), and correctly captures the known jump in the value of the master function at the location of the particle. Derivatives of Ψ are also required to calculate the flux, but can be easily obtained using the same procedure as for Ψ, and computing the derivative before taking the limit to the particle. The functions s 1 and s 2 , which contain all of the mmode dependence of Ψ, can be obtained in terms of the functionsF ,G defined in [23]. Accounting for the authors' convention of fully evaluating the source in terms of r rather than r * , expressions for s 1 and s 2 can be obtained, s 1 (t, r * ) = f −1G (t, r * ) + ∂ r * f f −2F (t, r * ) (36a) s 2 (t, r * ) = f −2F (t, r * ).
When restricted to a circular orbit, the functionsF and G are given byF odd m = 8πµ λ(λ + 1) where E and L denote the energy and angular momentum along the circular orbit, and Y , Y φφ and X φ denote the even sector scalar, tensor and odd sector vector harmonics, respectively, [22]. An overbar denotes complex conjugation. All instances of f and Λ in Eqs. (36a)-(37d) are evaluated at r , rather than r. Note that in this convention the ∂ r * s 2 (t, r * ) term in Eq. (35) vanishes. We compute second-order numerical values for the = 2 mode of the Regge-Wheeler and Zerilli Green functions, as well as all necessary derivatives (note, this computation requires mixed higher order derivatives). Only the third derivative ∂ r r r g 2, is computed using finite differencing, all others are computed directly. We thus obtain the flux values shown in Table II. We expect use of higher order schemes for the Green functions and all relevant derivatives will lead to improved accuracy in the flux calculation.

VI. CONCLUSIONS
In this work, we have implemented and demonstrated an efficient means of computing the retarded Green func-tion for the Regge-Wheeler and Zerilli equations. We derive, for the first time, initial data inside the lightcone to arbitrary order in the distance from the lightcone. This initial data can be used to seed standard numerical approaches.
We also implement methods of removing the direct part of the scalar Green function and its derivatives for the first time, resulting in the computation of accurate self-field and self-force values for both a circular and eccentric orbit. In future works, this method can be trivially applied to compute the self-force along more complex worldlines, such as a hyperbolic encounter with the central black hole.
The improvements over previous methods are apparent, with the number of modes required approximately halved, without a significant loss of accuracy. The numerical scheme and initial data may also be systematically extended to any numerical order in the grid-spacing h, allowing for potentially further gains in efficiency, via coarser grids which yield the same accuracy. Additionally, a large portion of the calculations detailed here were performed in Mathematica. Converting these to C may take time, and require the use of third-party libraries such as the GNU Scientific Library. The gain in computational speed, however, would make this a worthwhile effort.
Currently, the greatest outstanding issue is the means of removing the direct part of the Green function in the gravitational case. With this done, the results could be applied to compute full metric perturbations in the Regge-Wheeler gauge, and gauge invariants used to draw comparison with other approaches to the two-body problem.
Finally, many of the methods developed in this work can be applied to other problems, such as the Lorenz gauge wave equation for metric perturbations of Schwarzschild spacetime, or even the the Teukolsky equation. The efficiency of the code also makes other numerically intense applications more feasible, such as a selfconsistent evolution scheme.