Frequencydomain calculation of the self force:
The highfrequency problem and its resolution
Abstract
The modesum method provides a practical means for calculating the self force acting on a small particle orbiting a larger black hole. In this method, one first computes the sphericalharmonic mode contributions of the “full force” field , evaluated at the particle’s location, and then sums over subject to a certain regularization procedure. In the frequencydomain variant of this scheme the quantities are obtained by fully decomposing the particle’s self field into Fourierharmonic modes , calculating the contribution of each such mode to , and then summing over and for given . This procedure has the advantage that one only encounters ordinary differential equations. However, for eccentric orbits, the sum over is found to converge badly at the particle’s location. This problem (reminiscent of the familiar Gibbs phenomenon of Fourier analysis) results from the discontinuity of the timedomain field at the particle’s worldline. Here we propose a simple and practical method to resolve this problem. The method utilizes the homogeneous modes of the self field to construct (rather than the inhomogeneous modes, as in the standard method), which guarantees an exponentiallyfast convergence to the correct value of , even at the particle’s location. We illustrate the application of the method with the example of the monopole scalarfield perturbation from a scalar charge in an eccentric orbit around a Schwarzschild black hole. Our method, however, should be applicable to a wider range of problems, including the calculation of the gravitational selfforce using either Teukolsky’s formalism, or a direct integration of the metric perturbation equations.
I Introduction
The problem of calculating the gravitational selfforce MST (); QW (); GW () acting on a pointlike test particle as it moves in orbit around a black hole is attracting considerable attention in recent years Poisson:2003nc (); Lousto:2005an (). Within this context, various authors have been studying also the analogous problem of the scalarfield selfforce Q () acting on a particle endowed with a scalar charge, which proved to be a useful toy model. The electromagnetic (EM) selfforce, acting on an electric point charge, was also studied by various authors following the seminal work by DeWitt and Brehme DeWittBrehme (). A practical algorithm, commonly used for computing the self force in all three cases (gravitational, EM and scalar), is the modesum method Barack:1999wf (); Barack:2001gx (); Barack:2002mh (). This method requires as input the multipole modes of the full (retarded) perturbation fields, along with their derivatives, evaluated at the particle’s location. These multipoles can be calculated using either frequencydomain methods (as in, e.g., Burko:2000xx (); Detweiler:2002gi (); Keidl:2006 ()) or timedomain numerical evolution (as in, e.g., Barack:2000zq (); Barack:2002ku (); Haas:2007kz (); Barack:2007tm ()). In both computational approaches one sets off by writing down the appropriate set of perturbation equations, modeling the source term associated with the point particle (namely, the energymomentum, the electric fourcurrent, or the scalar charge) as a deltafunction distribution. In the frequencydomain approach one then decomposes the inhomogeneous perturbation equations into Fourierharmonic modes (“ modes”) and proceeds by solving the resulting ordinary differential equations (ODEs) with suitable boundary conditions at spatial infinity and at the event horizon. In the alternative, timedomain approach, one refrains from decomposing the field into frequency modes, and instead tackles the partial differential equations for each directly using time evolution.
Each of the above two approaches has its strengths and weaknesses. The timedomain approach has the advantage that one only deals with a single field for each , whereas in the frequencydomain approach one has to sum over the various modes. On the other hand, the latter approach has the obvious advantage that one only faces ODEs. Despite the fact that timedomain methods are winning growing popularity in recent years, frequencydomain calculations remain an appealing option for some range of orbital parameters Barton (). Also, it turns out that the nonradiative multipoles of the gravitational perturbation in Schwarzschild (i.e., the modes ) are difficult to analyze in the time domain, due to instabilities Barack:2002ku (); Barack:2007tm (), and one resorts to a frequencydomain calculation in this case BGSprep (). Working in the frequency domain, however, brings about a technical issue which, to the best of our knowledge, has not been addressed so far in the current context.
To illustrate the problem, is it instructive to refer to the simple case of (minimallycoupled and massless) scalarfield perturbations from a pointlike scalar charge orbiting a Schwarzschild black hole. In this case, the scalar field can be decomposed in spherical harmonics , yielding the multipolar mode functions . Here are the standard Schwarzschild coordinates. Let us denote the value of the particle’s location at time by . With suitable boundary/initial conditions, a unique solution is obtained for (for each ), which is continuous along . However, the derivatives and will generally suffer a finite jump discontinuity across , which reflects the presence of a source “shell” resulting from decomposing the point charge in spherical harmonics. In particular, if the orbit is eccentric, the derivatives of will generally be discontinuous functions of at a fixed value of along the orbit.
Now imagine trying to reconstruct (for some fixed along the orbit) as a sum over its Fourier components:
(1) 
Since, for an eccentric orbit, is only a function of at the particle’s worldline, it follows from standard Fourier theory James () that the Fourier sum in Eq. (1) will only converge there like . The actual situation is even worse, however, because for selfforce calculations we need not only but also its derivatives. For instance, to calculate the component of the self force we need to evaluate (onesided limits of) as approaches . Suppose again that we want to reconstruct this quantity from its Fourier components (namely the functions ). Since is a discontinuous function of , we will inevitably face here the well known “Gibbs phenomenon” Gibbs (). Namely, the Fourier sum will fail to converge to the right value at . (The problematic behavior of the Fourier sum is simply a consequence of our attempt to construct a discontinuous function—or a nonsmooth function in the case of the field itself—from a sum over smooth harmonics.)
From a practical point of view this would mean that (i) at the coincidence limit the sum over modes would fail to yield the correct onesided values of , however many modes are included in the sum; and (ii) if we reconstruct at a point off the worldline, then the Fourier series should indeed converge; Alas, the number of modes required for achieving a prescribed precision would grow unboundedly as approaches , making it extremely difficult to evaluate at the coincidence limit.
This technical difficulty is rather generic, and will show also in calculations of the local EM and gravitational fields. Consider, as a second example, the gravitational perturbation from a point mass moving in an eccentric orbit in Schwarzschild: In suitable gauges (like the Lorenz gauge, often applied in selfforce calculations) the multipole modes of the physical metric perturbation ^{1}^{1}1Formally, the relevant multipoles in this case are not the standard spherical harmonics but rather the 2ndrank tensorial harmonics. Similarly, in the discussion below regarding Teukolsky fields in Schwarzschild, the relevant multipoles are the spinweighted spherical harmonics. These technical details do not in any way affect the discussion here. are again functions of and along the orbit, and their derivatives are generally discontinuous there. Attempting to construct them naively from a sum over frequency modes would encounter the same difficulty as in the scalar case: A poor convergence of the (modes of the) metric perturbations, and lack of convergence for their derivatives.
For orbits in Kerr spacetime the situation is basically similar though more subtle. The Kerr variant of the modesum method Barack:2002mh () requires, just as in its Schwarzschild counterpart, the spherical harmonic modes of the perturbation field (as well as their derivatives) as input.^{2}^{2}2Although the full separation of the field equation in Kerr is based on the spheroidal harmonics, the mode sum method (at least in its present form Barack:2002mh ()) is based on the spherical harmonics. Note that the latter harmonics do provide a valid mode decomposition even in Kerr, even though the field equation couples modes of different sphericalharmonic . One of the possible practical ways to construct a Fouriersphericalharmonic mode is by solving for the separable spheroidalharmonic modes for various , and then summing over their contributions to , as discussed in Ref. Barack:2002mh (). For given , the sphericalharmonic decomposition of the point charge will again result in a functiontype source term distributed over a shell, which in turn renders the derivatives of discontinuous. Therefore an attempt to construct (and, more crucially, its derivatives) through a naive summation over its Fourier modes will lead to the same difficulties as in the Schwarzschild case.
The problem discussed here takes an even more extreme form when considering EM or gravitational perturbations via the Teukolsky formalism: Here, the modes of the perturbation fields (now the NewmanPenrose fields or ) are not even continuous at the particle’s orbit—a consequence of the fact that the source term for Teukolsky’s equation involves derivatives of the electric fourcurrent or the energymomentum tensor associated with the particle (a single derivative in the EM case; a second derivative in the gravitational case).^{3}^{3}3Note also that the gravitational selfforce calculation requires the local metric perturbation, and the construction of the latter requires one to apply certain differential operators to the Teukolsky fields or nowski (). The situation with the EM selfforce is similar, as it again requires applying certain differential operators to or . Again, a naive attempt to construct these multipoles as a sum over their modes will be hampered by the Gibbs phenomenon, and the associated lack of convergence.
In this article we propose a way around the above problem, which is both elegant and extremely simple. In our method we use the homogeneous radial functions (extended all the way through to the particle’s worldline), instead of the actual inhomogeneous functions. The Fourier sum of these homogeneous radial functions is found to converge exponentiallyfast, and to yield the correct values of the perturbation multipoles (and their derivatives) along the particle’s worldline. We shall focus in this paper on the scalarfield case. We justify our new method using simple analytical arguments, and then demonstrate the validity of the method (and the exponential convergence) with an explicit numerical calculation in the case . The same method should be applicable, however, for any of the other problems mentioned above: EM and gravitational perturbations using Teukolsky’s equation (or Sasaki–Nakamura’s equation), as well as metric perturbations in the Lorenz gauge. A forthcoming paper BGSprep () will report on the computation of the local monopole and dipole modes of the Lorenzgauge perturbation (for eccentric orbits in Schwarzschild), facilitated by the new method suggested here.
This paper is structured as follows. In Sec. II we set up the physical scenario—a pointlike scalar charge in a bound orbit around a Schwarzschild black hole—and review the formalism commonly used in this case to construct the scalarfield multipoles and the scalar selfforce. Section III demonstrates how the naive sum over frequencies leads to the Gibbs phenomenon and to the associated problematic convergence at the particle’s location. Then in Sec. IV we present our new method of extended homogeneous solutions, and show how it cures the problematic behavior of the Fourier sum. We provide the theoretical justification to this method, as well as numerical verification in the monopole () case. In Sec. V we highlight the advantages of the new method and discuss foreseeable applications. Appendices A–C give details of the methods used for our numerical illustrations, and App. D contains some technical details relating to the formal justification of our new method.
Throughout this work we use standard geometrized units (with ) and metric signature .
Ii Preliminaries
ii.1 Physical setup and scalarfield equation
Consider a pointlike particle which moves on an eccentric, bound geodesic orbit around a Schwarzschild black hole with mass parameter . The particle’s worldline is denoted , where is the proper time. The particle’s trajectory is bounded within the range for certain and . Without loss of generality we shall take the orbit to be equatorial, namely .
Assume now that the particle carries a scalar charge . This charge couples to a massless, minimallycoupled scalar field , satisfying the field equation
(2) 
Here is the scalar charge density, which takes the form of a function along the particle’s worldline:
(3) 
where is the metric determinant, and hereafter the vectorial indices of and are suppressed for brevity.
Since is timelike (hence monotonic) we can use it instead of to parametrize the orbit. In the plane the orbit is then denoted by . Transforming the integration variable in Eq. (3) from to and substituting , we find
(4) 
where . Note that only depends on : We have , where is a constant of motion.
ii.2 Sphericalharmonics decomposition
We now separate the field equation (2) by decomposing in spherical harmonics, in the form
(5) 
Here are the standard (complexvalued) normalized spherical harmonics given by , where are the associated Legendre polynomials and are (real) normalization constants. The factor is introduced for later convenience. The charge density in Eq. (4) is decomposed in a similar manner:
(6) 
where
(7) 
Here , is an element of solid angle, and an asterisk denotes complex conjugation. The field equation (2) is now separated, and for each the function satisfies the partial differential equation
(8) 
where , is the tortoise radial coordinate, and
(9) 
The function is determined for each by Eq. (8), supplemented with suitable boundary conditions at null infinity (no incoming waves) and at the event horizon (no outgoing waves). Since the source term has a function form, the function must be continuous at , but it will generally fail to be differentiable there: Its derivative (and also its derivative, except at the two orbital turning points) will suffer a discontinuity along the orbit. On each side of the worldline , however, the function satisfies the homogeneous part of Eq. (8), and, since the homogeneous field equation and the curve are both analytic, we expect to be analytic (in both and ) anywhere off the worldline. We shall assume this analyticity here, although we are not aware of a proof. The alternative option appears highly unlikely, because it would mean that the actual field produced by the point charge somehow develops irregularities in the vacuum region off the particle.
ii.3 Self force via mode sum
Once the functions have been determined (e.g. numerically), the self force acting on the scalar charge may be constructed by the modesum method. This procedure is described in detail in Refs. Barack:1999wf (); Barack:2001gx (); Barack:2002mh (). Here we outline it briefly, in order to provide some perspective on how the quantities are incorporated in the construction of the self force.
Let denote the contribution of an individual to :
(10) 
The component of the fullforce field, , is then obtained by applying a certain linear differential operator to . [ is the same differential operator which determines the force that a smooth, “nonself” field would exert on a test charge.] In the scalarfield case we have , and therefore ^{4}^{4}4Strictly speaking, has been defined in Ref. Barack:1999wf () to be the multipole of , rather than . Here, however, the component of the self force vanishes as we consider an equatorial orbit; and for the three remaining components the two definitions coincide.
(11) 
The quantities are to be evaluated at the particle’s location. To be more specific, let us denote by the event (on the particle’s worldline) at which the self force is to be evaluated. Then the quantities are to be evaluated at (a somewhat rough statement which will be refined shortly). From Eq. (11) it is obvious that depends linearly on the following functions of and (for each ): , , and . We use the symbol () as an abbreviated notation for these three key functions.
As was discussed above, the quantities and are not truly defined on the curve , and in particular at . Instead, each of these functions has two welldefined (but generally different) onesided limits, corresponding to approaching the worldline point from the range or . We denote these two onesided limits as and , respectively. Correspondingly, the quantities will each have two onesided limits, and , defined by
(12) 
The self force at may now be derived from either set of quantities, or , via the modesum formula Barack:1999wf ()
(13) 
where and are certain parameters (“regularization parameter”) which Refs. Barack:2001gx (); Barack:2002mh () determine analytically. Note that the two onesided limits in Eq. (13) yield the same value of .
ii.4 Frequencydomain analysis
The partial differential equation (8), which determines the functions , may be tackled in either the time domain or the frequency domain. In timedomain calculations, one directly integrates this equation numerically, using timeevolution on a twodimensional grid. In the frequencydomain method, on the other hand, one first further separates this equation into Fourier frequency modes, using
(14) 
and
(15) 
Equation (8) then reduces to the ordinary differential equation
(16) 
Since only has support on the curve , it follows that is only supported within the range .
From Eq. (7) it is evident that only depends on through and . For an eccentric geodesic is periodic, and also has its inherent periodicity. It then follows (see App. D) that is 2periodic in , namely it has a discrete spectrum of the form . Here and are the two fundamental frequencies associated with the particle’s radial and azimuthal motions, respectively. The integrals in Eqs. (14) and (15) thus reduce to summation over . In particular,
(17) 
where in principle the summation is over all integer values of .
The physicallyacceptable solutions of Eq. (16) are those satisfying the appropriate boundary conditions at both edges and , which correspond to pure outgoing waves at spatial infinity and pure incoming waves at the event horizon. The standard procedure for constructing the desired physical solution, for given , begins with the construction of a basis of two independent homogeneous solutions, and . These two homogeneous solutions satisfy the required boundary conditions at, respectively, and (note that there is no nontrivial homogeneous solution which satisfies the required boundary conditions at both edges). One then utilizes the standard Wronskianbased formula for generating inhomogeneous solutions to secondorder linear differential equations. Transforming the integration variable from to using , and recalling the bounded support of , this formula takes the form
(18)  
where
(19) 
is the Wronskian. In the regions and this formula reduces to the homogeneous solutions
(20) 
where the coefficients and are given by
(21) 
We conclude this section by explicitly writing the frequencydomain expressions for the three key functions , in the particle’s neighborhood:
(22) 
(23) 
(24) 
Since the particle resides in the range , the radial functions involved in these expressions are truly inhomogeneous [unlike the functions defined in Eq. (20), which are homogeneous].^{5}^{5}5Equation (23) should be viewed here as the Fourier decomposition of , rather than the result of a termbyterm differentiation of Eq. (22). The same applies to Eq. (24).
Iii The highfrequency problem
iii.1 Statement of the problem
The functions , whose onesided values are required for the selfforce calculation, are perfectly (onesided) smooth, even at the coincidence limit . Owing to this smoothness, one may naturally expect that there ought to be a way to calculate the required onesided quantities in the frequency domain without referring to large values. Such a method indeed exists, as we explain in the next section. However, with a straightforward application of the standard frequencydomain method, based on Eqs. (22)–(24), one finds that the Fourier series either fails to converge to the correct values (for and ) or converges very slowly (for ) as .
To demonstrate this convergence problem we consider first the Fourier sum (23) for , which is required for calculating . [Essentially the same argument applies to Eq. (24) for .] Suppose that we attempt to evaluate at a point on the worldline, with coordinates and . The values of the radial functions at are just the Fourier components of the function
(25) 
Since is discontinuous at the worldline, is discontinuous at [as well as at any other value for which ]. We therefore encounter here the Gibbs Phenomenon Gibbs (): If a function is discontinuous, its Fourier sum will fail to converge to the correct value at the discontinuity. (Away from the discontinuity the Fourier sum will converge, but rather slowly and only conditionally: The th order term will behave as .)
Consider next the convergence properties of the Fourier sum for in Eq. (22). Defining , we observe that is continuous, yet its derivative is discontinuous at . Standard Fourier theory James () then has it that the Fourier sum will indeed converge (to the correct value) at , but this convergence will be rather slow: The th term of the Fourier series is expected to behave as .
iii.2 Numerical illustration: Scalarfield monopole
It is instructive to illustrate the above problem with an explicit calculation. For this, we consider the example of the monopole mode . The spectrum in this case becomes simply (with integer ), and the Fourier sum (17) takes the form
(26) 
We hereafter use the notation , , etc. to represent the various monopole quantities. For a given orbit, the inhomogeneous mode radial functions can be computed numerically based on Eq. (18). The relevant numerical procedure is rather standard, and we relegate its description to App. C. In the following we present sample results and discuss their significance.
Figures 1 and 2 display numerical solutions for the sample orbital parameters and . (This corresponds to “semilatus rectum” and “eccentricity” , both quantities defined in App. A.) In Fig. 1 we plot the (realvalued) partial sums
(27) 
and
(28) 
as functions of at the fixed time when the particle is located at , for a sample of values. For comparison, we also plot the full monopole solution (and its derivative), which we obtain using a timedomain numerical evolution code similar to that developed by Haas in Ref. Haas:2007kz (). Evidently (and as expected), the convergence of the mode sum for both and seems very fast at and , but deteriorates significantly in the domain , where “Gibbs waves” dominate the behavior. Figure 2 illustrates the convergence properties of the partial sums and at the very location of the particle (on the same time slice as in Fig. 1). The data on the left panel suggest that the partial sum for the field converges at the particle as —in accordance with theoretical expectation James (). The results shown in the right panel of Fig. 2 demonstrate that the partial sum for the derivative fails to converge to the correct (onesided) values. They also (loosely) suggest that this partial sum in fact converges to the twoside average value of at the particle. This, indeed, would again accord with theoretical prediction James ().
iii.3 Practical implications of the highfrequency problem
The above numerical example serves to illustrate the following: From a practical point of view, it would seem very difficult to extract the correct values of the key functions at the particle’s location, based plainly on a naive summation over frequency modes as in Eqs. (22)–(24). The partial Fourier sum for the field itself would converge very slowly (as ) and its evaluation would hence be computationally expensive. Worse, the partial sums for the derivatives and would simply fail to yield the desired onesided values at the particle, even if one could sum over infinitely many modes.
Having stated the above, we should also point out that Eqs. (23) and (24) should not be deemed entirely useless for the purpose of calculating and at the particle. In principle, one could pick a point close to , yet not quite at , and calculate (say) there. The Fourier sum will converge at this point, although rather slowly. One could then pick a series of values which approach (say, from the ‘+’ side), calculate at each of these points, and then evaluate the desired sidedlimit of through extrapolation. This procedure, however, is cumbersome and is hardly likely to be computationally tractable.
An alternative implementation strategy could make use of the fact that the sums in Eqs. (23) and (24) actually converge to the twoside averages of and at the particle (recall the discussion relating to Fig. 2). These average values (along with itself) could be used to construct the average force modes , which could then be directly implemented in a “twoside averaged” version of the mode sum formula: . Although this method is likely to be by far more efficient than the extrapolation method mentioned earlier, it would still present a computational challenge, as this method, too, involves the evaluation of slowlyconverging Fourier sums.
Iv Method of extended homogeneous solutions
iv.1 Formulation of method
In this section we describe our alternative method for frequencydomain construction of the key quantities required for calculation of the self force. This method, to which we refer as the method of extended homogeneous solutions, completely avoids the highfrequency problem described above and ensures exponentiallyfast convergence of the Fourier series.
We begin by extending the definition of the homogeneous functions in Eq. (20) to the entire domain :
(29) 
We then define the two timedomain extended homogeneous solutions and to be the outcome of replacing in Eq. (22) by the homogeneous solutions or , respectively:
(30) 
We emphasize that each of the fields and is defined in the entire domain .
The convergence properties of the sum in Eq. (30) are dictated by the large asymptotic behavior of the coefficients . This highfrequency asymptotic behavior can be examined using a WKBtype analysis, which we carry out in App. D. This analysis shows that, at least within the leadingorder WKB approximation, the terms on the righthand side of Eq. (30) decay (at least) exponentially in . This exponential decay is uniform in and , throughout . Also, one naturally expects that the contribution from higherorder terms in the large WKB expansion will converge even faster than the leadingorder contribution, and hence will not affect this uniform exponential decay.
The exponential convergence of the sum (30) is extremely convenient for numerical applications, as illustrated in the next subsection. But it also has important mathematical consequences. Since the homogeneous radial functions are analytic, the uniform exponential decay of the individual terms in the above sum implies that the overall sum—namely the extended homogeneous solution —is an analytic function of and throughout .
We now argue that on each side of the curve the actual timedomain function coincides with one of these extended homogeneous solutions, namely
(31) 
For concreteness, let us present our argument explicitly referring to the first of these equalities: (i) In the domain this equality obviously holds because and coincide in that domain. (ii) As was already mentioned in Sec. II, we assume that the function is analytic throughout the range [as well as in the other range, ], because the alternative option appears unreasonable. (iii) As was just discussed above, the highfrequency analysis in App. D strongly suggests that the extended homogeneous functions are analytic throughout . (iv) Since both functions and are analytic throughout the domain [from (ii) and (iii)], and coincide at [from (i)], they must coincide throughout . (v) By continuity of both functions and , they coincide at as well. Obviously, the same line of argument applies to the second of the equalities (31) as well.
In the rest of this subsection we describe the utility of the extended homogeneous fields defined above in calculations of the selfforce via the modesum method. Recall from Eqs. (11)–(13) that this method requires as input (either of) the onesided limits of at the particle, which we now denote
(32) 
In terms of , the various components of the quantities [as defined through Eqs. (11) and (12) and used in the modesum formula (13)] are expressed directly as
(33) 
(along with ).
For concreteness, let us focus first on one of the quantities , say . By definition, the limit in Eq. (32) only samples the range . Using Eqs. (31) and (30) we may reexpress as
(34) 
Since the sum over here converges uniformly, we may interchange the limit and summation. However, the functions and are analytic, so we can now omit the limit and instead simply evaluate these functions at . The final outcome from these manipulations is stated in (the ‘+’ case of) Eq. (35) below.
The above treatment is equally applicable to , and we obtain a similar formula for constructing out of the extended homogeneous modes [the ‘’ case of Eq. (35)]. Moreover, the same treatment also applies to and . The six key quantities can all be constructed from the extended homogeneous radial functions (and their derivatives) in the form
(35) 
(36) 
(37) 
iv.2 Numerical illustration: Scalarfield monopole revisited
Let us revisit the calculation of the scalarfield monopole—this time using the method of extended homogeneous solutions. The homogeneous basis solutions are constructed numerically in just the same manner as in the standard approach (see App. B). Then, however, instead of calculating the actual inhomogeneous modes as in Sec. III, we construct the extended homogeneous solutions as they are defined in Eq. (29), with the coefficients calculated through Eq. (59) of App. C. The timedomain extended fields and their derivatives are then approximated by the (realvalued) partial sums
(38) 
(39) 
with sufficiently large .
We point out the following matters relating to the implementation of Eqs. (38) and (39): (i) In the new approach, the computation of (or ) for all and involves the (numerical) evaluation of only two integrals—the ones in Eq. (59) of App. C; In contrast, the standard approach requires the evaluation of two integrals—the ones in Eq. (56)—separately for each value of between and . (ii) The full scalar monopole is continuous at ; however, the contributions to the extended functions and from each individual mode do not match continuously along this curve. Consequently, for any finite , the partial sums for these extended functions do not match at . The amplitude of this mismatch is expected to decrease rapidly with growing , as the results below indeed demonstrate.
Figures 3–5 display numerical solutions obtained based on Eqs. (38) and (39). Our goal here is to assess the performance of the new method against the standard method, and to this end we have chosen for our numerical experiment the same orbital parameters as in Figs. 1 and 2 of Sec. III. For clarity, we only show the ‘’ and ‘’ fields in their respective relevant domains, i.e., for the former and for the latter.
Our numerical illustration serves to demonstrate the following: (i) The sum over ‘’ and ‘’ extended modes converges quickly to the correct, full solution everywhere in the respective domains and . (ii) In particular, the mismatch between the values of the ‘’ and ‘’ partial sums at the particle’s location quickly converges to zero with growing . (iii) The convergence of the extended mode sum is exponential everywhere—even in the region ; in particular, it is exponential at the very location of the particle. This applies to both the field and its derivatives. (iv) The new scheme completely circumvents the Gibbs effect which disrupts the convergence of the inhomogeneous modes in the standard approach.
V Discussion
The basic relation underpinning our new method is expressed in Eq. (31), which describes the construction of the multipole of the scalar field from extended homogeneous solutions. Eqsuation (29) and (35)–(37), in conjunction with Eqs. (33) and (13), constitute our new prescription for constructing the self force in the frequency domain. The following list highlights the advantages of this new formulation.

