Effective-Source and Puncture Methods

Scrollable formula. Use Left and Right arrow keys to reveal hidden content.

Scrollable table. Use Left and Right arrow keys to reveal additional columns.

Scrollable figure. Use Left and Right arrow keys to inspect hidden content.

An effective-source calculation removes a point particle’s distributional singularity before numerical evolution. A finite-order puncture approximates the Detweiler–Whiting singular field near the worldline; subtracting it defines a residual field, and applying the wave operator analytically converts the delta distribution into an extended source of controlled differentiability. The gain is decisive in time-domain and multidimensional solvers, but the output is trustworthy when the implemented puncture, localization, boundary treatment, and numerical errors have independent checks appropriate to the claimed accuracy.

Required background. PDEs and Green Functions for Gravitation supplies distributions, Green functions, and hyperbolic boundary-value problems. Helpful background. The Einstein Equation supplies the nonlinear parent equation; Weak Metric Perturbations and Gauge Freedom supplies perturbative gauge control; Presymplectic Potentials and Currents supplies the variational interpretation of the force; Black-Hole Perturbation Theory supplies practical wave operators and radiative boundary conditions.

Green-function bridge. If $L\Phi=-4\pi\rho$ with $\rho=q\int\delta_4(x,z(\tau))d\tau$ , the equality is distributional: after multiplying by a smooth compactly supported test function and integrating, derivatives may be transferred by parts. A puncture is not a new Green function. It is a local function whose distributional image under $L$ reproduces the delta source and enough nearby singular terms that their difference from the retarded solution has a well-defined derivative on the worldline.

From a puncture to a regular numerical equation

Start with the scalar model

L\Phi_{\mathrm{ret}} \equiv(\Box-\xi R)\Phi_{\mathrm{ret}} =-4\pi\rho.

The exact Detweiler–Whiting decomposition is $\Phi_{\mathrm{ret}}=\Phi_{\mathrm S}+\Phi_{\mathrm R}$ , with $L\Phi_{\mathrm R}=0$ near the particle. In practice, construct a local covariant expansion $\Phi_{\mathrm P}^{[n]}$ such that, at spatial distance $s$ from the worldline,

\Phi_{\mathrm P}^{[n]}-\Phi_{\mathrm S}=O(s^n).

The order label must be defined by this remainder, because competing papers and codes number punctures differently. Introduce a smooth window $W$ that equals unity to sufficient order at the worldline and vanishes outside a neighborhood, and set

\Phi_{\mathrm{res}} =\Phi_{\mathrm{ret}}-W\Phi_{\mathrm P}^{[n]}.

Then

L\Phi_{\mathrm{res}} =S_{\mathrm{eff}} \equiv-4\pi\rho-L\!\left(W\Phi_{\mathrm P}^{[n]}\right).

The two terms on the right are individually distributional, but their delta pieces cancel analytically. They must never be evaluated separately as floating-point arrays. The implemented source is the already-cancelled local expression. Wardell, Vega, Thornburg, and Diener derive the generic scalar puncture and explain why direct numerical differentiation near the particle fails in Wardell, Vega, Thornburg, and Diener 2012, §§ II–III, especially eqs. (12)–(18).

For a first-order Lorenz-gauge metric perturbation, replace $L$ by

E_{\mu\nu}[\bar h] =\Box\bar h_{\mu\nu} +2R_\mu{}^\rho{}_\nu{}^\sigma\bar h_{\rho\sigma}

on a vacuum background. With a tensor puncture $\bar h_{\mu\nu}^{\mathrm P}$ ,

E_{\mu\nu}[\bar h^{\mathrm{res}}] =-16\pi G T_{\mu\nu} -E_{\mu\nu}[W\bar h^{\mathrm P}] \equiv S_{\mu\nu}^{\mathrm{eff}}.

If the puncture contains enough local orders, $\bar h^{\mathrm{res}}$ and its first derivative agree with the regular field on the worldline, while outside the window support the residual field equals the retarded field. These identities and the required conditions $W=1$ , $\nabla W=0$ , and $\nabla\nabla W=0$ on the worldline are stated in Wardell and Warburton 2015, § III B, pp. 7–8, eqs. (3.16)–(3.19).

Puncture order and finite differentiability

