Teukode: Kerr Waveform Solver

• Teukode is a time-domain solver that evolves 2+1 m-mode Teukolsky equations for gravitational perturbations of Kerr black holes.
• By implementing horizon-penetrating hyperboloidal coordinates with spatial compactification, it achieves direct waveform extraction at future null infinity.
• It robustly handles diverse source models—including non-spinning, eccentric, and spinning bodies—providing benchmarks for frequency-domain analyses.

Teukode is a time-domain solver of the Teukolsky equation for gravitational perturbations of Kerr black holes, designed to compute gravitational radiation from point-particle sources in the large–mass–ratio limit with direct waveform extraction at future null infinity. Its distinctive numerical structure is a horizon-penetrating, hyperboloidal foliation with spatial compactification, so that both the event horizon and future null infinity lie on the finite computational domain. Across its documented applications, Teukode evolves $2+1$-dimensional $m$-mode decompositions of the Teukolsky equation, reconstructs gravitational-wave strain from $\Psi_4$, computes energy and angular-momentum fluxes to infinity and, where needed, into the horizon, and serves both as a production code for inspiral–plunge–merger calculations and as a benchmark for frequency-domain perturbation theory (Harms et al., 2014).

1. Historical development and problem domain

Teukode was introduced as a framework for waveform generation from point-particle perturbations of the Kerr spacetime, with the original implementation and algorithmic foundations described by Harms, Bernuzzi, Nagar, and Zenginoglu in 2014. In that formulation, the code targeted equatorial inspiral–plunge–merger motion in the large–mass–ratio limit, where the smaller black hole is modeled as a point particle and the dynamics can be driven by an effective-one-body-resummed analytical radiation reaction. The same work presented direct extraction of waveforms at future null infinity, horizon-absorbed fluxes reconstructed in the time domain, and a self-consistent iterative procedure for computing gravitational-wave fluxes at leading order in the particle’s mass (Harms et al., 2014).

Subsequent work expanded the source modeling and orbit classes handled by the code. For eccentric equatorial motion of a spinning secondary represented in the pole–dipole Mathisson–Papapetrou–Dixon approximation, Teukode was used as the time-domain counterpart to a frequency-domain Teukolsky formalism and provided cross-checks of asymptotic fluxes and waveforms (Skoupý et al., 2021). Later, for spinning bodies on generic Kerr orbits, including off-equatorial configurations, Teukode served as the time-domain benchmark against which a new frequency-domain solver was validated at linear order in the small body’s spin (Skoupý et al., 2023). In more recent ringdown studies, it supplied the $\psi_4$-based numerical inspiral–merger–ringdown waveforms for eccentric equatorial test-mass dynamics used to calibrate closed-form merger–ringdown ansätze (Albanesi et al., 19 Mar 2026).

This trajectory places Teukode within black-hole perturbation theory as both a numerical evolution engine and a validation instrument. A plausible implication is that its role is not confined to one waveform model or one source prescription, but extends across several perturbative descriptions provided the source can be encoded consistently in the Teukolsky formalism.

2. Geometric formulation and evolution system

The defining geometric feature of Teukode is its use of horizon-penetrating hyperboloidal coordinates with radial compactification. In the 2014 formulation, the general scri-fixing transformation is written as

$$t = \bar{t} + h(\bar{r})\ ,\qquad r = \frac{\bar{r}}{\Omega(\bar{r})}\ ,$$

with a specific ${\rm HH}_S$ family

$$h(\rho) = \frac{\rho}{\Omega} - \rho - 4 M \ln \Omega\ ,\qquad \Omega(\rho) = 1 - \frac{\rho}{S}\ ,$$

so that $\rho=S$ corresponds to $\mathscr{I}^+$. In this construction, outgoing radial characteristic speeds vanish at the horizon and ingoing ones vanish at $\mathscr{I}^+$, eliminating the need for artificial boundary conditions while allowing direct extraction at null infinity (Harms et al., 2014).

In its standard $m$0-mode form, Teukode exploits axial symmetry by decomposing the field into azimuthal modes, reducing the original $m$1 problem to a set of $m$2 evolutions, one for each $m$3. For gravitational perturbations this evolution is performed for spin weight $m$4. A representative PDE form reported for the 2014 implementation is

$m$5

where the coefficient functions encode the Kerr geometry and the hyperboloidal foliation. Regularity is enforced through the rescaling

$m$6

consistent with the peeling behavior at the horizon and null infinity (Harms et al., 2014).

The generic-orbit spinning-body work describes the same underlying strategy in a slightly different notation. There, Teukode evolves the Teukolsky equation for

