Point-Source Integration and Perturbative Waveforms

A point particle supplies a distributional source to a black-hole master equation. Numerically replacing the delta distribution by an arbitrary narrow Gaussian can erase the derivative jump, slow convergence, and bias fluxes while still producing a smooth-looking waveform. A controlled calculation either integrates the distribution analytically into mode amplitudes or enforces the exact jumps at the moving worldline. This page develops that interface and its evidence; self-force regularization and the physical extreme-mass-ratio problem remain outside its numerical scope.

Required background. Ordinary Differential Equations supplies Green functions, matching, and boundary-value normalization.

Helpful background. Formulations of General Relativity fixes the geometric variables, Black-Hole Perturbation Theory derives the sourced master equations, and Gravitational Radiation and Waveforms supplies flux and waveform conventions.

Distributional sources impose jumps

For a schematic frequency-domain mode in a coordinate $x$,

$$\left(\frac{d^2}{dx^2}+\omega^2-V\right)X =G\,\delta(x-x_p)+F\,\delta'(x-x_p),$$

take $F,G$ to be coefficients evaluated at the fixed source point and define $[X]_p=X(x_p^+)-X(x_p^-)$. Distributional differentiation gives $X''=(X'')_{\rm regular}+[X']_p\delta+[X]_p\delta'$, hence $[X]_p=F$ and $[X']_p=G$. If the written source is instead $F(x)\delta'$, first use $F(x)\delta'=F(x_p)\delta'-F'(x_p)\delta$; omitting this term changes the derivative jump.

More generally, for smooth coefficients near $x_p$ and $a_p\ne0$,

$$aX''+bX'+cX=G\delta+F\delta' \quad\Longrightarrow\quad \begin{cases} a_p[X]_p=F,\\ a_p[X']_p+(b_p-a'_p)[X]_p=G. \end{cases}$$

The term $-a'_p[X]_p$ comes from multiplying $\delta'$ by the variable leading coefficient. Discontinuous operator coefficients require a separately specified weak formulation. A different radial coordinate also changes delta distributions through its Jacobian. The safe procedure is to derive jumps from the implemented equation, then test both one-sided limits. Martel and Poisson 2005, §§IV–V give the even- and odd-parity Schwarzschild source terms from which such implemented jumps must be derived.

Let $X_{\rm in}$ and $X_{\rm up}$ be horizon-ingoing and infinity-outgoing homogeneous solutions. Variation of parameters constructs the retarded field from integrals of the source divided by their Wronskian. For a delta and delta derivative, those integrals reduce to values and derivatives of the homogeneous solutions at the particle. The normalization of the delta, tortoise Jacobian, harmonic, master variable, and Fourier transform all enter the amplitude. Hopper and Evans 2010, §II.B and §§III.A–C gives explicit jump conditions and a method that avoids Gibbs contamination from a moving distribution.

For the unit-leading-coefficient stationary equation above, $W=X_{\rm in}X'_{\rm up}-X'_{\rm in}X_{\rm up}$ is constant. Direct integration against its Green function gives

$$X(x)=\begin{cases} X_{\rm in}(x)\,[G X_{\rm up}(x_p)-F X'_{\rm up}(x_p)]/W,&x\lt x_p,\\ X_{\rm up}(x)\,[G X_{\rm in}(x_p)-F X'_{\rm in}(x_p)]/W,&x\gt x_p. \end{cases}$$

The minus signs multiplying $F$ follow from $\int f(x)\delta'(x-x_p)dx=-f'(x_p)$. Subtracting the two limits recovers $[X]=F$ and $[X']=G$, independently checking the amplitude formula. Moving time-domain sources have additional velocity and acceleration terms because time derivatives also act on the moving jump; the stationary formulas must not be applied to them unchanged.

Frequency and time-domain strategies

For periodic bound motion, decompose the source into discrete frequencies and solve one radial boundary-value problem per harmonic. Sum modes only after applying consistent angular and waveform normalizations. High-frequency amplitudes can result from cancellation between large source terms, so working precision and summation order require tests. The horizon and infinity amplitudes determine fluxes; a balance relation with the orbital energy loss is a powerful independent check when the physical approximation supports it.

In the time domain, one may evolve separate homogeneous fields on either side of the worldline and communicate the analytic jumps across the moving interface. Extended homogeneous solutions reconstruct smooth one-sided fields and can converge exponentially in the harmonic sum where a direct global representation suffers Gibbs oscillations. Hopper and Evans 2010, §§IV.B–D, Fig. 4 and Tables I–II compare jumps and fluxes across implementations. Smearing is useful only when its width is treated as an additional approximation parameter and extrapolated; it is not equivalent to the distributional problem at fixed width.

Benchmark the complete normalization

Use the eccentric Schwarzschild benchmark with dimensionless semilatus rectum $p=10$, eccentricity $e=0.1$, and $\ell_{\max}=20$. Let $M_{\rm geom}=G M_{\rm phys}$ and $\mu_{\rm geom}=G\mu_{\rm phys}$; the source calculation uses $G=c=1$, so these geometric masses are the $M$ and $\mu$ printed in the source table. The target fluxes are

$$\begin{aligned} \frac{M_{\rm geom}^{2}}{\mu_{\rm geom}^{2}}\dot E^\infty &=6.31752474718\times10^{-5}, & \frac{M_{\rm geom}^{2}}{\mu_{\rm geom}^{2}}\dot E^H &=1.53365819446\times10^{-8},\\ \frac{M_{\rm geom}}{\mu_{\rm geom}^{2}}\dot L^\infty &=1.95274165241\times10^{-3}, & \frac{M_{\rm geom}}{\mu_{\rm geom}^{2}}\dot L^H &=4.48832141611\times10^{-7}. \end{aligned}$$

Hopper and Evans 2010, §IV.D and Table II report these four values and truncate both compared calculations at $\ell_{\max}=20$. Record the source coefficients, homogeneous normalizations, Wronskian, integration tolerance, frequency-mode ranges, and flux formula. Demand the source’s stated agreement at approximately the $10^{-9}$ level only after reproducing the same conventions. Also evaluate

$$\epsilon_J=\max_t\max\left( \frac{\lvert[X]_{p}-F_{\rm impl}\rvert}{X_*}, \frac{L_*\lvert[X']_{p}-G_{\rm impl}\rvert}{X_*} \right)$$

with the correct operator- and worldline-dependent expressions. Here $X_*$ is a fixed nonzero field scale, $L_*=M_{\rm geom}$ is a reference length, and the prime differentiates the declared radial coordinate. This makes both residuals dimensionless; taking an unscaled maximum of field and derivative jumps would mix units. Choose the scales before comparing resolutions, and do not divide by an instantaneous jump that may vanish.

Refine radial resolution or ODE tolerance, boundary-series order, working precision, and harmonic truncation independently. Move the numerical endpoints and matching point. Flux convergence without jump convergence is insufficient: a normalization error may cancel in one integrated observable. Conversely, a small jump residual at one time does not ensure accurate accumulated phase or mode sum.

For eccentric motion, turning points stress source interpolation and frequency sums. Demonstrate convergence over orbital phase and include modes until a controlled tail estimate is below the flux tolerance. The method’s domain is first-order perturbation theory for the declared worldline; it does not establish the worldline’s self-consistent inspiral or the validity of the point-particle approximation.

Environment and rerun record

Pin the background and orbit parameters, source and harmonic generator, ODE or evolution library, delta convention, radial coordinate, precision, boundary expansions, mode set, flux formula, and code revision. A clean interface could be point-source-wave --config eccentric-p10-e01.toml --output run/. It should emit source coefficients, jumps, homogeneous solutions, Wronskians, amplitudes, fluxes, mode-tail estimates, refinement sweeps, and environment metadata. This is the executable record required for reproduction, not a claim that such code is provided here.

Reject a result when the analytic jumps fail, the Wronskian drifts with matching radius, horizon plus infinity flux disagrees with an available balance relation, or the answer depends on Gaussian width without extrapolation. A common failure mode is a wrong delta Jacobian; others include inconsistent spherical-harmonic phase, mode cancellation at insufficient precision, and apparent convergence produced by truncating all resolutions at the same inadequate harmonic cutoff.

Diagnostic

Two implementations agree on the infinity flux to eight digits, but both replace $\delta'$ by the derivative of a Gaussian of fixed width $0.2M_{\rm geom}$. The field jump converges to zero under grid refinement although the analytic jump is nonzero. Is the flux agreement sufficient?

Solution

No. Both codes solve the same smoothed problem, and the vanished field jump proves it is not the intended distributional equation. Enforce the analytic field and derivative jumps, or vary the Gaussian width independently and extrapolate to zero while resolving it. Then compare horizon and infinity fluxes, jump residuals, and a balance relation against an implementation that integrates the delta terms analytically.

References

• Hopper, S., and Evans, C. R. (2010). “Gravitational perturbations and metric reconstruction: Method of extended homogeneous solutions applied to eccentric orbits on a Schwarzschild black hole.” Physical Review D 82, 084010. arXiv.
• Martel, K., and Poisson, E. (2005). “Gravitational perturbations of the Schwarzschild spacetime: A practical covariant and gauge-invariant formalism.” Physical Review D 71, 104003. arXiv.