Applying a second-order wave operator lowers the local order by two. If $\Phi_{\mathrm P}-\Phi_{\mathrm S}=O(s^n)$ in a differentiable local asymptotic expansion, and $W=1$ near the worldline, then

S_{\mathrm{eff}}=O(s^{n-2}), \qquad \Phi_{\mathrm{res}}-\Phi_{\mathrm R}=O(s^n),

up to angular dependence and possible logarithms. A low-order source can therefore be finite but only continuous, or bounded but direction dependent. Finite differentiability can limit convergence: a nominally tenth-order stencil does not establish tenth-order accuracy when it crosses a point where the source has only one continuous derivative. Domain decomposition, local correction terms, and the observable being extracted can change the resulting convergence rate, which must be measured.

If the window is not identically one nearby, its error also matters. Because $\Phi_{\mathrm P}=O(s^{-1})$ , preserving the displayed $O(s^n)$ remainder requires $W-1=O(s^{n+1})$ with controlled derivatives. Merely requiring the first two derivatives of $W$ to vanish preserves the force limit but need not preserve an arbitrarily high puncture order.

The necessary puncture order is set by the observable. A field value requires a continuous residual; a force requires a unique first derivative; curvature invariants require two derivatives. A computation must demonstrate the relevant worldline limit rather than infer it from source boundedness. Higher local order usually improves convergence but enlarges symbolic expressions, worsens cancellation, and can make evaluation slower than mesh refinement. The best order is consequently empirical within a mathematically sufficient range.

The puncture is local and approximate. Calling it “the singular field” globally is incorrect. Changing $W$ changes the residual field away from the worldline and changes $S_{\mathrm{eff}}$ , but it must not change the regular field or a derived invariant on the worldline. This controlled nonuniqueness is a strong numerical test.

Window functions and worldtubes

There are two common localization strategies. A window method uses a smooth $W$ with compact or effectively compact support. It permits one global residual equation, but derivatives of $W$ contribute shell-like source structure. The transition width must be resolved, and the outer boundary must lie where the residual has already become the retarded field.

A worldtube method uses the puncture only inside a tube surrounding the orbit. One solves for a residual variable inside and the retarded variable outside, with matching data at the tube boundary:

\Phi_{\mathrm{ret}}^{\mathrm{out}} =\Phi_{\mathrm{res}}^{\mathrm{in}}+\Phi_{\mathrm P}, \qquad n^\alpha\nabla_\alpha\Phi_{\mathrm{ret}}^{\mathrm{out}} =n^\alpha\nabla_\alpha \left(\Phi_{\mathrm{res}}^{\mathrm{in}}+\Phi_{\mathrm P}\right).

The interface must not inject reflections or lower the scheme’s convergence order. Moving the tube boundary while holding the mesh fixed probes matching error; refining the mesh at fixed tube probes discretization error. Varying both together cannot distinguish them. Window and worldtube alternatives, including their relation in frequency-domain constructions, are compared in Wardell and Warburton 2015, § IV B, pp. 11–12, especially eq. (4.23).

Reproduction protocol and benchmark

An exact windowed source before curved spacetime

A static charge at the origin of flat space supplies a complete example of delta cancellation. Let $r=|\mathbf x|$ , take $\Phi_{\mathrm P}=q/r$ , and use the rapidly decaying window $W=\exp[-(r/w)^4]$ with width $w\gt0$ . The exact residual and its Poisson source are

\begin{aligned} \Phi_{\mathrm{res}}&=\frac{q}{r}(1-W),\\ S_{\mathrm{eff}}&=-\frac{q}{r}W''(r) =q e^{-(r/w)^4}\left(\frac{12r}{w^4}-\frac{16r^5}{w^8}\right). \end{aligned}

The residual is defined continuously by $\Phi_{\mathrm{res}}(0)=0$ and approaches $q/r$ at infinity; together with regularity these conditions fix the static solution.