$m$7

in hyperboloidal coordinates and extracts fluxes at future null infinity by computing the stress-energy of the wave. The hyperboloidal, horizon-penetrating formulation is explicitly identified as the mechanism by which asymptotic behavior at the horizon and at $m$8 is handled naturally in the computational domain (Skoupý et al., 2023).

3. Source construction and dynamical inputs

Teukode’s source term is a point-particle stress-energy projected into the Newman–Penrose formalism. In the original large–mass–ratio inspiral application, the particle is non-spinning and follows equatorial motion driven by an effective-one-body Hamiltonian and resummed analytical radiation reaction. The point-particle stress–energy is inserted into the Teukolsky source through contractions such as $m$9, $\Psi_4$0, and $\Psi_4$1, together with Newman–Penrose differential operators, after transformation to the hyperboloidal evolution coordinates (Harms et al., 2014).

For spinning secondaries, Teukode incorporates a pole–dipole source consistent with the Mathisson–Papapetrou–Dixon equations under the Tulczyjew–Dixon spin supplementary condition. In the equatorial eccentric spinning-body study, the governing equations are

$\Psi_4$2

with

$\Psi_4$3

The stress–energy in pole–dipole form contains a monopole term and a derivative term involving $\Psi_4$4, and Teukode supports source representations with $\Psi_4$5-function derivatives up to third order (Skoupý et al., 2021).

The generic-orbit spinning-body validation study makes this source structure more explicit at linear order in the secondary spin. There the frequency-domain source is written in terms of tetrad projections

$\Psi_4$6

with monopole, dipole, and derivative pieces determined by $\Psi_4$7, $\Psi_4$8, $\Psi_4$9, and the Newman–Penrose spin coefficients. Teukode’s time-domain source is constructed to be consistent with the same MPD $\psi_4$0 TD SSC framework, which is why it can be used as a benchmark for frequency-domain fluxes derived from the same physical model (Skoupý et al., 2023).

This source flexibility is central to Teukode’s later uses. In the 2026 eccentric equatorial ringdown study, the source is again a non-spinning point particle, now on an eccentric equatorial trajectory parameterized by semilatus rectum $\psi_4$1, eccentricity $\psi_4$2, and relativistic anomaly $\psi_4$3, with

$\psi_4$4

There the particle worldline is supplied by an EOB dynamics including non-circular corrections, and Teukode is used to generate the $\psi_4$5 waveform data against which ringdown models are fitted (Albanesi et al., 19 Mar 2026).

4. Waveforms, strain reconstruction, and flux observables

The principal radiative variable in Teukode is the Newman–Penrose scalar $\psi_4$6. The asymptotic relation used in the spinning-body studies is

$\psi_4$7

with $\psi_4$8, while the 2014 and 2026 descriptions use the equivalent statement

$\psi_4$9

Operationally, Teukode evolves curvature perturbations, extracts the field at future null infinity, and obtains strain modes by time integration (Skoupý et al., 2021).

For fluxes at infinity, the 2014 formulation gives

$$t = \bar{t} + h(\bar{r})\ ,\qquad r = \frac{\bar{r}}{\Omega(\bar{r})}\ ,$$0

$$t = \bar{t} + h(\bar{r})\ ,\qquad r = \frac{\bar{r}}{\Omega(\bar{r})}\ ,$$1

The same work also gives time-domain formulae for horizon absorption in terms of the $$t = \bar{t} + h(\bar{r})\ ,\qquad r = \frac{\bar{r}}{\Omega(\bar{r})}\ ,$$2 field $$t = \bar{t} + h(\bar{r})\ ,\qquad r = \frac{\bar{r}}{\Omega(\bar{r})}\ ,$$3 evaluated at $$t = \bar{t} + h(\bar{r})\ ,\qquad r = \frac{\bar{r}}{\Omega(\bar{r})}\ ,$$4, following Poisson’s formalism, and reconstructs absorbed mass and angular-momentum fluxes through the auxiliary functions $$t = \bar{t} + h(\bar{r})\ ,\qquad r = \frac{\bar{r}}{\Omega(\bar{r})}\ ,$$5 (Harms et al., 2014).

In the frequency-domain studies for spinning bodies, Teukode is used to benchmark mode-resolved asymptotic fluxes derived from amplitudes $$t = \bar{t} + h(\bar{r})\ ,\qquad r = \frac{\bar{r}}{\Omega(\bar{r})}\ ,$$6. For generic Kerr orbits, the averaged fluxes are written as

$$t = \bar{t} + h(\bar{r})\ ,\qquad r = \frac{\bar{r}}{\Omega(\bar{r})}\ ,$$7