In the standard frequencydomain scheme, the Fourier sum over modes suffers from Gibbstype irregularities near the particle’s location. In particular, the Fourier sum fails to correctly recover the derivatives of the field’s multipoles at the particle (which are essential input in selfforce calculations). In the new method one constructs the local field’s multipoles as Fourier sums of globally homogeneous modes. These sums converge uniformly, circumventing the above complication. This is the essential and most crucial merit of the new method. [Note also that the uniform convergence (which also applies to the derivatives of the extended homogeneous mode functions) allows one to obtain the derivatives of the field’s multipoles using a termbyterm differentiation of the individual Fourier components. This again leads to Eqs. (36) and (37) above.]

The sum over the homogeneous modes converges exponentially everywhere, and even at the very location of the particle. This is extremely convenient from a practical point of view. It should be noted that, within our new scheme, it becomes “as easy” to obtain the local mode field near the particle as it is to obtain the same field in the farzone—in sharp contrast with the situation in both the standard frequencydomain method and the timedomain method.

In the standard scheme, each of the inhomogeneous modes is computed via Eq. (18). In practice, this involves the (numerical) evaluation of two integrals for each value of between and . In the new scheme, the extended homogeneous modes are obtained via Eq. (29), which requires the same two integrals (21) for all values of . Thus, remarkably, the new scheme does not only perform better mathematically—it is also much simpler to implement.

