Teukolsky Sources and Point-Particle Coupling

A localized body sources Kerr curvature through its stress–energy tensor, but the Teukolsky source is not $T_{ab}$ itself. It is a convention-dependent differential projection onto a principal null tetrad. The projection produces derivatives of distributions, and its mode amplitudes depend on every Fourier, harmonic, tetrad, and radial normalization chosen upstream. This page constructs that chain for a test-mass worldline and states precisely where flux calculations end and self-consistent inspiral dynamics begins.

Required background. Frames, Coframes, and Anholonomy supplies null-tetrad projections and their transformation weights.

Helpful background. Black-Hole Geometry and Mechanics supplies Kerr geodesics and horizon flux; PDEs and Green Functions for Gravitation supplies distributions and retarded solutions; Formulations of General Relativity supplies sourced linearized equations, gauge, and constraints.

Distributional stress–energy

Let a body of physical mass $\mu_p\ll M$ follow a prescribed timelike Kerr geodesic $z^a(\tau)$ with four-velocity $u^a=dz^a/d\tau$. Write $M_g=GM$ for the hole’s geometric mass, reserving $\mu_{\mathrm{NP}}$ for the spin coefficient. The Kerr functions are

$$\Sigma=r^2+a^2\cos^2\theta,\qquad \Delta=r^2-2M_gr+a^2,\qquad r_+=M_g+\sqrt{M_g^2-a^2}.$$

Use the radial page’s tortoise coordinate, $dr_*/dr=(r^2+a^2)/\Delta$. At leading order in the mass ratio, the stress–energy is

$$T^{ab}(x)= \mu_p\int u^au^b \frac{\delta^{(4)}(x-z(\tau))}{\sqrt{-g}}\,d\tau.$$

In Boyer–Lindquist coordinates, $\sqrt{-g}=\Sigma\sin\theta$. Using $t$ as the worldline parameter gives

$$T^{ab} =\frac{\mu_p u^au^b}{\Sigma\sin\theta\,u^t} \delta(r-r_p(t)) \delta(\theta-\theta_p(t)) \delta_{2\pi}(\phi-\phi_p(t)).$$

Here $\delta_{2\pi}$ is the periodic delta, normalized over one azimuthal circle. Integrating against the invariant four-volume checks the measure and the factor $1/u^t$. Covariant conservation, $\nabla_aT^{ab}=0$, holds distributionally when the worldline is geodesic. A forced trajectory without the force-carrying stress–energy violates the linearized Bianchi identity and cannot consistently source the Einstein equation.

The construction is first order in $\mu_p/M$ on a fixed Kerr background. Radiation reaction, orbital evolution, finite-size moments and a self-consistent second-order source require additional two-body dynamics, not just this source definition.

Null-tetrad projection and source preparation

For the $s=-2$ equation based on $\rho_{NP}^{-4}\delta\Psi_4$, form the Kinnersley-tetrad components

$$T_{nn}=T_{ab}n^an^b, \qquad T_{n\bar m}=T_{ab}n^a\bar m^b, \qquad T_{\bar m\bar m}=T_{ab}\bar m^a\bar m^b.$$

Use the connection coefficients defined on Kinnersley Tetrads and Newman–Penrose Curvature. Set $\rho=\rho_{NP}$, $\varpi=-\rho=(r-ia\cos\theta)^{-1}$ and $N=n^a\nabla_a=\Delta_{\mathrm{NP}}$. Bars denote complex conjugation. The complete preparation can be written in two stages:

$$\begin{aligned} U_1={}&(\bar\delta-2\bar\tau+2\alpha)T_{n\bar m} -(N+2\gamma-2\bar\gamma+\bar\mu_{\mathrm{NP}})T_{\bar m\bar m},\\ U_2={}&(N+2\gamma+2\bar\mu_{\mathrm{NP}})T_{n\bar m} -(\bar\delta-\bar\tau+2\bar\beta+2\alpha)T_{nn},\\ \mathcal Q[T]={}&(N+3\gamma-\bar\gamma+4\mu_{\mathrm{NP}} +\bar\mu_{\mathrm{NP}})U_1\\ &+(\bar\delta-\bar\tau+\bar\beta+3\alpha+4\pi)U_2,\\ T_{-2}={}&2G\varpi^{-4}\mathcal Q[T]. \end{aligned}$$