with per-mode contributions

$$t = \bar{t} + h(\bar{r})\ ,\qquad r = \frac{\bar{r}}{\Omega(\bar{r})}\ ,$$8

$$t = \bar{t} + h(\bar{r})\ ,\qquad r = \frac{\bar{r}}{\Omega(\bar{r})}\ ,$$9

Teukode does not compute these amplitudes in the frequency-domain form; instead, it provides direct time-domain fluxes at scri, which are then averaged and, when necessary, linearized numerically in the dimensionless secondary spin parameter ${\rm HH}_S$0 for comparison with the ${\rm HH}_S$1 theory (Skoupý et al., 2023).

A recurring practical point is that Teukode’s direct extraction at ${\rm HH}_S$2 removes the need for finite-radius extrapolation in production use. The 2014 paper nevertheless discusses a polynomial ${\rm HH}_S$3 extrapolation formula as a comparison tool when scri extraction is not used, but explicitly identifies direct scri extraction as preferable (Harms et al., 2014).

5. Numerical implementation and accuracy properties

Teukode is implemented as a method-of-lines evolution scheme with high-order finite differences in space and explicit Runge–Kutta time stepping. The 2014 paper reports that the ${\rm HH}_S$4 Teukolsky equation is cast first order in time and second order in space for reduction variables ${\rm HH}_S$5, with spatial derivatives tested up to eighth order and production runs typically using sixth order, while time integration uses fourth-order Runge–Kutta (Harms et al., 2014).

The spinning-generic-orbit validation paper states the time-domain scheme in the same high-order configuration most relevant for benchmarking: sixth-order finite differences in space and fourth-order Runge–Kutta in time, with ${\rm HH}_S$6-mode resolved fluxes output at future null infinity. That work also emphasizes that Teukode is robust for off-equatorial, generic Kerr motion and that trajectories driving the source are obtained from full MPD integrations, with the difference between those trajectories and linearized ones being ${\rm HH}_S$7 (Skoupý et al., 2023).

A numerically delicate aspect is the particle ${\rm HH}_S$8-function representation. In the 2014 implementation, two classes were tested: narrow Gaussian approximations and discrete ${\rm HH}_S$9-point stencils enforcing integral properties. For moving sources during inspiral–plunge, Gaussian smoothing in the radial direction together with a discrete representation in $$h(\rho) = \frac{\rho}{\Omega} - \rho - 4 M \ln \Omega\ ,\qquad \Omega(\rho) = 1 - \frac{\rho}{S}\ ,$$0 was found to yield stability and accuracy (Harms et al., 2014). The 2021 spinning-body paper extends this discussion to sources with $$h(\rho) = \frac{\rho}{\Omega} - \rho - 4 M \ln \Omega\ ,\qquad \Omega(\rho) = 1 - \frac{\rho}{S}\ ,$$1-function derivatives up to third order. It reports that piecewise polynomials are typically more accurate and faster for circular equatorial orbits, while Gaussian approximations are more stable when the particle moves in $$h(\rho) = \frac{\rho}{\Omega} - \rho - 4 M \ln \Omega\ ,\qquad \Omega(\rho) = 1 - \frac{\rho}{S}\ ,$$2 or $$h(\rho) = \frac{\rho}{\Omega} - \rho - 4 M \ln \Omega\ ,\qquad \Omega(\rho) = 1 - \frac{\rho}{S}\ ,$$3, especially for highly eccentric motion where rapid changes of the source relative to the grid increase numerical noise (Skoupý et al., 2021).

The code’s documented performance reflects its geometric setup. The 2014 work reports that $$h(\rho) = \frac{\rho}{\Omega} - \rho - 4 M \ln \Omega\ ,\qquad \Omega(\rho) = 1 - \frac{\rho}{S}\ ,$$4 is approximately twice as fast as $$h(\rho) = \frac{\rho}{\Omega} - \rho - 4 M \ln \Omega\ ,\qquad \Omega(\rho) = 1 - \frac{\rho}{S}\ ,$$5 while avoiding artificial dissipation, and that strong-field inspirals at $$h(\rho) = \frac{\rho}{\Omega} - \rho - 4 M \ln \Omega\ ,\qquad \Omega(\rho) = 1 - \frac{\rho}{S}\ ,$$6 can reach $$h(\rho) = \frac{\rho}{\Omega} - \rho - 4 M \ln \Omega\ ,\qquad \Omega(\rho) = 1 - \frac{\rho}{S}\ ,$$7 in approximately two weeks on a desktop in serial GNU C (Harms et al., 2014). The 2026 eccentric ringdown study reports typical resolutions $$h(\rho) = \frac{\rho}{\Omega} - \rho - 4 M \ln \Omega\ ,\qquad \Omega(\rho) = 1 - \frac{\rho}{S}\ ,$$8, with $$h(\rho) = \frac{\rho}{\Omega} - \rho - 4 M \ln \Omega\ ,\qquad \Omega(\rho) = 1 - \frac{\rho}{S}\ ,$$9 runs used for convergence checks (Albanesi et al., 19 Mar 2026).

