Self-force via Green functions and worldline integration

1.2 MB·

by  Barry Wardell

Self-force via Green functions and worldline integration
Barry Wardell,1,2 Chad R.Galley,3 Anıl Zengino˘glu,3 Marc Casals,4 Sam R.Dolan,5 and Adrian C.Ottewill1
1School of Mathematical Sciences and Complex & Adaptive Systems Laboratory,
University College Dublin, Belfield, Dublin 4, Ireland
2Department of Astronomy, Cornell University, Ithaca, NY 14853, USA
3Theoretical Astrophysics, California Institute of Technology, Pasadena, California USA
4Department of Cosmology, Relativity and Astrophysics (ICRA),
Centro Brasileiro de Pesquisas F´ısicas, Rio de Janeiro, CEP 22290-180, Brazil.
5Consortium for Fundamental Physics, School of Mathematics and Statistics,
University of Sheffield, Hicks Building, Hounsfield Road, Sheffield S3 7RH, United Kingdom.
A compact object moving in curved spacetime interacts with its own gravitational field. This
leads to both dissipative and conservative corrections to the motion, which can be interpreted as a
self-forceactingontheobject. Theoriginalformalismdescribingthisself-forcereliedheavilyonthe
Greenfunctionofthelineardifferentialoperatorthatgovernsgravitationalperturbations. However,
because the global calculation of Green functions in non-trivial black hole spacetimes has been an
4
openproblemuntilrecently,alternativemethodswereestablishedtocalculateself-forceeffectsusing
1
sophisticatedregularizationtechniquesthatavoidthecomputationoftheglobalGreenfunction. We
0
presentamethodforcalculatingtheself-forcethatemploystheglobalGreenfunctionandistherefore
2
closely modeled after the original self-force expressions. Our quantitative method involves two

n stages: (i) numerical approximation of the retarded Green function in the background spacetime;
a (ii) evaluation of convolution integrals along the worldline of the object. This novel approach
J can be used along arbitrary worldlines, including those currently inaccessible to more established