The outer derivatives act on the inner coefficients as well as on the stress distributions. This is the ordered projection from Teukolsky 1973, § II, Eq. (2.15), p. 638, PDF, translated together with the site’s curvature signs. The $G$ converts physical stress–energy to a geometrized source; the factor $\varpi^{-4}$ matches the rescaled master field. Neither can be dropped while preserving the equation $\mathcal T_{-2}[\varpi^{-4}\delta\Psi_4]=-4\pi\Sigma T_{-2}$.

For a Fourier mode with the chapter’s phase, the derivatives are explicit:

$$\begin{aligned} K&=(r^2+a^2)\omega-am, &J_+&=\partial_r+iK/\Delta,\\ L_j&=\partial_\theta+m\csc\theta-a\omega\sin\theta+j\cot\theta,\\ N&=-\frac{\Delta}{2\Sigma}J_+, &\bar\delta&=\frac{\varpi}{\sqrt2}L_0. \end{aligned}$$

This is also a practical route to a distributional source: substitute these operators, transpose their derivatives onto test functions and only then evaluate on the orbit. It avoids differentiating a sampled delta function. The coordinate operator and source preparation are given in Sasaki–Tagoshi 2003, § 2.1, Eqs. (4)–(7), pp. 7–8, PDF. Their printed Eq. (6) suppresses the $\varpi^{-4}$ needed for the rescaled field; their full radial projection, Eq. (13), contains the corresponding power. The explicit definition above fixes that ambiguity.

Under a constant spin-boost, the three raw components and the differential operators transform together. Replacing the prepared source by one component loses derivative terms and the required spin/boost weight. A position-dependent tetrad transformation also changes the connection; it does not leave this fixed Kinnersley coordinate formula unchanged.

Fourier–spheroidal projection

Normalize the polar functions by $\int_0^\pi|S_{s\ell m\omega}|^2\sin\theta\,d\theta=1$. The basis function

$$B_{s\ell m\omega}(\theta,\phi) =S_{s\ell m\omega}(\theta)e^{im\phi}$$

has norm squared $2\pi$ on the sphere. A unit-normalized full angular mode is ${}_s\mathcal Y_{\ell m\omega}=B_{s\ell m\omega}/\sqrt{2\pi}$. Use the Fourier pair

$$f(t)=\int_{-\infty}^{\infty}\frac{d\omega}{2\pi} \widetilde f(\omega)e^{-i\omega t}, \qquad \widetilde f(\omega)=\int_{-\infty}^{\infty}dt\,f(t)e^{i\omega t}.$$

Keep the same unnormalized azimuthal basis as the derivation and Reference pages. With $\mathcal T_s\psi_s=-4\pi\Sigma T_s$, define the positive radial coefficient by

$$4\pi\Sigma T_s =\int_{-\infty}^{\infty}\frac{d\omega}{2\pi} \sum_{\ell,m}T_{s\ell m\omega}(r) S_{s\ell m\omega}(\theta)e^{-i\omega t+im\phi},$$

or equivalently

$$T_{s\ell m\omega}(r) =\frac{1}{2\pi}\int dt\,d\Omega\; e^{i\omega t-im\phi}S_{s\ell m\omega}^*(\theta) \,4\pi\Sigma T_s(t,r,\theta,\phi).$$

If the unit-normalized ${}_s\mathcal Y$ is used instead, its source coefficient is $\sqrt{2\pi}\,T_{s\ell m\omega}$ and the radial field coefficient scales by the same factor. The displayed inverse projection fixes the Fourier factors directly; compare a source formula in another convention only after reconstructing its full field. In particular, Sasaki and Tagoshi’s Eq. (33) includes a $1/\sqrt{2\pi}$ not displayed in their Eq. (8) Sasaki–Tagoshi 2003, § 2.1, pp. 8–9, 12, PDF. After angular integration, a circular point source has the radial form

$$T_{s\ell m\omega}(r) =A_0(r)\delta(r-r_p) +A_1(r)\delta'(r-r_p) +A_2(r)\delta''(r-r_p)$$