6. Validation, scientific applications, and present scope

Teukode has been validated against both frequency-domain black-hole perturbation calculations and independent time-domain solvers. In the original 2014 study, circular-orbit fluxes at $\rho=S$0 and into the horizon matched frequency-domain results to $\rho=S$1 at $\rho=S$2, radial plunge waveforms in Schwarzschild agreed with a $\rho=S$3 Regge–Wheeler–Zerilli code in amplitude, phase, and tail decay, and low-resolution self-convergence showed the expected fourth-order rates (Harms et al., 2014).

For eccentric equatorial spinning bodies, the 2021 study used Teukode to cross-check a newly developed frequency-domain formalism for asymptotic fluxes. Because the frequency-domain calculation uses spin-weighted spheroidal harmonics whereas the time-domain decomposition in Teukode uses spin-weighted spherical harmonics, single $\rho=S$4 fluxes retain a projection mismatch. The paper therefore stresses that agreement should be assessed after summing over $\rho=S$5 at fixed $\rho=S$6, in which case the relative differences converge toward zero with increasing resolution. It also documents that the choice of $\rho=S$7-representation becomes increasingly important as eccentricity rises (Skoupý et al., 2021).

In the 2023 generic-orbit spinning-body work, Teukode plays a more specialized validation role. The authors drive the code with trajectories obtained from the full nonlinear-in-spin MPD equations, then numerically differentiate the resulting fluxes with respect to $\rho=S$8 to isolate the linear-in-spin contribution. For nearly spherical orbits with $\rho=S$9, $\mathscr{I}^+$0, and $\mathscr{I}^+$1, the relative differences in the linear-in-spin energy flux per $\mathscr{I}^+$2-mode are reported as typically at or below a few $\mathscr{I}^+$3, with a maximum observed discrepancy $\mathscr{I}^+$4. For generic eccentric and inclined cases such as $\mathscr{I}^+$5, the discrepancies are typically between $\mathscr{I}^+$6 and $\mathscr{I}^+$7, and remain around $\mathscr{I}^+$8 even for higher eccentricity examples (Skoupý et al., 2023).

That same study assigns Teukode a broader interpretive significance. By comparing time-domain fluxes for generic spin orientations against frequency-domain fluxes computed with the small-body spin aligned to the orbital angular momentum, it shows that at linear order in spin only the $\mathscr{I}^+$9 harmonics contribute to the fluxes, while the $\mathscr{I}^+$0 harmonics associated with the perpendicular spin component contribute only at quadratic order. This supports the claim that aligned-spin linear-in-spin fluxes capture the complete $\mathscr{I}^+$1 asymptotic energy and angular-momentum fluxes needed for adiabatic EMRI modeling even when the spin is not aligned (Skoupý et al., 2023).

In the 2026 eccentric equatorial ringdown study, Teukode is used not primarily as a flux benchmark but as the waveform generator underlying phenomenological ringdown modeling. There the post-merger model is anchored at a time closely related to the light-ring crossing, identified operationally by the peak of the orbital frequency $\mathscr{I}^+$2, with $\mathscr{I}^+$3. The Teukode waveforms support modeling of all multipoles with $\mathscr{I}^+$4 up to $\mathscr{I}^+$5, together with $\mathscr{I}^+$6, $\mathscr{I}^+$7, $\mathscr{I}^+$8, and $\mathscr{I}^+$9, including spherical-spheroidal mode mixing and beating between co-rotating and counter-rotating quasi-normal modes (Albanesi et al., 19 Mar 2026).

The present scope of Teukode, as documented in these works, remains that of a perturbative large–mass–ratio code. The original implementation assumes equatorial motion, later studies extend source modeling to eccentric equatorial spinning bodies and to generic Kerr motion for validation purposes, and ringdown applications focus on eccentric equatorial test-mass trajectories. This suggests that Teukode’s main scientific niche is high-accuracy Kerr perturbation theory with direct asymptotic extraction, especially where one needs either time-domain waveform generation beyond strictly periodic motion or an independent benchmark for frequency-domain constructions.