7 computational techniques. Furthermore, it yields geometrical insight into the contributions to self-
interactionfromcurvedgeometry(back-scattering)andtrappingofnullgeodesics. Wedemonstrate

] themethodonthemotionofascalarchargeinSchwarzschildspacetime. Thistoymodelretainsthe
c
physicalhistory-dependenceoftheself-forcebutavoidsgaugeissuesandallowsustofocusonbasic
q
principles. Wecomputetheself-fieldandself-forceformanyworldlinesincludingacceleratedcircular
-
r orbits, eccentric orbits at the separatrix, and radial infall. This method, closely modeled after the
g
original formalism, provides a promising complementary approach to the self-force problem.
[

1
v I. INTRODUCTION the motion of a compact object in a curved background
6 spacetime. Thisproblemhasitsrootsinthemotionofan
0 electrically-charged particle in flat spacetime [3], though
Gravitationalwavesemittedbybinarysystemsfeatur-
5 thephysicsofthegravitationalself-forceandelectromag-
ing black holes contain a wealth of information about
1
netic radiation reaction are quite different. Astrophysi-
gravity in the strong field regime. Of particular interest
.
1 cally, the gravitational self-force is important because it
is the case of a compact body of mass m (e.g. a neutron
0 drives the inspiral of a binary with an extremely small
star or black hole) in orbit around a larger black hole
4 mass ratio (∼ 10−5–10−7). Such systems are expected
of mass M, such that m/M (cid:28) 1. In this scenario, the
1
gravitationalwavesourcesforspace-baseddetectors(e.g.,
: radiation reaction timescale is much longer than the or-
v eLISA[4,5])thatwillprovidedetailedinformationabout
bital period, and the system undergoes many cycles in
i theevolutionofgalacticmergers,generalrelativityinthe
X the vicinity of the innermost stable circular orbit (typi-
strong-field regime, and possibly the nature of dark en-
cally, ∼ M/m in the final year before coalescence). The
r
a gravitational wave signal provides a direct probe of the ergy, among others [6].
An expression for the gravitational self-force was for-
spacetime of the massive black hole.
mulated first in 1997 by Mino, Sasaki & Tanaka [7] and
In the self-force interpretation, the compact body’s
Quinn & Wald [8]. Working independently, they ob-
motionisassociatedwithaworldlinedefinedontheback-
tainedtheMiSaTaQuWa(pronouncedMee-sah-tah-kwa)
ground spacetime of the larger partner. In the test-body
equation: anexpressionfortheself-forceatfirst-orderin
limit m = 0, that worldline is a geodesic of the back-
the mass ratio (in recent years, a more rigorous foun-
ground. Forafinitemass,theworldlineisacceleratedby
dation for this formula has been established [9–11] and
a gravitational self-force as the compact body interacts
second-orderextensionsinm/M havebeenproposed[12–
withitsowngravitationalfield. Theself-forcehasadissi-
15]). The MiSaTaQuWa force can be split into “instan-
pativeandaconservativepart,whichdriveaninspiraling
taneous” and “history-dependent” terms,
worldline.
Thecalculationoftheself-forceingravitationalphysics Fα =Fα +Fα . (1)
isanimportantproblemfrombothafundamentalandan inst hist
astrophysical point of view (see reviews [1, 2]). Funda- The instantaneous terms encapsulate the interaction be-
mentally, knowing the gravitational self-force gives in- tween the source and the local spacetime geometry, and
sight into the basic physical processes that determine appear at a higher than leading order in the mass ra-
2
tio. The history-dependent term involves the worldline respectively, is the computation of the RGF along the
integration of the gradient of the retarded Green func- entire past worldline of the particle.
tion, which indicates a dependence on the past history In addition to its role in computing the self-force the
of the source’s motion. The integral extends to the in- Greenfunctionalsoplaysanessentialpartinunderstand-
finite past and is truncated just before coincidence at ing wave propagation. When the Green function to a
τ −τ(cid:48) = (cid:15) → 0+. The dependence on the history of the linear partial differential operator, such as the wave op-
source is due to gravitational perturbations that were erator, is known, any concrete solution with arbitrary
emitted in the past by the source. initial data and source can be constructed via a simple
Describing and understanding gravitational self-force convolution. Naturally, much effort has been devoted
effects is complicated due to the gauge freedom and to the study and the calculation of Green functions in
the computational burdens from the tensorial structure curved spacetimes, but a sufficiently accurate, quantita-
of gravitational perturbations. An important stepping tive description of its global behavior has remained an
stone for computing gravitational self-force has histori- open problem until recently.
cally been scalar models. These models avoid the tech- There are various difficulties regarding the calculation
nical complications of the gravitational problem but still of Green functions in curved spacetimes. In four dimen-
capturesomeessentialfeaturesofself-force,inparticular, sional flat spacetime the RGF is supported on the light-
the history-dependence. cone only (when the field is massless), whereas in curved
TheretardedGreenfunction(RGF)forthescalarfield spacetimestheRGFhasextendedsupportalsowithin the
is defined as the fundamental solution G (x,x(cid:48)) which lightcone (here “lightcone” refers to the set of points x
ret
yields causal solutions of the inhomogeneous scalar wave connectedtox(cid:48)vianullgeodesics). Further,thelightcone
equation generically self-intersects along caustics due to focusing
caused by spacetime curvature.
(cid:3) Φ(x)=−4πρ(x) (2)
x These rich features of the RGF make it difficult to
through computeglobally. Alocalapproximationisnotsufficient
toevaluatetheintegralin(6)becausetheamplitudeata
(cid:90)
Φ(x)= G (x,x(cid:48))ρ(x(cid:48))(cid:112)−g(x(cid:48))d4x(cid:48). (3) pointontheworldlinedependsontheentirepasthistory
ret
ofthesource. Earlyattempts[17–21]werefoundedupon
the Hadamard parametrix [22], which gives the RGF in
Here,G (x,x(cid:48))isadistributionwhichactsontestfunc-
ret the form
tions ρ(x(cid:48)) and also depends on the field point x, with
thepropertythatitssupportisx(cid:48) ∈J−(x),whereJ−(x) G (x,x(cid:48))=Θ (x,x(cid:48))(cid:2)U(x,x(cid:48))δ(σ)−V(x,x(cid:48))Θ(−σ)(cid:3).
ret −
is the causal past of x. This can be expressed as
(7)
(cid:3) G (x,x(cid:48))=−4πδ (x,x(cid:48)), (4)
x ret 4 Here U and V are smooth biscalars, Θ(·) is the Heavi-
side function, Θ (x,x(cid:48)) is unity when x(cid:48) is in the past of
along with appropriate boundary conditions. Here, −
δ (x,x(cid:48)) = δ (x−x(cid:48))/√−g is the four-dimensional in- x (and zero otherwise), and σ = σ(x,x(cid:48)) is the Synge
4 4
world-function (i.e., one-half of the squared distance
variantDiracdeltadistributionandg isthedeterminant
along the geodesic connecting x and x(cid:48)). The Hadamard
of the metric g .
µν parametrix is only valid in the region in which x and x(cid:48)
It was shown by Quinn [16] that the self-field of a par-
are connected by a unique geodesic (more precisely, if x
ticle with scalar charge q depends on the past history of
lies within a causal domain, commonly referred to using
the particle’s motion through
the less precise term normal neighborhood, of x(cid:48) [23]). In
Φ (zλ)=q lim (cid:90) τ−(cid:15)G (zλ,zλ(cid:48))dτ(cid:48). (5) stronglycurvedspacetimesthecontributionfromoutside
hist ret the normal neighbourhood cannot be neglected. There-
(cid:15)→0+ −∞
fore, the Hadamard parametrix is insufficient for com-
Here, zλ(τ) describes the particle’s worldline parameter- puting worldline convolutions.1 A method of matched
ized by proper time τ and a prime on an index means expansions was outlined [25, 26], in which a quasi-local
that it is associated with the worldline parameter value expansion for G , based on the Hadamard parametrix,
ret
τ(cid:48) (e.g., zα(cid:48) =zα(τ(cid:48))). Likewise, Quinn showed that the wouldbematchedontoacomplementaryexpansionvalid
self-force experienced by the particle can be written as in the more distant past.
Fα = Fα + Fα where Fα is the “instantaneous” Despite initial promise, this idea for the global eval-
inst hist inst
force, which vanishes for geodesic motion in a vacuum uation of the RGF proved difficult to implement, and
spacetime. The history-dependent term is other schemes for computing self-force came to promi-
(cid:90) τ−(cid:15) nence: themode-sumregularizationmethod[27]andthe
Fα (zλ)=q2gαβ lim ∇ G (zλ,zλ(cid:48))dτ(cid:48). (6)
hist β ret
(cid:15)→0+ −∞
Thecovariantderivativeistakenwithrespecttothefirst 1 See[24]foraproposaltocalculatetheglobalRGFusingconvolu-
argument of the RGF. The main difficulty in comput- tionsoftheHadamardparametrix,althoughdifficultforpractical
ing the self-field and the self-force through (5) and (6), applicationsinblackholespacetimes.
3
effective source method [28, 29]. A practical advantage lar time intervals, exponentially-decaying caustic echoes.
of these methods is that they work at the level of the Both the arrival intervals and the exponential decay are
(sourced) field rather than requiring global knowledge determined by the properties of null geodesics around
of RGFs. These methods require a regularization pro- thephotonsphere(theirtraveltimeandLyapunovexpo-
cedure based on a local analysis of the Green function nent). This work also provided a physical understanding
by a decomposition into the singular and regular fields of the four-fold sequence in caustic echoes: the trapping
given by Detweiler & Whiting [30]. They have yielded causes the wavefront to intersect itself at caustics as en-
impressive results [31–56], including the computation of capsulated by a Hilbert transform [71], which is equiv-
self-consistent orbits at first-order in the mass ratio [57– alent to the −π/2 phase shift observed previously. A
59]. See [1] for a recent review of numerical self-force two-foldcyclewasalsoobservedandexplainedin[24,71]
computations. wheneverthefieldpointwasintheequatorialplane(with
The steady progress in the computation of the self- thesource)atanazimuthalangleofϕ=nπfornaninte-
force compared with the lack of quantitative results on ger. Eventually, the echoes diminish below the late-time
global Green functions may have led to the notion that tail.
thecomputationoftheself-forceviatheintegrationofthe The second development is the culmination of previ-
RGF was not practicable. Nevertheless, there has been ous semi-analytical efforts [60–63] in the calculation of
recent work toward the global calculation of the RGF. the global RGF in Schwarzschild spacetime using the
The method of matched expansions for history-integral method of matched expansions [72]. The RGF was cal-
evaluation was revisited in [60], where it was shown that culated semi-analytically in the more distant past using
a distant-past expansion for the GF could be used to a Fourier-mode decomposition and contour-deformation
compute the self-force in the Nariai spacetime (a black in the complex frequency plane, which also allowed for
holetoymodel). Subsequently, stepsweretakentowards the computation of the self-force on a scalar charge in
overcoming the technical difficulties with applying the Schwarzschild spacetime.
analysis on realistic black-hole spacetimes [61–63]. In this paper, we combine a numerically computed
An interesting by-product of these investigations has global approximation to the RGF based on [71] with
been an improved understanding of the structure of the semi-analyticapproximationsatearlyandlatetimes[72]
RGF in spacetimes with caustics. While it was already to compute the self-force for arbitrary motion via the
known (e.g., [64, 65]) that singularities in the global worldline convolution integral, Eq. (6). To demonstrate
RGFoccurwhenthetwospacetimepointsareconnected thebasicprinciples,wefocusontheself-forceonascalar
via a null geodesic, the specific form of the singulari- charge in a Schwarzschild background spacetime.
tiesbeyondthenormalneighborhoodwasnotpreviously Our calculations consist of two parts. First, we con-
known within general relativity. Ori [66] noted that the struct quantitative global approximations of the Green
singular part of the RGF undergoes a transition each function in Schwarzschild spacetime by replacing the
timethewavefrontencountersacaustic,typicallycycling delta-distribution by a narrow Gaussian [71]. We aug-
through a four-fold sequence δ(σ)→1/(πσ)→−δ(σ)→ ment the numerical calculation with analytical approx-
−1/(πσ)→δ(σ) (there are some exceptions to this four- imations through quasi-local expansions at early times
fold cycle). Such four-fold structure can be understood ([20, 73, 74]) and through branch-cut integrals at late
as the wavefront picking up a phase −π/2 every time times [63, 75]. Second, we directly evaluate the convolu-
it crosses a caustic [60]. In fact, this property has first tionintegralsfortheself-field(5)andtheself-force(6)at
been discovered in optics in the late 19th century [67] the desired point along the given worldline. We discuss
andisrelatedtotheso-calledMaslovindex[68,69]. The a variety of orbital configurations, including accelerated
phasetransitionswereconfirmedviaanalyticmethodson circular orbits, eccentric geodesic orbits on the separa-
Schwarzschild spacetime [61] and black hole toy model trix, and radial infall.
space times [24, 60]. A general mathematical analysis Beyond its direct relation to the original formalism,
of the effect of caustics on wave propagation in general there are two additional motivations for computing the
relativity was given in [70]. self-force through worldline convolutions of the Green
Two recent developments along this line of research function. Thefirstmotivationisconceptual: themethod
play an essential role for this paper. The first one is allows us to answer questions such as, how does the self-
the numerical approximation of the global Green func- forcedependonthepast-history,andhowfarbackinthe
tion in Schwarzschild spacetime [71]. By replacing the historymustoneintegratetoobtainanaccurateestimate
delta-distribution source in (4) by a narrow Gaussian oftheconvolution. Thesecondmotivationiscomplemen-
andevolvingthescalarwaveequationnumerically,itwas tarity: the method is well-suited to computing the self-
possible to provide the first global approximation of the force along aperiodic (or nearly aperiodic) trajectories,
RGF [71]. The source in (4) generates “echoes” of it- such as unbound orbits, highly-eccentric or zoom-whirl
self due to trapping at the photon sphere, which were orbits, and ultra-relativistic trajectories. These motions
called “caustic echoes” in [71]. An ideal observer at null are difficult for or inaccessible to existing methods.
infinity first measures the direct signal from the delta- We believe this work challenges a perception that
distributionsource, andlater encounters, at nearlyregu- the original formulation by MiSaTaQuWa is ill-suited to
4
practical calculations. We aim to demonstrate that the versus (3+1)-dimensions] and in the role of the approxi-
new method has the potential to make self-force calcu- mated delta-distribution: as a source to the wave equa-
lations via worldline convolutions of Green functions as tion in [71] and as initial data in our approach.
routine and as accurate as (time-domain) computations
via the more established methods of mode-sum regular-
ization and effective source. 1. RGF from impulsive initial data
In Sec. II we describe methods for constructing the
RGF with numerical (IIA) and semi-analytic (IIB) ap- Given Cauchy data on a spatial hypersurface Σ, the
proaches. In Sec. III we validate our method (IIIA) KirchhoffrepresentationintermsoftheRGFcanbeused
and present a selection of results for three examples todeterminethesolutionofthehomogeneouswaveequa-
of motion: accelerated circular orbits (IIIB), eccentric tion,2
geodesics (IIIC), and radial infall (IIID). We conclude
with a discussion in Sec. IV. Throughout we choose con- (cid:3)Ψ=0 (9)
ventionswithG=c=1andmetricsignature(−+++).
at an arbitrary point x(cid:48) in the past of Σ,
II. CONSTRUCTION OF THE RGF Ψ(x(cid:48))= − 1 (cid:90) (cid:2)G (x,x(cid:48))∇µΨ(x)
4π ret
Σ
The most difficult technical step in our computation −Ψ(x)∇µGret(x,x(cid:48))(cid:3)dΣµ. (10)
of self-force is the construction of the global RGF. In √
Sec.IIAweapproximatetheGreenfunctionnumerically Here, dΣµ =nµ hd3x is the future-directed surfaceele-
byreplacingthedelta-distributionbyanarrowGaussian ment on Σ, where nµ =−αδµ0 is the future-directed unit
in the initial data using the Kirchhoff formalism. We ex- normal to the surface, α is the lapse, and h is the deter-
ploit the spherical symmetry of the background by per- minant of the induced metric. Using reciprocity, there
forming a spherical harmonic decomposition and solving is also an equivalent representation in terms of the ad-
(1+1)dimensionalwaveequationsfortheGreenfunction vanced Green function but we are interested in the RGF
usingstandardnumericalmethodsalongwithdoublehy- in this paper.
perboloidal layers. In Sec. IIB we augment the numer- We choose a time coordinate t (not necessarily
ical solution at early and late times with semi-analytic Schwarzschild time) so that the surface Σ corresponds
approximations based on the quasi-local expansion and to t=t0. Setting initial data
late-time behavior.
Ψ(x)| =Ψ(t ,x)=0 (11a)
Σ 0
n ∇µΨ(x)| =n ∇µΨ(t ,x)=−4πδ (cid:0)x,x (cid:1), (11b)
µ Σ µ 0 3 0
A. Numerical computation of the RGF
Eq. (10) reduces to
TheglobalconstructionoftheRGFin[71]usedanar- (cid:90) √
row Gaussian to approximate the 4-dimensional Dirac Ψ(x(cid:48))= G (t ,x;x(cid:48))δ (x,x ) hd3x
ret 0 3 0
delta-distribution source in (4) as, Σ
=G (t ,x ;x(cid:48)). (12)
ret 0 0
1 (cid:34) (cid:88)3 (xα−x(cid:48)α)2(cid:35)
δ (x−x(cid:48))≈ exp − . (8) Inotherwords,byevolvingthehomogeneouswaveequa-
4 (2πε¯2)2 2ε¯2
α=0 tionbackwardsintimewithappropriateimpulsiveinitial
data, we obtain the RGF with the point (t ,x ) fixed on
Here, ε¯isthewidthoftheGaussiancenteredatthebase 0 0
the hypersurface Σ for all values of x(cid:48) to the past. Thus,
pointx(cid:48). Thenumericalconstructionin[71]isperformed
with just one calculation we obtain the RGF for fixed
in (3+1) dimensions and does not rely on (or exploit)
point (t ,x ) and all possible source points z(τ(cid:48)) for tra-
any symmetries. In this paper we use the Gaussian ap- 0 0
jectories passing through (t ,x ) appearing in the self-
proximation not in the source but in initial data. We 0 0
force history-integral. Note that we have been careful
use a Kirchhoff representation of the solution of the ini-
to describe this as fundamentally evolving backwards in
tial value problem corresponding to a given Cauchy sur-
time. Thealternative—evolvingCauchydataforwards in
face in terms of the RGF and then set Gaussian initial
time—would have instead yielded G (x(cid:48);t ,x ), which
datatorecoveranapproximationtotheRGF.TheKirch- ret 0 0
would be impractical for worldline convolutions as it
hoffrepresentationcanalsobeusedin(3+1)-dimensions.
However, in a subsequent step we perform an (cid:96)-mode
decomposition to exploit the spherical symmetry of the
background, and reduce the wave equation to (1+1) di-
2 WeuseΨtodenotethefieldassociatedwithournumericalap-
mensions, rather than (3+1) dimensions. proximationoftheRGFtodistinguishitfromtheapproximated
Tosummarize,ourapproachdiffersfrom[71]inthedi- inhomogeneous solution to the wave equation, Φ, that results
mensionality of the time-domain wave equation [(1+1)- fromconvolvingΨwithagenericsource.
5
would require a separate simulation foreachpoint in the where Ψ ≡ Ψ . Likewise, we decompose the Green
(cid:96) (cid:96)m=0
convolution integral. For time-reversal-invariant back- function and its derivative by
groundssuchastheSchwarzschildspacetimethetwopro-
ceduresare,infact,equivalent,buttoemphasizethegen- 1 (cid:88)∞
G (t ,x ;x(cid:48))= (2(cid:96)+1)G (t ,r ;t(cid:48),r(cid:48))P (cosγ),
eral applicability of the method we do not exploit that ret 0 0 r r(cid:48) (cid:96) 0 0 (cid:96)
0
property here. (cid:96)=0
(19)
WecanalsocomputederivativesoftheGreenfunction
and
with respect to x(cid:48) by simply differentiating the evolved
field Ψ(x(cid:48)). However, in order to compute derivatives
∂ G(t ,x ;x(cid:48))=
with respect to the source point x, as is demanded by a r 0 0
∞ (cid:26) (cid:20) (cid:21)(cid:27)
self-forcecalculation,itisnecessarytomodifythescheme 1 (cid:88) 1
(2(cid:96)+1)P (cosγ) ∂ G (t ,r;t(cid:48),r(cid:48)) .
by setting r(cid:48) (cid:96) r r (cid:96) 0
(cid:96)=0 r=r0
(20)
Ψ(x)| =Ψ(t ,x)=0 (13a)
Σ 0
n ∇µΨ(x)| =n ∇µΨ(t ,x)=4π∂ δ (cid:0)x,x (cid:1). (13b)
µ Σ µ 0 i 3 0 Choosing the standard Schwarzschild time t and the
tortoisecoordinater∗ ≡r+2Mln(r/2M−1),theKirch-
Integrating by parts, we have
hoff formula (10) reduces to
(cid:90) √
Ψ(x(cid:48))= ∂ G (t ,x;x(cid:48))δ (x,x ) hd3x (cid:90)
i ret 0 3 0 Ψ (t(cid:48),r(cid:48))= G (t,r;t(cid:48),r(cid:48))∂ Ψ (t,r)dr∗. (21)
Σ (cid:96) (cid:96) t (cid:96)
=∂ G (t ,x ;x(cid:48)), (14) Σ
i ret 0 0
and we obtain the spatial derivatives of the Green func-
(cid:90) ∞
tion, as required. Finally, for the time derivative we set Ψ (t(cid:48),r(cid:48))= (cid:2)G (t,r;t(cid:48),r(cid:48))∂ Ψ (t,r)
(cid:96) (cid:96) t (cid:96)
−∞
−Ψ (t,r)∂ G (t,r;t(cid:48),r(cid:48))(cid:3)dr∗
(cid:0) (cid:1) (cid:96) t (cid:96)
Ψ(x)| =Ψ(t ,x)=4πδ x,x (15a)
Σ 0 3 0 (22)
n ∇µΨ(x)| =n ∇µΨ(t ,x)=0, (15b)
µ Σ µ 0
Setting initial data
to obtain
(cid:90) √ Ψ(cid:96)(t,r)|Σ =0, ∂tΨ(cid:96)(t,r)|Σ =δ(r∗−r0∗), (23)
Ψ(x(cid:48))= n ∇µG (t ,x;x(cid:48))δ (x,x ) hd3x
µ ret 0 3 0
Σ yields
=n ∇µG (t ,x ;x(cid:48)). (16)
µ ret 0 0
(cid:90) ∞
Ψ (t(cid:48),r(cid:48))= G (t ,r;t(cid:48),r(cid:48))δ(r∗−r∗)dr∗
(cid:96) (cid:96) 0 0
−∞
2. Decomposition for 1+1 dimensions
=G (t ,r ;t(cid:48),r(cid:48)). (24)
(cid:96) 0 0
In spherically symmetric spacetimes, we exploit the
Thus, by evolving the homogeneous 1+1 dimensional
spherical symmetry of the background to develop an ef-
wave equation (see Eqs. (31) and (32) below) with this
ficient scheme for solving (9) with initial data given by
initialdata,weobtainthe(cid:96)modesoftheretardedGreen
(11),(13),and(15). BydecomposingtheGreenfunction
functionwiththeretardedpointr fixedontheinitialhy-
0
into spherical harmonics and applying the Kirchhoff for-
persurface. Similarly, we can compute the radial deriva-
mula to the part depending on time and the radius, we
tive of the Green function by choosing initial data
significantly improve the efficiency of our calculations.
haCrmonosniidce(r(cid:96),thme)dmecoodmesp,osition of the field into spherical Ψ(cid:96)(t,r)|Σ =0, ∂tΨ(cid:96)(t,r)|Σ =−fr∂r(cid:2)f−1δ(r∗−r0∗)(cid:3)
(25)
∞ (cid:96)
(cid:88) (cid:88) 1 with f ≡1−2M/r. This gives
Ψ(t,r,θ,ϕ)= Ψ (t,r)Y (ϑ,ϕ). (17)
r (cid:96)m (cid:96)m
(cid:96)=0m=−(cid:96) Ψ (t(cid:48),r(cid:48))=−(cid:90) ∞ fG (t ,r;t(cid:48),r(cid:48))∂ (cid:2)f−1δ(r∗−r∗)(cid:3)dr∗
Thesphericalsymmetryoftheproblemallowsustoplace (cid:96) −∞ r (cid:96) 0 r 0
the base point x0 on the axis, in which case only the =∂ (cid:20)1G (t ,r;t(cid:48),r(cid:48))(cid:21) . (26)
m = 0 modes contribute. Then the spherical harmonic r r (cid:96) 0
decomposition reduces to r=r0
Finally, to compute the t derivative we set initial data
∞
(cid:88)1
Ψ(t,r,γ)= (2(cid:96)+1)Ψ (t,r)P (cosγ), (18)
r (cid:96) (cid:96) Ψ (t,r)| =δ(r∗−r∗), ∂ Ψ (t,r)| =0, (27)
(cid:96)=0 (cid:96) Σ 0 t (cid:96) Σ
6
which yields 10-1
ε=0.1M
(cid:90) ∞ 10-2 ε=0.2M
Ψ (t(cid:48),r(cid:48))= [∂ G (t,r;t(cid:48),r(cid:48))] δ(r∗−r∗)dr∗
(cid:96) −∞ t (cid:96) t=t0 0 Gr10-3 εε==00..46MM
=[∂tG(cid:96)(t,r0;t(cid:48),r(cid:48))]t=t0. (28) 3M 10-4
We obtain G and its t and r derivatives by solving the 10-5
(cid:96)
corresponding initial value problem. The γ derivatives 10-1
aretriviallygivenbyusing(24)in(19)anddifferentiating ‘cut=40
the Legendre polynomials with respect to γ. 10-2 ‘cut=29
Gr10-3 ‘‘cut==1172
3M 10-4 cut
3. Gaussian approximation to the Dirac delta distribution
and the smooth angular cutoff 10-5
44 45 46 47 48 49 50 51
WeapproximatetheDiracdeltadistributionintheini- ∆t/M
tial data (see Eqs. (23), (25) and (27)) on our numerical
grid by a Gaussian of finite width, ε. In the limit of zero FIG. 1. Top: The radial derivative of the Green function
width, the delta distribution is recovered, for the e = 0.5, p = 7.2 eccentric orbit case computed using
Gaussians of varying widths and a fixed smooth-sum cutoff
δ(r∗−r0∗)=εl→im0+ (2πε12)1/2e−(r∗−r0∗)2/2ε2 (29) onfea(cid:96)rcutthe=nu40ll.geAosdetshieccwriodstshingisadreocurnedas∆edt,≈th4e7.s6hMarparfeeabteutrteesr
resolved. Awayfromthenullgeodesiccrossings,theimprove-
This replacement effectively limits the shortest length mentwithdecreasingwidthismuchlesssignificant. Bottom:
scales which can be represented by the spherical har- Since a finite width Gaussian also implies an effective maxi-
mumnumberof(cid:96)modes,averysimilarplotisobtainedwhen
monicexpansionofthefield—inthiswaythereisadirect
the number of (cid:96) modes included is increased (with a fixed
correspondence between the use of a finite-width Gaus-
Gaussian width of ε=0.1M).
sian and a finite number of (cid:96) modes. Higher (cid:96) modes
oscillate more and more rapidly both in space and time.
Once the scale of these oscillations becomes compara-
to: i) the “smearing” in the angular coordinate of the
ble to the width of the Gaussian, any higher (cid:96) modes
distributional features of the Green function, and ii) the
are not faithfully resolved. The advantage of using a
replacement of the δ -source in (4) by a “narrow” Gaus-
smoothGaussian,however,isthatthiscutoffhappensin 4
sian distribution.
a smooth fashion.
The effect of varying ε and (cid:96) on the Green function
We cut off the formally infinite sum over (cid:96) in (19) and cut
can be seen in Fig. 1. The sharp features near caustic
(20) at some finite (cid:96) . A sharp cutoff in a spectral ex-
max echoes are resolved only for small ε and large (cid:96) ; but
pansion yields highly oscillatory, unphysical features in cut
away from caustic echoes the Green function is well re-
the result. The replacement of a delta distribution with
solved even for large ε and small (cid:96) .
a smooth Gaussian in our numerical scheme mitigates cut
Because the regularized self-field is a smooth func-
this somewhat; for a fixed Gaussian width ε there exists
tionofspacetime, thespectralrepresentationintermsof
a finite (cid:96) which is sufficient to eliminate the unphys-
max
spherical-harmonic modes converges exponentially with
ical oscillations (for a detailed discussion of the relation
the number of (cid:96) modes and the small-(cid:96) modes seem suf-
between ε and (cid:96) see Sec. IV A 3 of [72]). However,
max
ficient. There are, however, some cases where this ap-
forverysmallGaussianwidthsthis(cid:96) maybeunprac-
max
proximation is not very useful and large (cid:96) features be-
tically large.
come important. For example, for ultra-relativistic or-
In this work, as previously in [60, 72], we have em-
bits close to the Schwarzschild light ring at r = 3M,
ployed a smooth sum method which is very effective in
the exponentially-convergent regime is deferred to very
eliminating the high-frequency oscillations while main-
large (cid:96) [77]. In such cases, it is likely that alternative
taining the physically-relevant low (cid:96) contribution impor-
methods which capture the large-(cid:96) structure are more
tant for computing the self-force. The basic idea to to
suitable,includingasymptoticexpansionsofspecialfunc-
replacethesumin(18)(or,equivalently,in(19)and(20))
tions[24,78], WKBmethods[61], thegeometricaloptics
with a smooth cut off at large (cid:96) [76],
approximation [71], and frozen Gaussian beams [79].
Ψ(t,r,γ)=(cid:96)(cid:88)maxe−(cid:96)2/2(cid:96)2cut1r(2(cid:96)+1)Ψ(cid:96)(t,r)P(cid:96)(cosγ), souWrciethofouerrrochroinseonunrucmalecruiclaatliopnarcaommeetserfsr,omthetwdoomreilnaatendt
issues: thefactthatweonlysumoverafinitenumberof(cid:96)
(cid:96)=0
(30) modesandthatweuseaGaussianoffinitewidthε. Both
where we empirically choose (cid:96) = (cid:96) /5. The intro- approximations limit the minimum length scale (both in
cut max
duction of such a smoothing factor can be related [72] space and time) of the features that we can resolve. As
7
ε/M far from the worldline using the double layer construc-
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 tion from [82].
Wechooseacompactifyingradialcoordinateρdefined
1.40
throughr =ρ/Θ(ρ), whereΘisafunctionthatisunity
∗
1.38 in a compact domain and smoothly approaches zero on
both ends of the domain. The hyperboloidal time trans-
40 1.36
1× formation is determined such that t + r∗ = τ + ρ for
Fr1.34 r∗ →−∞, and t−r∗ =τ −ρ for r∗ →∞, which implies
2q 1.32 that level sets of the new time coordinate τ are horizon-
/
2M penetrating near the black hole and hyperboloidal near
1.30
null infinity. See [82] for details of this double hyper-
1.28 boloidal layer construction.
Substituting the spherical-harmonic decomposition of
1.26
the field, Eq. (18), into the homogeneous wave equation,
5 10 15 20 25 30 35 40 45
Eq. (9), we obtain 1+1D equations for each of the Ψ
‘ (cid:96)
cut which can be written in first order in time form,
FcuIrGv.e)2w. henExcotrmappoultaintigonthienra(cid:96)cduital(cboomttpoomnecnutrovfe)thaensdelεf-f(otrocpe ∂τΨ(cid:96) =1±Π(cid:96)H − H1∂±ρΨH(cid:96) , (31)
for the e = 0.5, p = 7.2 eccentric orbit. The convergence
(cid:104) ∂ Ψ HΠ (cid:105)
inbothcasesisapproximatelyquadratic,asindicatedbythe ∂ Π =∂ ρ (cid:96) − l
curves,whicharealeastsquaresfitofaquadraticfunctionto τ (cid:96) ρ 1±H 1±H
the data. The dashed line indicates the high-accuracy value (cid:104) 2MΘ(cid:105)(cid:104)(cid:96)((cid:96)+1) 2MΘ(cid:105)
−(Θ−ρΘ(cid:48)) 1− + Ψ ,
from [41]. ρ ρ2 ρ3 (cid:96)
(32)
showninFig.2theerrorfromtheseapproximationscon-
where a prime denotes differentiation with respect to ρ,
verges away quadratically in ε and 1/(cid:96) . By repeating
cut Π ≡(1±H)∂ Ψ +H∂ Ψ , andthefunctionH isgiven
cthane ceaxltcrualpaotiloanteftoor εa =ser0iesanodf v(cid:96)caultue=s o∞f εusainndg (cid:96)smtaanxdoanrde by(cid:96)H ≡±ddr∗(rτ∗−(cid:96)ρ)=±ρ (cid:16)(cid:96)1− Θ−Θρ2Θ(cid:48)(cid:17),wheretheunde-
Richardson extrapolation. terminedsignsarepositivenearnullinfinityandnegative
The fact that the extrapolation in both parameters near the black hole [82].
can be done reliably allows for significant improvements Weevolve(31)and(32)for200(cid:96)modesinthetimedo-
in the accuracy of the numerical results. However, doing mainusingthemethod-of-linesformulation. Wecompute
so requires multiple numerical simulations, one for each spatial derivatives using eighth-order accurate centered
value of ε. Since the focus of this paper is not on achiev- finite differencing. At the boundaries (corresponding to
ing maximal accuracy but rather on demonstrating the the black-hole horizon and null infinity) we use eighth-
feasibility of the method, we do not extrapolate in ε for order asymmetric stencils without imposing any bound-
theresultsinSec.IIIandinsteadrelyonextrapolationin ary conditions because all characteristics are purely out-
(cid:96) alone. Although this does not yield the same accuracy going in our hyperboloidal coordinates. We include a
as additionally extrapolating in ε (which gives an addi- small amount of Kreiss-Oliger [83] numerical dissipa-
tional factor of 10 improvement in accuracy in our test tion to damp out high-frequency noise. We evolve for-
cases),itdoesneverthelessimprovetheaccuracyrelative wards in time using standard fourth order Runge-Kutta
to the un-extrapolated result by a factor of 10. time integration with a fixed step size equal to half of
the spatial grid spacing. We use an equidistant spatial
grid with a resolution of ∆r∗ = 0.01M in the Regge-
4. Numerical evolution Wheeler-Zerilli coordinate region. We choose a domain
of ρ ∈ [−150M,150M] with layer interfaces located at
For the Schwarzschild spacetime, it is convenient to ±135M, so that all worldlines we consider are moving in
choosecoordinatesthatlocallymatchtheRegge-Wheeler the Regge-Wheeler-Zerilli coordinate region.
form in the vicinity of the worldline, but which also map We set Gaussian initial data of width ε = 0.1M. We
future null infinity and the future event horizon to finite evolve for a total time of t = 400M, so that the re-
values. This keeps the construction of initial data and maining portion of the history integral has a negligible
theinterpretationofresultssimple,whileimprovingcom- contribution to the self-force and its contribution to the
putational efficiency by eliminating problems with outer self-field is well approximated by the late-time asymp-
boundaries. To this end, we employ hyperboloidal com- totics described in Sec. IIB2. We extract values on the
pactification[80]intheformofalayer[81]withthestan- worldlines by simultaneously evolving the equations of
dard Schwarzschild time and tortoise coordinates in the motionusingtheosculatingorbitsframework[84]forthe
vicinity of the worldline and hyperboloidal coordinates circularandeccentricorbitsandusing (42)fortheradial
8
infall case. We interpolate the values on the worldline 2. Late-time behavior
using eighth order Legendre polynomial interpolation.
It is known since the work of Price [85, 86] that an
initial field perturbation of a Schwarzschild black hole
decays at late times with a leading power-law behav-
B. Analytical approximations to the retarded
ior: t−2(cid:96)−3 for the field mode with multipole (cid:96). This
Green function
late-time behaviour is related to the form at large ra-
dius of the potential in the ordinary differential equa-
In this section we describe the two analytical approx- tion satisfied by the radial part of the perturbation, i.e.,
imations to the RGF, one valid at early times and one the Regge-Wheeler equation for general integer spin of
valid at late times, which we use as a substitute for the the field. The form of the late-time tail can be de-
numericalsolutionintheregionswhereitcannotbeused. rived from low-frequency asymptotics along the branch
cut that the Fourier modes of the Green function pos-
sess in the complex-frequency plane [87, 88]. In [63], the
exact coefficients in the expansion at late times of the
1. Quasi-local expansion
Green function modes at arbitrary radius were given up
to the first four orders. In particular, a new behavior
In the quasi-local region, the spacetime points x and t−2(cid:96)−5log(t/M) was obtained in the third-order term,
x(cid:48) are assumed to be sufficiently close together that the thus showing a deviation from a pure power-law expan-
RGF is uniquely given by the Hadamard parametrix, sion. In the present paper we have used the method
Eq. (7). The term involving U(x,x(cid:48)) does not contribute presented in [63], which is described in detail in [75].
to the integrals in Eqs. (5) and (6) since, within those This method is based on the ‘MST formalism’ [89–91].
integrals, it has support only when x=x(cid:48), and the inte- It consists of expressing the solutions to the radial ordi-
grals exclude this point. We will therefore only concern nary differential equation as series of special functions,
ourselvesinthequasi-localregionwiththecalculationof namely ordinary hypergeometric functions and conflu-
the function V(x,x(cid:48)). ent hypergeometric functions. These series have the ad-
The proximity of x and x(cid:48) implies that an expansion vantage that they are particularly amenable to expan-
of V(x,x(cid:48)) in the separation of the points can give a sions for low frequency. We have used these series in or-
goodapproximationwithinthequasi-localregion. Refer- dertoobtainlow-frequencyasymptoticexpansionsupto
ence[74]usedaWKBmethodtoderivesuchacoordinate next-to-leading order of the Green function modes along
expansionandgaveestimatesofitsrangeofvalidity. Re- the branch cut. Finally, by integrating these expansions
ferringtotheresultstherein,wehaveV(x,x(cid:48))asapower alongthebranchcutweobtainthelate-timeasymptotics
series in (t−t(cid:48)), (1−cosγ), and (r−r(cid:48)), of the full Green function, which we will to as the late-
time tail.
∞
(cid:88)
V(x,x(cid:48))= v (r) (t−t(cid:48))2i(1−cosγ)j(r−r(cid:48))k,
ijk
III. RESULTS
i,j,k=0
(33)
where γ is the angular separation of the points and the In this section we present self-force computations for
v are computable analytic functions of r. Up to an a variety of orbits at the point where the orbits pass
ijk
overall minus sign, Eq. (33) gives the quasi-local contri- through r = 6M. One of the strengths of the Green
0
bution to the RGF. It is straightforward to take partial functionapproachisthatoncetheGreenfunctionandits
derivativesoftheseexpressionsateitherspacetimepoint derivativesareavailable, theself-forceforallorbitspass-
to obtain the derivative of the Green function. ingthroughonepointcanbecomputedbysimpleworld-
As proposed in Ref. [74], we use a Pad´e re-summation line evaluations (see Fig. 3). The entire set of results in
of Eq. (33) in order to increase the accuracy and ex- this section are obtained from just three numerical time-
tend the domain of validity of the quasi-local expansion. domain evolutions: for the Green function and its time
While not essential, this increases the region of overlap andradialderivatives(becauseofthespherical-harmonic
betweenthequasi-localandnumericalGaussiandomains. decomposition, angular derivatives act only on the Leg-
Since Pad´e re-summation is only well defined for series endrepolynomialsandthenumericaldatarequiredisthe
expansions in a single variable, it is necessary to first same as that for the undifferentiated Green function).
expand r(cid:48) and γ in a Taylor series in t−t(cid:48), using the Here, we focus on three particular families which best
equations of motion to determine the higher derivatives demonstrate the flexibility of the method:
appearing in the series coefficients. Then, with V(x,x(cid:48))
• Accelerated circular orbits, including ultra-
written as a power series in t−t(cid:48) alone, a standard di-
relativistic and static limits.
agonal Pad´e approximant of order (t−t(cid:48))26 in both the
numerator and denominator is constructed and used to • Geodesic eccentric orbits, including those near
represent the Green function in the quasi-local region. the separatrix between stable and unstable bound
9
Eccentricorbitsalongtheseparatrix component of the self-force for the eccentric orbit,
6 a higher resolution simulation with ∆r∗ = 0.005M
P
was used to estimate that the relative error from
4 numerical discretization is ∼1×10−5.
2 • Taylorseriesapproximationatearlytimes. Theac-
curacyoftheTaylorseriesapproximationdecreases
0
M as the field and source points are separated. By
(cid:144)y
repeating the calculation for a range of matching
(cid:45)2
times (between the Taylor series and the numer-
ical solution) in the region t ∈ [15M,21M], we
(cid:45)4 m
obtained relative differences in the radial self-force
for the eccentric case of ∼ 10−3. Although these
(cid:45)6
differences are comparable to those arising from a
finite cut-off in the sum over (cid:96), it turns out that
(cid:45)5 0 5 10
x(cid:144)M this is because the two issues are closely linked: by
using a larger cutoff in the (cid:96)-sum, the errors from
FIG. 3. A sample of bound geodesic orbits having different matching are also reduced. This connection also
eccentricities but passing through the same point P at r =
0 allows us to use the error from the finite cut-off in
6M. For any of these geodesics, the history-dependent part
the (cid:96)-sum as an accurate approximation for the er-
of the self-field and self-force at P can be computed from
ror from matching by choosing the matching time
fournumericaltime-domainevolutions(i.e.,G,∂ G,∂ G,and
t r
so that the error is no larger than that from the
∂ G) with P as the base point. The same is true for any
ϕ
worldline passing through P, even accelerated or unbound (cid:96)-sum.
ones.
• Asymptoticexpansionsforthelate-timetail. These
expansions become increasingly unreliable as field
geodesics. and source points are brought closer together. By
numerically evolving to t = 400M, we make sure
• Radial infall, reflecting the ability to handle un-
that the remainder of the history-integral con-
bound motion.
tributes a negligible amount to the self-force. In
Before discussing these cases, we first validate our meth- our test eccentric orbit case, the relative contribu-
ods and establish the accuracy of our results by compar- tion is ∼5×10−5.
ing against two cases for which accurate self-force values The contribution to the history integral from late
are available. timesismoresignificantforthescalarfieldbecause
theGreenfunctiondecaysmoreslowlyintimethan
its derivative. However, even in this case the mag-
A. Validation against existing results nitude of the late-time tail history integral for the
circular orbit is only ∼10−5 (compared to ∼10−2
for the entire history integral) and any error in our
Thefrequency-domainbasedmode-sumregularization
approximation is considerably smaller. We may
method is unparalleled in its ability to compute highly
therefore neglect errors from the late-time portion
accurate values for the self-force in cases where the re-
of the history integral. Unfortunately, this does
tarded field is given by a discrete spectrum involving a
not hold for generic orbits. For example, for the
small number of Fourier modes [35, 41, 92, 93]. Here, we
eccentric orbit test case, the late-time-tail-integral
considertwocaseswhichwerecomputedtohighaccuracy
error slightly dominates over other sources of error
using the methods in [41, 93] and which were considered
in the calculation; this dominance becomes more
in [72]: a circular geodesic orbit of radius r =6M, and
0
and more pronounced for orbits whose radial posi-
aneccentricorbitwitheccentricitye=0.5andsemilatus
tionvariesmoreandmoreunpredictablywithtime.
rectum p=7.2.
For this reason, when quoting errors for the self-
Our numerical calculations employ a number of ap-
field, we use the magnitude of the late-time tail
proximations:
integral as an estimate under the assumption that
• Numerical discretization of the 1+1 dimension it always provides an upper-bound (if somewhat
equation. This introduces a numerical error which over-conservative) approximation.
arisesfromerrorsinthetimeintegration,spatialfi-
nite differencing, and interpolation of values to the • Replacement of a delta distribution with a Gaus-
worldline, all of which converge away with increas- sian. This causes both a spurious pulse at early
ingresolution. Usingournumericalgridwithspac- times and a smoothing out of any sharp features.
ing ∆r∗ = 0.01M, these errors are all negligible. The initial pulse is eliminated by replacing the
Forexample,inthemostdifficultcaseoftheradial Green function at early times with its quasilocal
10
Computed value Rel. Err. Est. Err. B. Accelerated circular orbits
r M/qΦ −5.45517×10−3 6×10−5 3×10−3
ula M2/q2Ft 3.60779×10−4 4×10−4 2×10−3 We consider a particle in a circular orbit of radius r
Circ MM22//qq22FFr −15..6370846512××1100−−43 85××1100−−45 25××1100−−34 and constant angular velocity Ω, so that the azimutha0l
ϕ anglecoordinateisgivenbyϕ=Ωt. Forsuchorbits,the
c M/qΦ −7.70939×10−3 1×10−3 1×10−3
tri M2/q2F 6.65241×10−4 2×10−4 1×10−3 redshift factor is
n t
Ecce MM22//qq22FFr −71.2.38407838××1100−−43 84××1100−−45 45××1100−−34 z ≡ 1 =(cid:114)1− 2M −r2Ω2. (34)
ϕ ut r 0
0
TABLE I. Numerical results for circular and eccentric orbit
The three special cases Ω2 = {0,M/r3,(r −2M)/r3}
test cases, including estimated errors. 0 0 0
correspond to a static particle, circular geodesic, and
null orbit, respectively. Because the orbit is accelerated,
there is an additional instantaneous contribution to the
Taylorseriesexpansion. WithourchosenGaussian
self-force not present for a geodesic. This instantaneous
widthofε=0.1M,therelativeerrorcausedbythe
finite-width Gaussian is ∼1×10−3 (see Fig. 2). contribution is given by [8]
• Cutting off the infinite sum over (cid:96). Like the finite Finst = q2(g ν +u uν)Daν , (35)
width Gaussian, this has the effect of smoothing µ 3 µ µ dτ
out of any sharp features. We estimate this error
which, in the constantly accelerated circular orbit case,
by the difference between the self-force computed
has components
using (cid:96) = 40 and the value obtained by extrap-
cut
olating the curve in Fig. 2, which gives a relative q2M2(r −2M)(n2−1)(cid:20) (r −2M)(n2−1)(cid:21)
error ∼4×10−3. Finst = 0 1+ 0 ,
t 3r5z3 r z2
0 0
In our test cases, the dominant source of error therefore Finst =0, Finst =0,
r θ
comesfromthechoiceofεandthecutoffinthesumover
q2M3/2(r −2M)n(n2−1)(cid:20) M(n2−1)(cid:21)
(cid:96). Because the two sources of error are intimately con- Finst =− 0 1+ ,
nected (a choice of ε results in an effective maximum (cid:96) ϕ 3r7/2z3 r0z2
0
which can be resolved) we estimate the error in our re- (36)
sults by considering the errors in the sum over (cid:96) (for suf-
ficiently small ε). This estimate is conservative because where n ≡ Ω/Ω , and Ω = (M/r3)−1/2 is the geodesic
g g 0
it ignores accuracy improvements from extrapolating in frequency.
(cid:96) . The true accuracy of the computed self-force is up We compute the RGF and its derivatives on all con-
cut
to an order of magnitude better (see Table I). stantly accelerated circular orbits with radius r = 6M,
0
In Table I we compare the results of our numerical and plot the results in Fig. 4 as a function of the time
Green-function calculation with reference values com- ∆t/M from coincidence and of the orbital frequency Ω
puted using the frequency-domain mode-sum regulariza- relativetothegeodesicvalue. Thefigureshowsthecom-
tionmethod. Wealsogiveinternalerrorestimatesbased plex structure of the Green function that results from
on the assumption that the choice of a finite (cid:96) reason- theinterplayofwavepropagationwiththeorbitingparti-
cut
ablyreflectsthedominantsourceoferrorintheself-force cle;thelightshadedregionscorrespondtocausticechoes
and the finite integration time is the dominant source alongtheorbits. Thedottedcurvesindicatethetimesat
of error in the self-field. Our goal here is not to show which a given circular orbit crosses a point with ϕ=nπ
that the Green function method improves on the accu- for n a positive integer.
racyofexistingmethods. Frequency-domainbasedmeth- The left column in Fig. 4 shows the Green function it-
ods are the accuracy leaders for cases where they can be self(topleft)anditsderivatives(bottomthree)evaluated
used. But the Green function method provides a highly along the orbits. As the shading in the plots indicate,
flexible and complementary approach which gives good thelate-timetaildecayisslowestfortheGreenfunction,
results using modest computational resources. Further- faster for its time derivative, and fastest for the radial
more, it would not be difficult to significantly improve and azimuthal derivatives. We thus infer that late-time
on the accuracy of the results presented here through effects from back-scattering are more relevant for com-
either brute force methods (higher resolution, more (cid:96) puting the self-field than they are for the self-force. This
modes, smaller Gaussian, higher order Taylor series and isinagreementwithwhatwasobservedin[72]inthespe-
longerintegrationtimes)orthroughreadilyavailableim- cific instances of a circular geodesic and an eccentric or-
provementstothenumerics(betterhyperboloidalcoordi- bit with p=7.2 and e=0.5. Conversely, caustic echoes
nates, spectral methods for spatial derivatives, improved givenon-negligiblecontributionstotheself-forceatlater
time-integration schemes and analytic asymptotics for times than they would for the self-field because the late-
the large-(cid:96) modes [24, 78]). We leave the implementa- time tail contributes less. The contributions from the ef-
tion of such improvements for future work. fects of curved geometry (back-scattering) and trapping