with time-dependent generalizations for eccentric or inclined motion. Move derivatives between delta distributions and test functions only by distributional integration by parts. Below we derive all three coefficients in a product-derivative form, whose ordering makes that operation simple.

Circular orbit as a normalization check

Choose a prograde circular equatorial geodesic with $0\le a\lt M_g$, $r_0>r_+$ and timelike four-velocity. The source calculation does not require orbital stability, although stable orbits are the natural slow-inspiral application. Its coordinates and angular velocity are

$$r_p=r_0, \qquad \theta_p=\frac\pi2, \qquad \phi_p(t)=\Omega t+\phi_0, \qquad \Omega=\frac{\sqrt{M_g}}{r_0^{3/2}+a\sqrt{M_g}}.$$

Orbit parameters and tetrad contractions

At $(r_0,\pi/2)$ the needed metric components are

$$g_{tt}=-1+\frac{2M_g}{r_0},\qquad g_{t\phi}=-\frac{2M_ga}{r_0},\qquad g_{\phi\phi}=r_0^2+a^2+\frac{2M_ga^2}{r_0}.$$

Set $u^a=U(1,0,0,\Omega)$. Normalization and the Killing constants give

$$\begin{aligned} U&=[-g_{tt}-2\Omega g_{t\phi}-\Omega^2g_{\phi\phi}]^{-1/2},\\ E&=-U(g_{tt}+\Omega g_{t\phi}), &L_z&=U(g_{t\phi}+\Omega g_{\phi\phi}). \end{aligned}$$

The square root must be positive. Here $E=-u_t$ and $L_z=u_\phi$ are specific energy and angular momentum: the particle mass has not been absorbed into them. The radial geodesic equation $g_{tt,r}+2\Omega g_{t\phi,r}+\Omega^2g_{\phi\phi,r}=0$ checks the chosen $\Omega$.

Define $P(r)=E(r^2+a^2)-aL_z$ and $Q=L_z-aE$, holding $E,L_z,U$ at their orbit values. Equatorial contraction with the Kinnersley vectors gives

$$u\cdot n=-\frac{P}{2r^2},\qquad u\cdot\bar m=-\frac{iQ}{\sqrt2 r}.$$

Thus the three stress components are $\mu_p C_{ab}\delta(r-r_0)\delta(\theta-\pi/2) \delta_{2\pi}(\phi-\Omega t-\phi_0)/\sin\theta$, with

$$C_{nn}=\frac{P^2}{4r^6U},\qquad C_{n\bar m}=\frac{iPQ}{2\sqrt2 r^5U},\qquad C_{\bar m\bar m}=-\frac{Q^2}{2r^4U}.$$

The mixed component is imaginary; changing its sign changes the source. These follow directly from the invariant point-particle stress tensor Sasaki–Tagoshi 2003, § 2.1, Eqs. (21)–(25), pp. 10–11, PDF.

Discrete support and radial coefficients

Every projected tetrad component contains the phase $e^{-im\phi_p(t)}$. The time transform therefore gives

$$\int_{-\infty}^{\infty}dt\, e^{i\omega t-im\Omega t-im\phi_0} =2\pi e^{-im\phi_0}\delta(\omega-m\Omega).$$

An eternal circular source has support only at discrete harmonics $\omega=m\Omega$. A finite-duration or adiabatically evolving orbit broadens that distribution; replacing the delta by a narrow ordinary function is an additional approximation, not part of the exact circular source.

At a real frequency choose the polar eigenfunction $S$ real. For any other constant phase, replace each $S$ in the following source coefficients by $S^*$. Introduce the equatorial angular combinations

$$\begin{aligned} b&=a\omega-m,&S_0&=S(\pi/2),\\ H&=S'(\pi/2)+bS_0,\\ J&=S''(\pi/2)+2bS'(\pi/2)+(b^2-2)S_0. \end{aligned}$$

The angular equation at the equator gives $J=2bH-\lambda S_0$, where $\lambda=A_{-2,\ell m}+a^2\omega^2-2am\omega$. Write $k_r=K/\Delta$, distinguishing this radial coefficient from the horizon frequency $k_H$. Transposing the preparation operators against a radial test function gives the three complete coefficients