To obtain the second line, use $\nabla^2(q/r)=-4\pi q\delta_3$ and $\nabla^2 f=r^{-2}(r^2f')'$ for a radial function away from the origin. The delta pieces cancel exactly. Near the charge, $\Phi_{\mathrm{res}}=qr^3/w^4+O(r^7)$ and $S_{\mathrm{eff}}=12qr/w^4+O(r^5)$ . The source is continuous but not differentiable as a Cartesian function at the origin; nevertheless the residual has a unique vanishing gradient there, giving the correct zero self-force for every $w$ .

For $q=w=1$ , the one-sided radial derivative estimate $D_h=[\Phi_{\mathrm{res}}(h)-\Phi_{\mathrm{res}}(0)]/h$ has exact error $D_h=h^2+O(h^6)$ . This Python 3 standard-library example computes it without cancellation in $1-e^{-h^4}$ :

from math import expm1

for h in (.2, .1, .05, .025):
    derivative = -expm1(-h**4) / h**2
    print(h, f"{derivative:.9e}")

Sampling distance $h/w$	Derivative error in units $q/w^2$
0.2	$3.996801706\times10^{-2}$
0.1	$9.999500017\times10^{-3}$
0.05	$2.499992188\times10^{-3}$
0.025	$6.249998779\times10^{-4}$

Halving the distance reduces the error by approximately four. This is an exact-source and extraction example, not a numerical wave evolution or a Schwarzschild benchmark. It illustrates why source regularity, derivative extraction and window choice require separate checks.

A physical scalar wave evolution

Use $G=c=1$ in this benchmark. The downloadable effective-source calculation solves the Schwarzschild scalar wave equation for a unit charge on fixed circular geodesics at $r_0/M=6$ and $10$ . It uses the established public SelfForce-1D time-domain solver at a pinned source revision (Diener, Wardell, and Warburton 2026, SelfForce-1D). This is a wave evolution with a physical point-source puncture; the flat-space calculation above remains a simpler example of exact cancellation. The $r_0=6M$ geodesic is marginally stable. Here the source orbit is held fixed; the calculation does not evolve an inspiral through it.

The small download consists of the build helper, run driver, analysis script, source probe, and case matrix, plus the shared regularization parameters and independent reference modes. Keep them in the sibling folders shown in the README. After building the pinned solver in a scratch directory, run the main $r_0=10M$ case with

python3 run_suite.py fine_r10 --case --work "$WORK" --threads 2

Here WORK is the scratch directory used by the build helper; the solver and raw time histories remain there. The README gives the analysis and complete convergence-matrix commands.

Set $M=q=1$ and choose the equatorial orbit $\theta=\pi/2$ , $\phi=\Omega t$ , with $\Omega=\sqrt{M/r_0^3}$ . Normalize the complex spherical harmonics by $\int |Y_{\ell m}|^2d\Omega=1$ , and write $\Phi=\sum_{\ell m}\psi_{\ell m}Y_{\ell m}/r$ . With $f=1-2M/r$ and $dr_*/dr=f^{-1}$ , each harmonic uses

\mathcal L_\ell =-\partial_t^2+\partial_{r_*}^2-V_\ell, \qquad V_\ell=f\left[\frac{\ell(\ell+1)}{r^2}+\frac{2M}{r^3}\right].

Inside a worldtube the evolved variable is $\psi^{\mathrm{res}}_{\ell m}=\psi^{\mathrm{ret}}_{\ell m}-\psi^{\mathrm P}_{\ell m}$ ; outside it the variable is retarded. Puncture jumps enter the numerical fluxes at the tube boundary. The worldline sits at a shared element boundary, and the observer averages the two coincident face values. It evaluates the Schwarzschild covariant component, with the harmonics evaluated on the orbit,

F_r^{\mathrm{res},\ell} =q\sum_m\operatorname{Re}\left[ Y_{\ell m}\left( \frac{\partial_{r_*}\psi^{\mathrm{res}}_{\ell m}}{r-2M} -\frac{\psi^{\mathrm{res}}_{\ell m}}{r^2} \right)\right]_{r=r_0}.

The two radial factors follow by differentiating $\Phi_{\ell m}=\psi_{\ell m}/r$ and using $\partial_r=f^{-1}\partial_{r_*}$ . Conjugate $m$ pairs count twice, while $m=0$ counts once. This quantity is a coordinate component, not a gauge-invariant gravitational observable.

The solver uses nodal discontinuous Galerkin elements and a fourth-order Runge–Kutta integrator. Hyperboloidal compactification includes the actual future horizon and future null infinity, implementing the retarded endpoint conditions without a finite-radius outer boundary. The chosen second-order puncture is the code’s named construction, with harmonic radial terms through $\Delta r^4$ and rotated azimuthal indices $m'=0,\pm1,\pm2$ . Individual harmonic sources are continuous but have a derivative kink; the residual harmonics are $C^2$ . This harmonic statement does not assign the same regularity to the unsummed three-dimensional source, nor equate the code’s order label with the earlier remainder index $n$ . Heffernan, Ottewill, Warburton, Wardell, and Diener 2018, § IV.B, p. 10, PDF describes this implementation and why placing the kink at an element boundary permits rapid convergence of the smooth one-sided approximations.

The retained source probes approach the particle from both sides. Write $L=\ell+1/2$ , as on Mode-Sum Regularization. The probes recover the analytic $B_r$ term, which is independent of $\ell$ for each circular orbit, and the $2A_{r,+}L$ derivative jump to absolute errors below $4\times10^{-17}$ and $3\times10^{-16}$ , respectively. Applying the radial wave operator independently to the puncture agrees with the supplied source to $3\times10^{-11}$ absolute error in the sampled cases; that check uses finite differences for one derivative and has its own differentiation error.

Separate the numerical errors

For the representative $r_0=10M$ , $\ell=16$ force mode, keeping the physical time step fixed at $\Delta t/M\simeq0.00211225$ gives:

DG polynomial degree	Absolute difference from the independent frequency-domain mode
8	$6.05\times10^{-8}$
10	$6.53\times10^{-10}$
14	$5.93\times10^{-13}$
18	$1.12\times10^{-12}$

The first refinements remove spatial error. The final entries reach a roughly $10^{-12}$ time-discretization and relaxation floor; they do not support a claim of indefinitely improving spatial accuracy. At fixed degree 10, reducing $\Delta t/M$ from $0.033796$ to $0.016898$ and $0.008449$ reduces the difference from the much smaller-step result from $1.61\times10^{-8}$ to $8.46\times10^{-10}$ and $5.34\times10^{-11}$ . The finer reduction is close to the factor $16$ expected for fourth-order time integration.

The full calculation retains $\ell=0,\ldots,16$ . Its raw residual modes equal the averaged retarded modes minus $B_r$ . Analytic $F_{r[2]}^\ell$ , $F_{r[4]}^\ell$ , and $F_{r[6]}^\ell$ terms from Mode-Sum Regularization are then subtracted to improve angular convergence. Their complete infinite sums vanish; this postprocessing changes neither the wave solution nor the fixed puncture order.

Every force value below denotes the dimensionless combination $M^2F_r/q^2$ in $G=c=1$ units; the calculation sets $M=q=1$ .

$r_0/M$	Evolved and regularized sum through $\ell=16$	Independent wave solution at the same cutoff
6	$1.67733512\times10^{-4}$	$1.67726288\times10^{-4}$
10	$1.37850777\times10^{-5}$	$1.37844625\times10^{-5}$

The extra digits show numerical differences, not the final accuracy claim. At these equal cutoffs the absolute differences are $7.22\times10^{-9}$ and $6.15\times10^{-10}$ : relative differences of $4.31\times10^{-5}$ and $4.46\times10^{-5}$ . The selected runs use degree 14, with $\Delta t/M\simeq0.00743582$ at $6M$ and $0.00889829$ at $10M$ ; they finish at $t/M\simeq600$ and $500$ , respectively.

Three other errors must be kept separate:

Relaxation. During the final $50M$ , the summed forces vary by $1.16\times10^{-9}$ and $2.60\times10^{-10}$ . The retained histories show the slower approach of the unseeded monopole at $6M$ .
Localization and mapping. At $10M$ , doubling the worldtube half-width from $2.04829M$ to $4.09657M$ in tortoise radius changes the fine-step result by $1.38\times10^{-9}$ . At the coarser step the change is $2.59\times10^{-8}$ , showing that tube independence improves with refinement. Moving the hyperboloidal transition changes the fine result by $1.58\times10^{-13}$ . These variations were performed at $10M$ .
Angular truncation. The independent high-precision calculation gives remaining $\ell\gt16$ contributions of $2.05\times10^{-9}$ at $6M$ and $2.01\times10^{-11}$ at $10M$ after the same analytic subtraction. They are distinct from the equal-cutoff wave errors above and are not inserted into the time-domain output.

These checks support reproduction at roughly the $10^{-4}$ relative level. The table’s smaller comparison errors are not a complete uncertainty bound: settling, localization and the uncomputed angular tail also contribute, and their signs can partly cancel. No rigorous total continuum-error bound or higher-arithmetic-precision rerun of this double-precision PDE solver is claimed.

The retained modes and numerical settings separate the angular cutoff from the numerical wave error. The independent wave solves share the known analytic regularization coefficients, so their agreement does not independently validate those shared coefficients.

The final $r_0=10M$ run and its worldtube/mapping variants initialize the static monopole analytically; higher modes still relax from a smooth turn-on. The $r_0=6M$ runs start from zero. Consequently, comparisons between seeded and zero-start full runs combine startup and discretization errors. The matched single-mode tests above isolate spatial and temporal refinement. Worldline extraction is direct at the two element faces: no off-worldline extrapolation radius is introduced.

This calculation varies degree, time step, settling time, worldtube width, and the hyperboloidal transition position. It uses one fixed puncture order and a sharp worldtube; it does not test a family of smooth windows or vary puncture order. The result establishes this implementation’s scalar benchmark, not a tensor or second-order gravitational calculation.

Published comparisons must also use the same cutoff. Warburton and Wardell’s 2014 Table I uses thirty modes, $\ell=0,\ldots,29$ , and higher-order subtraction. The independent finite sums $1.677283080352924\times10^{-4}$ and $1.3784482233828664\times10^{-5}$ reproduce its mode-sum column on rounding. Its method agreement therefore does not claim that those finite sums equal the infinite angular sums to all displayed digits. Warburton and Wardell 2014, Table I and adjacent text, p. 10, PDF

A published gravitational benchmark

A complementary published gravitational comparison concerns a point mass $\mu$ on a circular Schwarzschild geodesic at $r_0=10M$ , in Lorenz gauge with $G=c=1$ . Use ingoing horizon and outgoing infinity conditions, tensor spherical harmonics, and the dimensionless conventions

\hat h_{uu}^{\mathrm R}=\frac{M}{\mu}h_{\mu\nu}^{\mathrm R}u^\mu u^\nu, \qquad \hat F^r=\left(\frac{M}{\mu}\right)^2F^r.

The published effective-source benchmark is

\hat h_{uu}^{\mathrm R}=-0.48925802(2), \qquad \hat F^r=1.3389466(3)\times10^{-2},

with independent mode-sum values differing by at most a few parts in $10^7$ Wardell and Warburton 2015, § V B–C, pp. 13–14, Table III. This is a sharper test than checking a plot because it fixes field, component, orbit, gauge, and normalization.

The radial force here is contravariant, as in the source’s eq. (3.1). In Schwarzschild coordinates $F^r=(1-2M/r_0)F_r$ , so at $r_0=10M$ the covariant benchmark would instead be $\hat F_r=\hat F^r/0.8$ . The contraction $h_{uu}^{\mathrm R}$ includes the physical first-order mass factor before the displayed normalization. The masses in this benchmark are geometrized; its radius is $10GM_{\mathrm{phys}}$ when units are restored.

The executable package above covers the scalar calculation. The tensor benchmark remains a published comparison: reproducing it requires a compatible metric-perturbation solver, tensor puncture and force-extraction routine. A scalar source generator cannot compute it. Record the gauge, low-multipole completion, basis normalization, boundary prescription, mode cutoff and resolutions together with the raw residual modes.

A gravitational frequency-domain reproduction should follow these stages:

Generate the puncture and apply $E_{\mu\nu}$ before numerical source evaluation. Sample the cancelled $S_{\mu\nu}^{\mathrm{eff}}$ along different directions toward the particle and check its predicted powers and directional limits.
Record the source and numerical environment sufficiently to regenerate the same source and modes. Independently verify a sample of the generated expressions rather than relying on identical output alone.
Solve the homogeneous inner and outer radial systems, verify their Wronskian or constraint residuals, and impose the declared retarded boundaries before adding the compact effective-source integral.
Retain raw residual modes and low-multipole completions, and sum with the declared tail model. Increase $\ell_{\max}$ until the force and field tails both meet the chosen error target.
Repeat with doubled radial resolution or increased spectral order, then with a different window width or worldtube radius, and finally with higher arithmetic precision. Compare each change separately with the benchmark.

Convergence and stability criteria

Suppose an observable computed at grid scale $h$ has expansion $Q(h)=Q_\star+Ch^p+o(h^p)$ . Three resolutions $h$ , $h/2$ , and $h/4$ should satisfy

p_{\mathrm{obs}} =\log_2 \left| \frac{Q(h)-Q(h/2)} {Q(h/2)-Q(h/4)} \right|.

Use Richardson extrapolation only after successive refinements support a single leading power and the differences remain above roundoff and other error sources. There is no universal tolerance on $p_{\mathrm{obs}}$ : finite differentiability, extraction and interfaces must be included in the convergence prediction. Compare the extrapolated values with the benchmark using both calculations’ uncertainties. Separately vary each worldtube boundary or window width over a resolved range; the extrapolated force should agree within its error estimate. For finite-domain solvers, vary boundary locations and treatments independently. With hyperboloidal compactification, test the transition and endpoint implementation instead. Small equation residuals or convergence with one setting held fixed do not bound these other errors.

For frequency-domain calculations, the high- $\ell$ residual modes must exhibit the inverse-power falloff predicted by puncture order. Wardell and Warburton report $h_{uu}^{\mathrm R,\ell}=O(\ell^{-4})$ and $F^{r,\ell}=O(\ell^{-2})$ for their working punctures Wardell and Warburton 2015, § V C, p. 14. For time-domain calculations, halve the timestep independently of spatial refinement and verify the expected temporal order. Arithmetic precision must also be tested when cancellation or roundoff is comparable to the claimed error; it cannot be inferred from mesh convergence alone.

Failure modes and their signatures

If the delta cancellation is performed numerically, the source develops enormous opposite-sign values or NaNs near the worldline. If the puncture is one order too low, the field may converge while its derivative approaches different directional limits. An under-resolved window produces a ring of high-frequency error; a poorly matched worldtube launches reflections each time it crosses a grid feature. A boundary placed too close creates a delayed drift at the particle, whereas a bad static monopole can create an apparently stationary offset from the beginning.

Constraint violation that converges more slowly than the field signals an inconsistent tensor source or gauge equation. Dependence on window shape after continuum extrapolation means the near-worldline conditions on $W$ were not satisfied or the worldline limit was extracted incorrectly. Agreement with the benchmark at one resolution is not decisive: tail fitting, tube error, and discretization error can cancel. Finally, a first-order puncture cannot be inserted unchanged into the quadratic second-order equation; nonlinear products require a distinct construction developed on Second-Order Self-Force and Nonlinear Regularization.

Exercises

Let $\Phi_{\mathrm P}-\Phi_{\mathrm S}=a(\Omega)s^4+O(s^5)$ and let $L$ be a second-order operator with smooth coefficients. What are the expected local orders of $S_{\mathrm{eff}}$ and $\Phi_{\mathrm{res}}-\Phi_{\mathrm R}$ ? If a second-order finite-difference solver is run at spacings $h$ , $h/2$ , and $h/4$ , what ratio should appear in the convergence formula above?

Solution

Applying two derivatives to the puncture error lowers its radial order by two, so

S_{\mathrm{eff}}=O(s^2), \qquad \Phi_{\mathrm{res}}-\Phi_{\mathrm R}=O(s^4).

Angular dependence in $a(\Omega)$ can limit differentiability at the origin even though the source vanishes there, so the required derivative should still be checked along several rays. For a genuinely second-order convergent observable,

Q(h)-Q(h/2)\simeq Ch^2\left(1-\frac14\right)

and

Q(h/2)-Q(h/4)\simeq \frac14Ch^2\left(1-\frac14\right).

Their absolute ratio is $4$ , hence $p_{\mathrm{obs}}=\log_2 4=2$ .

Suppose the large- $\ell$ residual force modes behave as $a_\ell\simeq C/(\ell+1/2)^2$ . Why is the last retained mode a poor estimate of the uncomputed tail? Use $a_{16}\simeq6.1\times10^{-7}$ to estimate its scale.