Green functions of the Regge-Wheeler and Teukolsky equations in Schwarzschild spacetime

Abstract: We present a calculation of the full retarded Green functions of the Regge-Wheeler and Teukolsky equations obeyed by gravitational field perturbations of Schwarzschild spacetime. We perform the calculations for spacetime points along: (i) a timelike circular geodesic (where null-separated points are not at caustics); and (ii) a static worldline (where null-separated points are at caustics). These Green functions show a 4-fold singularity structure away from caustics, and 2-fold at caustics (similarly to the case of scalar field perturbations, which we also reproduce). Physical oscillations near the singularities appear in the gravitational case, which were not present in the scalar case. We obtain our results by developing various numerical and analytical methods.

• The paper computes the first full retarded Regge-Wheeler and BPT Green functions for Schwarzschild by combining Hadamard expansions, characteristic evolution, frequency-domain methods, and mode-sum acceleration.
• The gravitational Green functions reproduce the scalar-field singularity cycles—fourfold away from caustics and twofold at caustics—with light-crossing times explicitly identified for circular and static worldlines.
• The results reveal genuine physical oscillations near singularities, stronger for BPT than Regge-Wheeler perturbations, while unresolved quasi-local BPT terms and Green-function regularization remain important open problems.

Overview

This paper presents the first calculation of the full retarded Green functions (GFs) of the Regge-Wheeler (RW) and Bardeen-Press-Teukolsky (BPT) equations for gravitational perturbations ($s=2$ and $s=-2$, respectively) of Schwarzschild spacetime (2603.07747). Prior work had computed only the $\ell=2$ multipolar mode of the spin-2 RW GF; here the authors sum over all $\ell$-modes to obtain the complete GF, evaluated at pairs of spacetime points along two timelike worldlines: a circular geodesic at $r_0=6M$, where null-separated points ("light crossings") are not at caustics, and a static worldline, where light crossings occur at caustics. The central findings are that the gravitational GFs inherit the same global singularity structure previously established for scalar perturbations — a 4-fold cycle away from caustics and a 2-fold cycle at caustics — while additionally exhibiting genuine physical oscillations near the singularities that are absent in the scalar case.

Global singularity structure

The retarded GF diverges whenever its two arguments are connected by a null geodesic. In Schwarzschild, null geodesics can orbit the photon sphere at $r=3M$ arbitrarily many times, crossing caustics at angular separations $\gamma=0$ or $\pi$. For scalar fields it is known that after $n$ caustic crossings the leading singularity cycles through $$\text{PV}(1/\sigma)\to-\delta(\sigma)\to-\text{PV}(1/\sigma)\to\delta(\sigma)\to\cdots$$ when the field point is not itself a caustic, while at caustics the cycle is 2-fold with enhanced strength, involving asymmetric distributions of the form $s=-2$0 carrying both Dirac-delta and $s=-2$1 behavior.

For the specific worldlines considered, the paper computes the light-crossing times explicitly: on the circular geodesic at $s=-2$2, crossings occur at $s=-2$3 following the 4-fold pattern; on the static worldline, caustic crossings occur at $s=-2$4 following the 2-fold pattern. The numerical results confirm that both gravitational GFs reproduce these structures, which is expected on general grounds since all these operators share the wave operator's principal part and hence admit Hadamard forms locally.

The Regge-Wheeler Green function

The RW operator is $s=-2$5, so the direct Hadamard coefficient $s=-2$6 (the van Vleck determinant) is spin-independent, while the tail coefficient satisfies $s=-2$7, vanishing only for $s=-2$8. The full GF is obtained by the method of matched expansions between a quasi-local (QL) region and a distant-past (DP) region.

Quasi-local region: The tail term $s=-2$9 is expanded in coordinate separations as $\ell=2$0, with coefficients computed via a generalization to arbitrary RW spin of the Hadamard-WKB method, truncated at order 26 (with code made publicly available). A new variant of Padé resummation is introduced: rather than applying the approximant directly to $\ell=2$1, the singularity at the first light crossing is isolated via $\ell=2$2 and the Padé approximant is applied to the regularized $\ell=2$3. Notably, for the scalar GF on the circular geodesic this scheme with $\ell=2$4 extends well beyond the end of the normal neighbourhood — better than $\ell=2$5 — although for the static setting and for $\ell=2$6 no such extension occurs, plausibly because the caustic singularity structure involves fractional powers of $\ell=2$7.

Distant past: Two independent methods compute the $\ell=2$8-modes $\ell=2$9: (i) time-domain evolution of characteristic initial data $\ell$0 on the null surfaces through coincidence, integrated to fourth order; and (ii) inverse Fourier integrals of radial Green functions built from ingoing/upgoing homogeneous solutions, obtained via Jaffé series (detailed in an appendix) and numerical integration from the Black Hole Perturbation Toolkit. The two methods agree to more than ten significant digits at early times for $\ell$1, degrading to about eight digits by $\ell$2 due to accumulated CID error; agreement drops to roughly five to six digits at $\ell$3. Because precision loss in the Fourier modes at large $\ell$4 hinders convergence for $\ell$5, the CID method is adopted for the full RW GF, with Fourier integrals serving as validation.