$$\begin{aligned} a_0=\frac1{2r^2U}\bigg[& -\frac{P^2}{\Delta^2}\left(J-\frac{2ia}{r}H\right) +\frac{2iPQ}{\Delta}\left(ik_r+\frac2r\right)H\\ &+Q^2S_0\left(-ik_r'-k_r^2+\frac{2ik_r}{r}\right)\bigg],\\ a_1=\frac1{r^2U}\bigg[&\frac{iPQ}{\Delta}H +Q^2S_0\left(ik_r+\frac1r\right)\bigg],\\ a_2={}&\frac{Q^2S_0}{2r^2U}, \qquad k_r'=\frac{2r\omega\Delta-K\Delta'}{\Delta^2}. \end{aligned}$$

For example, the $nn$ sector contributes the first term of $a_0$, the mixed sector contributes its second term and the first term of $a_1$, and the $\bar m\bar m$ sector supplies the remaining terms. The $2iaH/r$ term comes from differentiating $\varpi$ before setting $\theta=\pi/2$. Substituting equatorial values too early would lose it. This specializes the full separated preparation in Sasaki–Tagoshi 2003, § 2.1, Eqs. (26)–(29), p. 11, PDF to our orbit and Fourier normalization.

With $\delta_0=\delta(r-r_0)$ the result is

$$\begin{aligned} T_{-2,\ell m\omega}(r) &=2\pi\delta(\omega-m\Omega)\,\mathfrak T_{\ell m}(r),\\ \mathfrak T_{\ell m} &=G\mu_p e^{-im\phi_0}\Delta^2 \left[a_0\delta_0+\partial_r(a_1\delta_0) +\partial_r^2(a_2\delta_0)\right]. \end{aligned}$$

Each derivative acts on the whole product. This is not $\Delta^2[a_0\delta_0+a_1\delta_0'+a_2\delta_0'']$ with freely varying coefficients. Because $a_j\delta_0=a_j(r_0)\delta_0$, any smooth extension of $a_j$ with the same orbit value gives the same ordered distribution.

Outgoing amplitude and waveform

For $s=-2$ the exact Green-function Wronskian is

$$W=\Delta^{-1}\left[R^{\mathrm{in}}(R^{\mathrm{up}})' -R^{\mathrm{up}}(R^{\mathrm{in}})'\right] =2i\omega B^{\mathrm{inc}}C^{\mathrm{trans}}.$$

The first equality is constant by Abel’s identity; the second follows from the incoming $r^{-1}$ and outgoing $r^3$ branches at infinity. For $\omega=m\Omega\ne0$, away from thresholds or a zero of $W$, define the discrete outgoing amplitude by $R\sim Z^\infty_{\ell m}r^3e^{i\omega r_*}$. The radial Green weight $\Delta^{-2}$ cancels the source’s outer $\Delta^2$. Integration by parts then gives

$$Z^\infty_{\ell m} =\frac{G\mu_p e^{-im\phi_0}}{2i\omega B^{\mathrm{inc}}} \left[a_0R^{\mathrm{in}}-a_1(R^{\mathrm{in}})' +a_2(R^{\mathrm{in}})''\right]_{r_0}.$$

The factor $C^{\mathrm{trans}}$ cancels. Rescaling $R^{\mathrm{in}}$ also rescales $B^{\mathrm{inc}}$ and leaves the amplitude unchanged Sasaki–Tagoshi 2003, § 2.1, Eqs. (17)–(20), (30)–(31), pp. 9, 12, PDF. The horizon amplitude follows by replacing the numerator’s $R^{\mathrm{in}}$ by $R^{\mathrm{up}}$ and its prefactor $1/(2i\omega B^{\mathrm{inc}})$ by $B^{\mathrm{trans}}/W$. Converting that amplitude to horizon energy flux additionally requires the radial page’s Starobinsky factors.

First perform the inverse Fourier integral; do not square a frequency delta distribution. For the discrete radiative harmonics, with $u=t-r_*$,

$$\begin{aligned} \Psi_4&\sim\frac1r\sum_{\ell,m}Z^\infty_{\ell m} S_{\ell m}(\theta)e^{im\phi-i\omega u},\\ \mathcal F_E^\infty&=\sum_{\ell,m;\,\omega\ne0} \frac{|Z^\infty_{\ell m}|^2}{2G\omega^2}. \end{aligned}$$

This follows from $\Psi_4=\partial_u^2(h_+-ih_\times)/2$ and the physical flux $r^2\int\langle|\partial_u(h_+-ih_\times)|^2\rangle d\Omega/(16\pi G)$, using the full-sphere norm $2\pi$. In a unit-sphere basis, $Z_{\rm unit}=\sqrt{2\pi}Z$ gives the equivalent coefficient $|Z_{\rm unit}|^2/(4\pi G\omega^2)$. For each radiative harmonic, $\mathcal F_J=m\mathcal F_E/\omega$ Sasaki–Tagoshi 2003, § 2.1, Eqs. (33)–(36), p. 12, PDF.

A weak-field quadrupole check

An independent normalization test is available without solving a generic Kerr orbit. Set $a=0$, $v^2=M_g/r_0\ll1$ and take $\ell=m=2$, so $\omega=2v/r_0$. The horizon-regular static radial solution is $R_{\mathrm{static}}=\Delta^2$, approaching $r^4$ in the overlap $M_g\ll r\ll1/|\omega|$. The leading flat-space regular solution with that normalization is

$$\begin{aligned} R_{\mathrm{reg}} &=r^4e^{-i\omega r}\,{}_1F_1(5;6;2i\omega r)\\ &=\frac5{(2i\omega)^5r} \left[e^{i\omega r}P_4(2i\omega r)-24e^{-i\omega r}\right],\\ P_4(z)&=z^4-4z^3+12z^2-24z+24. \end{aligned}$$

Direct differentiation verifies $r^2R''-2rR'+(\omega^2r^2-4i\omega r-4)R=0$ and $R_{\mathrm{reg}}/r^4\to1$ as $r\to0$. Its incoming coefficient is $B^{\mathrm{inc}}=15i/(4\omega^5)$, hence $1/(2i\omega B^{\mathrm{inc}})=-2\omega^4/15$. The normalized polar harmonic is $S_{22}=\sqrt{5/32}(1+\cos\theta)^2$, giving $S_0=\sqrt{5/32}$, $H=-4S_0$ and $J=12S_0$. The leading source term is $a_0=-6S_0/r_0^2$; the $a_1,a_2$ contributions are higher order in $v$. The amplitude is therefore

$$Z^\infty_{22}=\frac{16}{\sqrt{10}} \frac{G\mu_p}{r_0^2}v^4e^{-2i\phi_0}[1+o(1)].$$

The $m=-2$ harmonic has the same leading magnitude. Including both gives

$$\mathcal F_E^\infty=\frac{32}{5}\frac{G\mu_p^2v^6}{r_0^2} =\frac{32}{5}\frac{G^4\mu_p^2M^3}{r_0^5} \quad\text{at leading order}.$$

Independently, the Newtonian orbit $\boldsymbol x=r_0(\cos\Omega t,\sin\Omega t,0)$ has STF quadrupole $I_{ij}=\mu_p(x_ix_j-\delta_{ij}r_0^2/3)$. Its directly differentiated quadrupole power is

$$\frac G5\dddot I_{ij}\dddot I_{ij} =\frac{32}{5}G\mu_p^2r_0^4\Omega^6,$$

which agrees. This tests the angular norm, Fourier factor, Green-function normalization and $G$ restoration together. Squared power cannot detect an overall sign, so compare the waveform too. At observer azimuth zero, let $c=\cos\theta$, set $\phi_0=0$ and denote the distant radius by $r_{\mathrm{obs}}$. Both the two curvature modes and the direct transverse STF contraction give

$$h_+-ih_\times=-\frac{G\mu_pv^2}{r_{\mathrm{obs}}} \left[(1+c)^2e^{-2i\Omega u}+(1-c)^2e^{2i\Omega u}\right].$$

For the direct check, use $G\bar e^i\bar e^j\ddot I_{ij}/r_{\mathrm{obs}}$ with $\bar{\boldsymbol e}=(c,-i,-\sin\theta)$. This screen vector is transverse and bilinearly null, so it already selects the transverse-traceless polarization. The negative sign agrees with $h_+-ih_\times=-2\Psi_4/\omega^2$ for each harmonic. This is a leading weak-field test, not a full Kerr absorption calculation or a local self-force calculation.

Flux balance and regularization boundary

The prescribed geodesic itself has constant $E$ and $L_z$: it does not already undergo radiation reaction. Its first-order field supplies quadratic radiative fluxes. In an additional adiabatic, dissipative approximation, those periodic-orbit fluxes determine the slow orbital loss,

$$\left\langle\frac{dE_{\mathrm{orb}}}{dt}\right\rangle =-\left( \mathcal F_E^\infty+\mathcal F_E^H \right),$$

with $E_{\mathrm{orb}}=\mu_p E$ to leading order and a parallel angular-momentum relation. Horizon flux may be negative in the superradiant regime. The explicit infinity and horizon conversions are given in Sasaki–Tagoshi 2003, § 2.1, Eqs. (33)–(36), p. 12, PDF and Sasaki–Tagoshi 2003, § 5.6, Eqs. (175)–(176), p. 41, PDF.

This balance is an independent end-to-end check: tetrad projection, Fourier support, radial normalization, asymptotic conversion, and orbital sign must agree. It is not a local self-force calculation. The retarded field diverges on the worldline, and a local force or regular metric requires singular-field subtraction and a specified gauge. Fluxes at the horizon and infinity can be finite without solving that local regularization problem.

Common pitfalls

Projecting before fixing normalization. A source coefficient changes when the harmonic, Fourier transform, master field, or radial solution is rescaled. Record all four before comparing amplitudes.

Dropping derivative distributions. The source-preparation operators differentiate the point-particle delta functions. Omitting $\delta'$ or $\delta''$ terms generally violates the Teukolsky equation.

Equating balance with inspiral. An orbit-averaged flux supplies dissipative evolution only after an adiabatic two-body approximation. It does not provide conservative self-force, transient resonances, finite-size effects, or second-order dynamics.

Exercises

For the circular equatorial worldline above, derive its frequency support and evaluate

$$I=\int_{r_+}^{\infty}F(r) \left[A_0\delta(r-r_0)+A_1\delta'(r-r_0)+A_2\delta''(r-r_0)\right]dr$$

when $r_0$ lies strictly inside the integration interval. State two checks that would detect a source-normalization error.

Solution

The angular phase in the projection is $e^{-im\phi_p}=e^{-im\Omega t-im\phi_0}$. Multiplication by the forward Fourier projection phase $e^{i\omega t}$ and integration over time gives $2\pi e^{-im\phi_0}\delta(\omega-m\Omega)$, so only $\omega=m\Omega$ contributes.

Distributional differentiation is defined by $\int F\delta'=-F'$ and $\int F\delta''=F''$. Hence

$$I=A_0F(r_0)-A_1F'(r_0)+A_2F''(r_0),$$

if the $A_j$ shown are constants. For smooth functions multiplying the distributions, the general answer is

$$I=\left[A_0F-(A_1F)'+(A_2F)''\right]_{r_0}.$$

This differs from pairing derivatives of the whole products $A_j\delta$; one must preserve the chosen order before moving any coefficient outside.

First, reconstruct the four-dimensional field and verify that rescaling the unit-normalized harmonic by a constant while inversely rescaling its coefficient leaves the field and flux unchanged. Second, compute the orbit-averaged energy loss independently and verify $\langle\dot E_{\rm orb}\rangle+\mathcal F_E^\infty+\mathcal F_E^H=0$ with the same orientation. A missing $2\pi$, tetrad weight, or delta-derivative sign will generally fail at least one of these checks.

References

• Sasaki, Misao, and Hideyuki Tagoshi. “Analytic Black Hole Perturbation Approach to Gravitational Radiation.” Living Reviews in Relativity 6 (2003): article 6. DOI. Open PDF.
• Teukolsky, Saul A. “Perturbations of a Rotating Black Hole. I. Fundamental Equations for Gravitational, Electromagnetic, and Neutrino-Field Perturbations.” Astrophysical Journal 185 (1973): 635–647. DOI.