Frequency-domain calculation of the self force: The high-frequency problem and its resolution

Frequency-domain calculation of the self force: The high-frequency problem and its resolution

Leor Barack, Amos Ori and Norichika Sago School of Mathematics, University of Southampton, Southampton, SO17 1BJ, United Kingdom
Department of Physics, Technion—Israel Institute of Technology, Haifa, 32000, Israel
July 15, 2019
Abstract

The mode-sum 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 spherical-harmonic -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 frequency-domain variant of this scheme the quantities are obtained by fully decomposing the particle’s self field into Fourier-harmonic 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 time-domain 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 exponentially-fast 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 scalar-field 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 self-force using either Teukolsky’s formalism, or a direct integration of the metric perturbation equations.

I Introduction

The problem of calculating the gravitational self-force 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 scalar-field self-force Q () acting on a particle endowed with a scalar charge, which proved to be a useful toy model. The electromagnetic (EM) self-force, 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 mode-sum 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 frequency-domain methods (as in, e.g., Burko:2000xx (); Detweiler:2002gi (); Keidl:2006 ()) or time-domain 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 energy-momentum, the electric four-current, or the scalar charge) as a delta-function distribution. In the frequency-domain approach one then decomposes the inhomogeneous perturbation equations into Fourier-harmonic 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, time-domain 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 time-domain approach has the advantage that one only deals with a single field for each , whereas in the frequency-domain 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 time-domain methods are winning growing popularity in recent years, frequency-domain calculations remain an appealing option for some range of orbital parameters Barton (). Also, it turns out that the non-radiative 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 frequency-domain 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 (minimally-coupled and massless) scalar-field 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:

 ϕlm(t,r)=∑ωRlmω(r)e−iωt. (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 self-force calculations we need not only but also its derivatives. For instance, to calculate the component of the self force we need to evaluate (one-sided 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 non-smooth 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 one-sided 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 self-force calculations) the multipole -modes of the physical metric perturbation 111Formally, the relevant multipoles in this case are not the standard spherical harmonics but rather the 2nd-rank tensorial harmonics. Similarly, in the discussion below regarding Teukolsky fields in Schwarzschild, the relevant multipoles are the spin-weighted 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 mode-sum 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.222Although 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 spherical-harmonic . One of the possible practical ways to construct a Fourier-spherical-harmonic mode is by solving for the separable spheroidal-harmonic modes for various , and then summing over their contributions to , as discussed in Ref. Barack:2002mh (). For given , the spherical-harmonic decomposition of the point charge will again result in a -function-type 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 Newman-Penrose 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 four-current or the energy-momentum tensor associated with the particle (a single derivative in the EM case; a second derivative in the gravitational case).333Note also that the gravitational self-force 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 self-force 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 exponentially-fast, 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 scalar-field 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 Lorenz-gauge 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 scalar-field multipoles and the scalar self-force. 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 AC 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 scalar-field 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, minimally-coupled scalar field , satisfying the field equation

 □Φ=−4πρ. (2)

Here is the scalar charge density, which takes the form of a -function along the particle’s worldline:

 ρ(x)=q∫∞−∞δ4[x−xp(τ)][−g(x)]−1/2dτ, (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

 ρ=q(r2ut)−1δ[r−rp(t)]δ[φ−φp(t)]δ[θ−π/2], (4)

where . Note that only depends on : We have , where is a constant of motion.

ii.2 Spherical-harmonics decomposition

We now separate the field equation (2) by decomposing in spherical harmonics, in the form

 Φ=∞∑l=0l∑m=−lϕlm(t,r)Ylm(θ,φ)/r. (5)

Here are the standard (complex-valued) 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:

 ρ=∞∑l=0l∑m=−l^ρlm(t,r)Ylm(θ,φ), (6)

where

 ^ρlm(t,r)=∫ρY∗lmdΩ=^clmq(r2ut)−1e−imφp(t)δ[r−rp(t)]. (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

 ∂2ϕlm∂r2∗−∂2ϕlm∂t2−Vl(r)ϕlm=−4πrf(r)^ρlm, (8)

where , is the tortoise radial coordinate, and

 Vl(r)=f(r)(2Mr3+l(l+1)r2). (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 mode-sum 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 :

 Φl(t,r,θ,φ)=l∑m=−lϕlm(t,r)Ylm(θ,φ)/r. (10)

The -component of the full-force field, , is then obtained by applying a certain linear differential operator to . [ is the same differential operator which determines the force that a smooth, “non-self” field would exert on a test charge.] In the scalar-field case we have , and therefore 444Strictly 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.

 Flμ=FμΦl=ql∑m=−l∂μ[ϕlm(t,r)Ylm(θ,φ)/r]. (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 well-defined (but generally different) one-sided limits, corresponding to approaching the worldline point from the range or . We denote these two one-sided limits as and , respectively. Correspondingly, the quantities will each have two one-sided limits, and , defined by

 F±lμ(xf)≡limx→xf±Flμ(x). (12)

The self force at may now be derived from either set of quantities, or , via the mode-sum formula Barack:1999wf ()

 Fselfμ=∞∑l=0[F±lμ(xf)∓(l+1/2)Aμ−Bμ], (13)

where and are certain parameters (“regularization parameter”) which Refs. Barack:2001gx (); Barack:2002mh () determine analytically. Note that the two one-sided limits in Eq. (13) yield the same value of .

ii.4 Frequency-domain analysis

The partial differential equation (8), which determines the functions , may be tackled in either the time domain or the frequency domain. In time-domain calculations, one directly integrates this equation numerically, using time-evolution on a two-dimensional grid. In the frequency-domain method, on the other hand, one first further separates this equation into Fourier frequency modes, using

 ϕlm(t,r)=∫dωRlmω(r)e−iωt (14)

and

 −4πrf(r)^ρlm(t,r)=∫dωZlmω(r)e−iωt. (15)

Equation (8) then reduces to the ordinary differential equation

 d2Rlmωdr2∗−[Vl(r)−ω2]Rlmω=Zlmω. (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 2-periodic 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,

 ϕlm(t,r)=∑nRlmωnm(r)e−iωnmt, (17)

where in principle the summation is over all integer values of .

The physically-acceptable 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 non-trivial homogeneous solution which satisfies the required boundary conditions at both edges). One then utilizes the standard Wronskian-based formula for generating inhomogeneous solutions to second-order linear differential equations. Transforming the integration variable from to using , and recalling the bounded support of , this formula takes the form

 Rlmω(r) = R+lmω(r)∫rrminR−lmω(r′)Zlmω(r′)Wf(r′)dr′+R−lmω(r)∫rmaxrR+lmω(r′)Zlmω(r′)Wf(r′)dr′ (18) ≡ Rinhlmω(r),

where

 W≡R−lmω(dR+lmω/dr∗)−R+lmω(dR−lmω/dr∗)=const (19)

is the Wronskian. In the regions and this formula reduces to the homogeneous solutions

 Rlmω(r)=⎧⎪ ⎪⎨⎪ ⎪⎩C−lmωR−lmω(r)≡~R−lmω(r),r≤rmin,C+lmωR+lmω(r)≡~R+lmω(r),r≥rmax, (20)

where the coefficients and are given by

 C±lmω=W−1∫rmaxrminR∓lmω(r)Zlmω(r)f(r)dr. (21)

We conclude this section by explicitly writing the frequency-domain expressions for the three key functions , in the particle’s neighborhood:

 ϕlm(t,r)=∑nRinhlmωnm(r)e−iωnmt, (22)
 ϕlm,r(t,r)=∑nddrRinhlmωnm(r)e−iωnmt, (23)
 ϕlm,t(t,r)=−i∑nωnmRinhlmωnm(r)e−iωnmt. (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].555Equation (23) should be viewed here as the Fourier decomposition of , rather than the result of a term-by-term differentiation of Eq. (22). The same applies to Eq. (24).

Iii The high-frequency problem

iii.1 Statement of the problem

The functions , whose one-sided values are required for the self-force calculation, are perfectly (one-sided) 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 one-sided 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 frequency-domain 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

 ϕflm,r(t)≡ϕlm,r(r=rf,t). (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: Scalar-field 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

 ϕ(r,t)=∞∑n=−∞Rinhn(r)e−inΩrt. (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 “semi-latus rectum” and “eccentricity” , both quantities defined in App. A.) In Fig. 1 we plot the (real-valued) partial sums

 ϕ(r,t;nmax)≡nmax∑n=−nmaxRinhn(r)e−inΩnt (27)

and

 ϕ,r(r,t;nmax)≡nmax∑n=−nmaxddrRinhn(r)e−inΩnt (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 time-domain 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 (one-sided) values. They also (loosely) suggest that this partial sum in fact converges to the two-side average value of at the particle. This, indeed, would again accord with theoretical prediction James ().

iii.3 Practical implications of the high-frequency 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 one-sided 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 sided-limit 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 two-side 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 “two-side 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 slowly-converging Fourier sums.

Iv Method of extended homogeneous solutions

iv.1 Formulation of method

In this section we describe our alternative method for frequency-domain 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 high-frequency problem described above and ensures exponentially-fast convergence of the Fourier series.

We begin by extending the definition of the homogeneous functions in Eq. (20) to the entire domain :

 ~R±lmω(r)≡C±lmωR±lmω(r),r>2M. (29)

We then define the two time-domain extended homogeneous solutions and to be the outcome of replacing in Eq. (22) by the homogeneous solutions or , respectively:

 ~ϕ±lm(t,r)≡∑n~R±lmωnm(r)e−iωnmt. (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 high-frequency asymptotic behavior can be examined using a WKB-type analysis, which we carry out in App. D. This analysis shows that, at least within the leading-order WKB approximation, the terms on the right-hand 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 higher-order terms in the large- WKB expansion will converge even faster than the leading-order 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 time-domain function coincides with one of these extended homogeneous solutions, namely

 ϕlm(t,r)=⎧⎪ ⎪⎨⎪ ⎪⎩~ϕ+lm(t,r),r≥rp(t),~ϕ−lm(t,r),r≤rp(t). (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 high-frequency 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 self-force via the mode-sum method. Recall from Eqs. (11)–(13) that this method requires as input (either of) the one-sided limits of at the particle, which we now denote

 ϕi±lm(xf)≡limx→xf±ϕilm. (32)

In terms of , the various components of the quantities [as defined through Eqs. (11) and (12) and used in the mode-sum formula (13)] are expressed directly as

 {F±lt,F±lr,F±lφ}=qrfl∑m=−l{ϕ±lm,t, ϕ±lm,r−ϕ±lm/rf, imϕ±lm}Ylm(π/2,φf) (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 re-express as

 ϕ+lm(xf)=limx→xf+∑n~R+lmωnm(r)e−iωnmt. (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

 ϕ±lm(xf)=∑n~R±lmωnm(rf)e−iωnmtf, (35)
 ϕ±lm,r(xf)=∑nddr~R±lmωnm(rf)e−iωnmtf, (36)
 ϕ±lm,t(xf)=−i∑nωnm~R±lmωnm(rf)e−iωnmtf. (37)

Equations (29) and (35)–(37), combined with Eqs. (33) and (13), constitute our new method of calculating the self force in the frequency domain. The high-frequency problem is entirely circumvented in this method, as the Fourier sum converges exponentially-fast for all functions .

iv.2 Numerical illustration: Scalar-field monopole revisited

Let us revisit the calculation of the scalar-field 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 time-domain extended fields and their derivatives are then approximated by the (real-valued) partial sums

 ~ϕ±(r,t;nmax)=nmax∑n=−nmax~R±n(r)e−inΩnt, (38)
 ~ϕ±,r(r,t;nmax)=nmax∑n=−nmaxddr~R±n(r)e−inΩnt, (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 35 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 frequency-domain scheme, the Fourier sum over -modes suffers from Gibbs-type 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 self-force 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 term-by-term 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 far-zone—in sharp contrast with the situation in both the standard frequency-domain method and the time-domain 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 scalar-field 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 spherical-harmonic 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 Gibbs-type 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 Lorenz-gauge 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 Lorenz-gauge 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 self-force 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 mode-sum 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 non-smoothness 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 spheroidal-harmonic functions depend on the frequency, one no longer has a strict notion of a time-domain ‘ mode’ in this case. One may (somewhat artificially) define the “spheroidal-harmonic -mode” by summing over all for given spheroidal-harmonic 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 high-frequency 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 mode-sum approach relies crucially on the fact that the individual mode contributions in Eq. (13) are well-defined quantities. This fact is a direct consequence of the perfect one-sided smoothness of the mode functions even at the limit . This smoothness, in turn, stems from the fact that for each spherical-harmonic 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 time-domain solutions on both sides of this shell—even at the immediate particle’s neighborhood. In spheroidal-harmonic decomposition this situation is changed, and the spheroidal-harmonic 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 spherical-harmonics decomposition is protected against this potential problem thanks to the combination of (i) the perfect off-shell homogeneity, and (ii) the independence of the harmonic—and hence of the -function distribution over the spherical shell—on .] Thus, in a spheroidal-harmonic decomposition, the high-frequency problem may take a much more severe form: It may endanger the very existence of (namely the spheroidal-harmonics analogs of ) as regular quantities, which would in turn render the (spheroidal-harmonics analog of the) mode-sum formula (13) meaningless.

The morphology of the spheroidal-harmonics smeared source shell still needs be investigated, and especially its structure near the point charge. The mode contributions may turn out to be well-defined after all, but this is far from obvious. In any case, a spheroidal-harmonics-based variant of the mode-sum method for Kerr has not been developed yet to the best of our knowledge. The existing Kerr mode-sum variant Barack:2002mh () incorporates the spherical-harmonics decomposition, and as such it exhibits the same high-frequency 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, four-velocity, 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 two-parameter 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 “semi-latus rectum”, , and “eccentricity”, , both analogous to their counterparts from Keplerian celestial mechanics. The parameter pairs and are related through

 p=2rmaxrminrmax+rmi