The mode sum is accelerated by subtracting the analytically known direct-part modes $\ell$6 (computed from the two-dimensional reduction of Schwarzschild, using either transport equations or order-20 coordinate expansions), yielding the non-direct part whose Hadamard form is simply $\ell$7 within the normal neighbourhood. Gaussian smoothing factors suppress truncation artifacts in both frequency and $\ell$8 sums.

Results: For points on the circular geodesic at $\ell$9, the spin-2 GF (with unphysical $r_0=6M$0 modes removed, since these cannot be associated with metric perturbation multipoles) exhibits the smeared 4-fold singularity sequence at the predicted light-crossing times, with a large matching region between $r_0=6M$1 and the DP calculation. The static setting reproduces the 2-fold structure, including the characteristic left-right asymmetry about each caustic divergence. An important subtlety is identified: since $r_0=6M$2, the non-direct part should not vanish at coincidence, yet the truncated mode sum does tend toward zero there. The authors attribute this to non-commutativity of the limits $r_0=6M$3 and $r_0=6M$4 for distributional quantities — the small-$r_0=6M$5 expansion of $r_0=6M$6 is not uniform in $r_0=6M$7. This issue exists but is far less visible in the scalar case, where $r_0=6M$8.

The BPT Green function

The BPT operator contains first-derivative terms absent from the wave operator, including a purely imaginary $r_0=6M$9 coefficient inherited from the complex Newman-Penrose tetrad; consequently the BPT GF is generally complex-valued (real for equatorial points, $r=3M$0, considered here). The paper derives the Hadamard form constraints for the BPT GF: transport equations for the direct bitensor $r=3M$1 and tail bitensor $r=3M$2, with initial condition $r=3M$3 and $r=3M$4 fixed by local curvature data. However, no practical method currently exists to evaluate these biscalars — the Hadamard-WKB approach does not transfer because Eq. for $r=3M$5 lacks spherical symmetry — so the QL calculation and matching are deferred to future work.

In the DP, the BPT $r=3M$6-modes are computed entirely in the frequency domain, exploiting the Chandrasekhar transformation: homogeneous BPT solutions for $r=3M$7 are generated by a first-order differential operator acting on spin-2 RW solutions, allowing reuse of the Jaffé-series and BHPT machinery together with known transmission coefficients and Wronskian ratios. Cross-validation against direct MST-based computation of the BPT solutions shows very good agreement up to $r=3M$8, beyond which the MST implementation degrades; some high-frequency MST points fail outright and are retained in the plots to illustrate this. Large-$r=3M$9 asymptotics derived in an appendix show that, unlike the RW case where $\gamma=0$0 decays exponentially, the BPT Fourier modes decay only polynomially ($\gamma=0$1 real, $\gamma=0$2 imaginary at radial coincidence) because the complex Chandrasekhar operator mixes real and imaginary parts. The real-part integral is therefore used numerically.

Two analytic checks validate the numerical $\gamma=0$3-modes: the ringdown contribution from the fundamental quasinormal mode ($\gamma=0$4, excitation amplitude $\gamma=0$5) and the late-time Price-law tail $\gamma=0$6 derived from the branch cut, both matching the computed $\gamma=0$7 — which, to the authors' knowledge, is calculated here for the first time.

Results: The full BPT GF $\gamma=0$8 and its radial derivative display the same 4-fold (circular) and 2-fold (static) singularity patterns, with physical oscillations near the divergences that are stronger in frequency and amplitude than in the RW case. Because the required smoothing parameter must be small ($\gamma=0$9–$\pi$0), smeared divergences initially become hard to distinguish from oscillation extrema; plots at increasing $\pi$1 confirm that peaks at the light-crossing times grow while off-crossing extrema remain fixed, and an alternative smoothing factor localized around $\pi$2 resolves the ambiguity. Values within roughly $\pi$3–$\pi$4 of coincidence are flagged as unreliable due to contamination from the unresolved coincidence divergence.

Limitations and open questions

The paper is explicit about several limitations. The BPT GF is computed only in the distant past: the Hadamard biscalars $\pi$5 and $\pi$6 are characterized but not evaluated, so no QL–DP matching is achieved for BPT, and values near coincidence cannot be trusted. The coordinate expansion for $\pi$7 is not guaranteed to converge everywhere inside the normal neighbourhood, and its truncation at order 26 was constrained by computational resources. The resolution of the non-commuting-limits issue for the non-direct part at coincidence remains conjectural. On applications, the regularization of the RW or BPT GFs needed to extract a gravitational self-force via metric reconstruction is not known and is left open, as is the computation of the full QL-region GF via the "Matsubara branch cut" contour-deformation technique recently proposed in the literature.

Conclusion

This work supplies the first complete retarded Green functions for gravitational perturbations of Schwarzschild, for both the RW and BPT formulations, across the entire time domain in the RW case. It establishes that the caustic-driven 4-fold and 2-fold singularity cycles known for scalar fields persist for spin-2 perturbations, and identifies new physical oscillations near the light-crossing divergences that intensify from scalar to RW to BPT. The combination of generalized Hadamard-WKB coefficients, a modified Padé resummation, fourth-order characteristic evolution, Jaffé-series frequency-domain solutions, Chandrasekhar-transformed BPT modes, and spectroscopic checks provides a validated computational toolkit whose most immediate outstanding task is extending the BPT calculation to the quasi-local region.