Finally, within the new scheme one is offered a convenient handle with which to monitor and control the large truncation error (i.e., the error caused by omission of the terms : The residual amount by which the partial sum (for any of the quantities ) fails to be consistent with the appropriate jump condition at is a faithful measure of the truncation error. One can therefore conveniently keep the latter below a set level by setting a threshold on the amount of residual inconsistency.
What is the application scope of the new method? In this paper we introduced the method as applied specifically for scalarfield perturbations on the Schwarzschild background. However, the basic idea is equally applicable to a wide range of other problems. It is clear from our discussion that the problematics of reconstructing the local multipole field as a sum over frequency modes has little bearing on the precise form of the underlying field equations. Rather, it is a feature of the sphericalharmonic decomposition: When we consider an individual mode, we effectively convert the physical problem of a point particle moving in an eccentric orbit into a problem of a source “shell” (which expands and contracts over time); the perturbation field is not smooth across the shell, which gives rise to Gibbstype complications when we attempt to express it as a sum over frequency modes. The same problem would occur in essentially any perturbation treatment in Schwarzschild which incorporates a spherical harmonic decomposition, such as the Teukolsky formalism (for EM and gravitational perturbations), the standard Regge–Wheeler/Zerilli/Moncrief formalism of gravitational perturbations, or the more direct Lorenzgauge formulation Barack:2005nr (). The method we propose here as a cure for the problem is directly applicable for any of these treatments. A forthcoming paper BGSprep () will report on the calculation of the monopole and dipole modes of the Lorenzgauge metric perturbation (for eccentric orbits in Schwarzschild), facilitated by the method of extended homogeneous solutions. This calculation is now being incorporated in a code which computes the total gravitational selfforce for generic orbits in Schwarzschild SBprep ().
To what extent is our new method relevant for Kerr perturbations? Here the situation is more subtle. The original modesum scheme for the self force Barack:2002mh () (which sets the main context for the current work) incorporates a decomposition in spherical harmonics even in the Kerr spacetime. This is, of course, technically possible (see the footnoted remark at the Introduction), although the resulting field equations then couple between different modes. Regardless of the latter fact, each of the modes in this decomposition would again be sourced by an expanding/contracting thin shell, the nonsmoothness of the perturbation across this shell would give rise to the Gibbs phenomenon, and our method would provide an efficient cure.
The situation is different if one tackles the Kerr problem by means of the more natural decomposition in spheroidal harmonics, which decouples the field equation in the frequency domain. Since the spheroidalharmonic functions depend on the frequency, one no longer has a strict notion of a timedomain ‘ mode’ in this case. One may (somewhat artificially) define the “spheroidalharmonic mode” by summing over all for given spheroidalharmonic numbers . However, in this case the effective geometric picture of a thin source shell would no longer apply: In the procedure of Fourier decomposition of the original point source (to obtain the source’s modes) followed by a Fourier summation over , the extra dependence of the spheroidal harmonics on will cause the reconstructed source’s mode function to deviate from a function in (for given ). Correspondingly, at a given the source’s mode function will most likely represent a “smeared” shell.
Based on the above discussion one might conclude that the highfrequency problem would not occur in the first place if spheroidal harmonics were used (as in this case there would be no type shell). We believe, however, that this may represent a false logic. The viability of the modesum approach relies crucially on the fact that the individual mode contributions in Eq. (13) are welldefined quantities. This fact is a direct consequence of the perfect onesided smoothness of the mode functions even at the limit . This smoothness, in turn, stems from the fact that for each sphericalharmonic the source term is confined to a function over a shell—and this function is distributed over the shell in a perfectly smooth manner. Note also that the functions are homogeneous timedomain solutions on both sides of this shell—even at the immediate particle’s neighborhood. In spheroidalharmonic decomposition this situation is changed, and the spheroidalharmonic mode function defined above will no longer be a homogeneous solution in the very neighborhood of the particle. It is conceivable that in the immediate particle’s neighborhood the smeared source shell will be dominated by large modes. [The sphericalharmonics decomposition is protected against this potential problem thanks to the combination of (i) the perfect offshell homogeneity, and (ii) the independence of the harmonic—and hence of the function distribution over the spherical shell—on .] Thus, in a spheroidalharmonic decomposition, the highfrequency problem may take a much more severe form: It may endanger the very existence of (namely the spheroidalharmonics analogs of ) as regular quantities, which would in turn render the (spheroidalharmonics analog of the) modesum formula (13) meaningless.
The morphology of the spheroidalharmonics smeared source shell still needs be investigated, and especially its structure near the point charge. The mode contributions may turn out to be welldefined after all, but this is far from obvious. In any case, a spheroidalharmonicsbased variant of the modesum method for Kerr has not been developed yet to the best of our knowledge. The existing Kerr modesum variant Barack:2002mh () incorporates the sphericalharmonics decomposition, and as such it exhibits the same highfrequency problem as in the Schwarzschild case. The method of extended homogeneous solutions then elegantly resolves this problem in the Kerr case as well.
Acknowledgements
LB and NS acknowledge support from PPARC/STFC through grant number PP/D001110/1.
Appendix A Eccentric geodesics in Schwarzschild
This appendix reviews the standard description of eccentric geodesics in Schwarzschild spacetime, and provides the necessary formulas for all orbital quantities (frequency, fourvelocity, etc.) needed for the numerical computations in Secs. III and IV.
As in the main text, we consider a pointlike test particle in a bound equatorial geodesic orbit around a Schwarzschild black hole with mass parameter . The radial location of the particle is bounded in the range , for some . Such geodesics constitute a twoparameter family. Each geodesic is uniquely characterized, for example, by the two “turning point” values and . An alternative parameterization (originally due to Darwin Darwin ()) employs the “semilatus rectum”, , and “eccentricity”, , both analogous to their counterparts from Keplerian celestial mechanics. The parameter pairs and are related through