FastEMRIWaveforms (FEW) Overview

• FastEMRIWaveforms (FEW) is a modular, open-source toolkit that generates and analyzes EMRI gravitational waveforms by separating intensive offline precomputations from rapid online synthesis.
• The framework features a layered architecture with dedicated modules for trajectory evolution, amplitude prediction via surrogate models, and efficient waveform summation using GPU/CPU acceleration.
• Its design leverages sparse sampling, advanced interpolation techniques, and mode selection to achieve sub-millisecond evaluations, ensuring high accuracy for EMRI data analysis.

FastEMRIWaveforms (FEW) is a modular, open-source software package for generating and analyzing gravitational waveforms from Extreme Mass Ratio Inspiral (EMRI) systems, with an initial focus on LISA-band sources (Katz et al., 2021). It was introduced as a collection of tools to build and analyze extreme mass ratio inspiral waveforms, and it was expanded from the Physical Review Letter that introduced the first fast and accurate fully-relativistic EMRI waveform template model (Katz et al., 2021). In its original formulation, the fully relativistic FEW model waveforms are for now limited to eccentric orbits in the Schwarzschild spacetime, while later work extends FEW to model eccentric equatorial inspirals into black holes with spin magnitudes $|a| \leq 0.999$ (Katz et al., 2021, Chapman-Bird et al., 11 Jun 2025). Across these implementations, the framework is characterized by separating computationally expensive offline precomputations from rapid online waveform generation via interpolation, mode selection, and hardware-accelerated summation (Chapman-Bird et al., 11 Jun 2025).

1. Origins, scope, and constituent modules

FEW is described as a highly modular design in which components can be used standalone or together, and the code structure is layered for flexibility and easy updating as waveform models improve (Katz et al., 2021). The main workflow modules are the Trajectory Module, the Amplitude Module, the Waveform Summation Module, and Utility Modules (Katz et al., 2021). The Trajectory Module computes the orbit-averaged motion using suitable ODE integration; the Amplitude Module provides for rapid prediction of all necessary harmonic mode amplitudes via reduced-order/deep-learning surrogates; the Waveform Summation Module efficiently sums all relevant harmonic modes across the EMRI trajectory; and the Utility Modules include tools for mode selection, separatrix finding, harmonic generation, frequency calculation, and mismatch analysis (Katz et al., 2021).

The package includes a Python API that allows users to easily switch between CPU and GPU computation by changing a single argument (Katz et al., 2021). For all computationally intensive routines, FEW uses C++/CUDA backends with smart device-agnostic wrappers for NumPy/CuPy in Python (Katz et al., 2021). This modular structure was later complemented by a frequency-domain implementation integrated into the FastEMRIWaveforms package for use with either graphics processing units (GPUs) or central processing units (CPUs) (Speri et al., 2023).

A plausible implication is that FEW was designed not as a single waveform approximant, but as an extensible analysis environment in which trajectory evolution, amplitude generation, and waveform synthesis can be updated independently as new self-force, Teukolsky, or surrogate ingredients become available.

2. Computational methodology and acceleration strategy

The core methodology relies on separating computationally expensive “offline” precomputations of physical data from the rapid “online” generation of specific waveforms via interpolation and integration (Chapman-Bird et al., 11 Jun 2025). In the offline phase, self-force data and spherical decomposition of mode amplitudes in the relevant parameter space are computed in advance; these involve expensive black hole perturbation theory or Teukolsky equation computations, but only need to be done once and can be parallelized (Chapman-Bird et al., 11 Jun 2025). These data are tabulated as dense multi-dimensional grids over parameters $(a, p, e)$, carefully parameterized to maximize interpolation accuracy, with specialized tapers and grids for efficiency (Chapman-Bird et al., 11 Jun 2025).

In the online phase, the system’s adiabatic trajectory is computed by interpolating the stored fluxes to drive coupled ODEs for orbital elements $(p, e, x_I)$ and phases, while mode amplitudes are rapidly interpolated (Chapman-Bird et al., 11 Jun 2025). Dense and smooth variations allow the use of spline interpolators, with tricubic for fluxes and bicubic+linear for amplitudes, leading to sub-millisecond evaluations per grid point (Chapman-Bird et al., 11 Jun 2025). Trajectory evolution, mode spectra, and waveform summation are then handled online, exploiting hardware acceleration for both interpolation and summation across many harmonics (Chapman-Bird et al., 11 Jun 2025).

Several acceleration mechanisms recur across FEW implementations. Sparse trajectory sampling exploits the slow evolution of inspirals, so that the ODE integrating the trajectory outputs only $\sim 100$ samples, which are efficiently upsampled by cubic spline during waveform synthesis (Katz et al., 2021). Efficient mode selection picks only modes contributing above a user threshold to total power, and symmetry relations imply that only $m \geq 0$ amplitudes need be directly computed (Katz et al., 2021). For the GPU implementation, spline interpolation and the large sum to construct the waveform are implemented as batched, parallelized GPU kernels, and the GPU summation is $\sim$2500x–3000x faster than CPU for waveform summation (Katz et al., 2021). In the spinning equatorial extension, GPU support on an NVIDIA A100 or similar gives more than 1000x speedup vs. CPU for large waveforms, and typical one-year FEW waveforms with thousands of harmonics can be generated in $\sim 100$ ms (Chapman-Bird et al., 11 Jun 2025).

This suggests that FEW’s speed derives less from a single surrogate component than from a pipeline architecture in which precomputation, interpolation, sparse trajectory sampling, mode truncation, and parallel reduction are all treated as first-class algorithmic objects.

3. Waveform representations and model families

The fully relativistic model for Schwarzschild eccentric orbits uses first-order gravitational self-force results (adiabatic inspiral) for motion and full Teukolsky-based amplitude computations for wave generation, restricted to Schwarzschild spacetime and thus to equatorial, eccentric orbits (Katz et al., 2021). In that case the waveform may be written as

$h = \frac{\mu}{d_L}\sum_{lmn} A_{lmn}(t)\,_{-2}Y_{lm}(\theta, \phi)e^{-i\Phi_{mn}(t)},$

with phase accumulation

$$\Phi_{\varphi,\theta,r}(t) = \int_0^t \Omega_{\varphi, \theta, r}\big(p(t'), e(t'), x_I(t')\big) dt'$$

(Katz et al., 2021).

The improved Augmented Analytic Kludge (AAK) model is a semi-relativistic kludge model designed for speed and flexibility, applicable to generic Kerr and generic inspirals (Katz et al., 2021). Its trajectory module uses robust 5PN-order inspiral ODEs for $(p, e, Y=\cos\iota)$, the frequency evolution uses actual Kerr fundamental frequencies $(a, p, e)$0, and the output is built in the LISA detector frame (Katz et al., 2021). Within FEW, the relativistic Schwarzschild model and the improved AAK model coexist as distinct waveform families with different domains of validity (Katz et al., 2021).

A separate development introduced the first ready-to-use Schwarzschild eccentric EMRI waveform implementation in the frequency domain for use with either graphics processing units (GPUs) or central processing units (CPUs) (Speri et al., 2023). In that implementation, the waveform is decomposed voice-by-voice,

$(a, p, e)$1

and the Fourier transform is generated using the stationary phase approximation (Speri et al., 2023). On GPUs, the frequency domain waveform takes in median $(a, p, e)$2 seconds to generate and is twice as fast to compute as its time domain counterpart when considering massive black hole masses $(a, p, e)$3 and initial eccentricities $(a, p, e)$4; on CPUs, the median waveform evaluation time is $(a, p, e)$5 seconds, and it is five times faster in the frequency domain than in the time domain (Speri et al., 2023).

A plausible implication is that FEW should be understood as a family of interoperable waveform-generation strategies rather than a single domain-specific approximant: time-domain and frequency-domain synthesis, fully relativistic and kludge amplitudes, and surrogate-augmented modules all appear within the same package architecture.

4. Extension to spinning eccentric equatorial inspirals

The 2025 extension addresses eccentric equatorial inspirals into black holes with spin magnitudes $(a, p, e)$6 (Chapman-Bird et al., 11 Jun 2025). The model supports eccentricities $(a, p, e)$7 and semi-latus recta $(a, p, e)$8, and includes trajectories very close to the separatrix with $(a, p, e)$9 (Chapman-Bird et al., 11 Jun 2025). It is designed for $(p, e, x_I)$0 but can be used for moderately lower $(p, e, x_I)$1, and it covers both prograde and retrograde spins, even beyond the so-called Thorne limit (Chapman-Bird et al., 11 Jun 2025). The type of orbits is eccentric equatorial inspirals, with the framework extendable to fully generic orbits once self-force data is available (Chapman-Bird et al., 11 Jun 2025).

For equatorial inspirals $(p, e, x_I)$2, the model is built around the multivoice sum

$(p, e, x_I)$3

with

$(p, e, x_I)$4

$(p, e, x_I)$5

$(p, e, x_I)$6

and

$(p, e, x_I)$7

(Chapman-Bird et al., 11 Jun 2025). The adiabatic evolution of orbital elements is written as

$(p, e, x_I)$8

$(p, e, x_I)$9

(Chapman-Bird et al., 11 Jun 2025).

Among the significant advancements listed for this version are extension to high spin and broad parameter space; fully relativistic amplitudes and harmonics incorporating thousands of $\sim 100$0 modes up to $\sim 100$1; advanced data gridding and interpolation; a Python-based, hardware-accelerated, modular codebase; robust validation; practical data analysis demonstrators; and immediate extensibility to inclined or generic orbits, merger-ringdown, post-adiabatic corrections, and environmental effects (Chapman-Bird et al., 11 Jun 2025).

5. Accuracy, mode content, and biases from reduced models

A persistent result across FEW studies is that only a modest number of harmonics is often required if those harmonics are chosen by cumulative power rather than by a fixed low-multipole truncation. In the 2021 package study, relativistic EMRI waveform templates can be generated with fewer harmonic modes $\sim 100$2–$\sim 100$3 without biasing signal extraction (Katz et al., 2021). The on-the-fly mode selector keeps all modes whose cumulative power exceeds $\sim 100$4, and only $\sim 100$510–100 harmonic modes are needed to capture $\sim 100$6 of emitted power for the purpose of data analysis (Katz et al., 2021). In the spinning equatorial extension, mode selection is based on cumulative amplitude-power up to a tolerance $\sim 100$7 (Chapman-Bird et al., 11 Jun 2025).

The same studies emphasize that reduced or kludge amplitude content is not interchangeable with full relativistic harmonic content. In the 2021 analysis, extraction of a relativistic injection with semi-relativistic amplitudes can lead to strong bias and anomalous structure in the posterior distribution for certain regions of parameter space (Katz et al., 2021). Relativistic, quadrupole-only templates have mismatch of order $\sim 100$8–$\sim 100$9, while semi-relativistic amplitudes can have mismatch up to $m \geq 0$0 at high eccentricity (Katz et al., 2021). The role of phase vs. amplitude is also explicitly differentiated: matching the phase is far more important than matching the amplitude for parameter estimation, while amplitude errors mostly manifest as mild biases unless gross (Katz et al., 2021).

In the spinning equatorial FEW model, the adiabatic error over most of parameter space is characterized by mismatches of $m \geq 0$1 with respect to error-free adiabatic models, while errors grow in regions of very high eccentricity $m \geq 0$2 or high spin near the separatrix, but remain $m \geq 0$3 for most cases (Chapman-Bird et al., 11 Jun 2025). By contrast, kludge models can overestimate the SNR by $m \geq 0$4 or underestimate by $m \geq 0$5, and induce marginal biases of up to $m \geq 0$6 in parameter recovery (Chapman-Bird et al., 11 Jun 2025). Using only the dominant $m \geq 0$7 multipoles induces similar biases, with inaccuracies scaling with spin and eccentricity (Chapman-Bird et al., 11 Jun 2025).

A common misconception is that a few harmonics are sufficient because the signal is effectively quadrupolar. The FEW results support a more specific statement: a few harmonics are often sufficient only when they are selected to cover the significant power in the signal; quadrupole-only or semi-relativistic truncations can be insufficient for precise source characterization (Katz et al., 2021).

6. Data-analysis applications, astrophysical reach, and related extensions

FEW was developed for millihertz gravitational-wave data analysis, and with the GPU-accelerated waveforms in hand, a variety of studies were performed including an analysis of EMRI mode content, template mismatch, and fully Bayesian Markov Chain Monte Carlo-based EMRI parameter estimation (Katz et al., 2021). The 2021 work states that fully Bayesian posteriors can be computed for two-year duration EMRI signals with fully relativistic templates, thanks to the speed of the FEW GPU implementation (Katz et al., 2021). The 2023 frequency-domain implementation further states that this enables, for the first time, EMRI parameter inference with fully relativistic waveforms on CPUs (Speri et al., 2023).

The spinning equatorial extension quantifies the astrophysical reach of FEW-based I/EMRI templates for LISA. For a detection threshold SNR$m \geq 0$8 and for high spin, $m \geq 0$9, the accurate waveform model gives EMRIs $\sim$0 to redshift $\sim$1 and IMRIs $\sim$2 to redshift up to $\sim$3, with only mildly dependent $\sim$4 variation on eccentricity (Chapman-Bird et al., 11 Jun 2025). For SNR $\sim$5, $\sim$6 and $\sim$7; source-frame masses are measured to $\sim$8, luminosity distance to $\sim$9, and 99% sky localization area down to $\sim 100$0 for EMRIs (Chapman-Bird et al., 11 Jun 2025). Full 12D parameter estimation confirms near-Gaussian posteriors and minimal parameter degeneracies for EMRIs, whereas IMRIs require including many harmonics because biases become significant if too few modes are kept (Chapman-Bird et al., 11 Jun 2025).

A related extension, "A Deep Learning Framework for Amplitude Generation of Generic EMRIs" (Zeng et al., 9 Mar 2026), addresses the major challenge for the fast generation of such waveforms: the generation of the Teukolsky amplitudes for generic Kerr orbits. It introduces a convolutional encoder-decoder architecture for a fast and end-to-end global fitting of the Teukolsky amplitudes, with transfer learning from Schwarzschild circular orbits to generic Kerr orbits (Zeng et al., 9 Mar 2026). Within this framework, the full harmonic amplitudes can be generated within milliseconds, while the median mode-distribution error for generic orbits is approximately $\sim 100$1, and the surrogate can be used within adiabatic waveform frameworks such as FEW (Zeng et al., 9 Mar 2026). This suggests a pathway from the eccentric equatorial Kerr extension of FEW toward generic $\sim 100$2 Kerr waveform generation once suitable amplitude datasets are available.

In aggregate, FEW occupies a specific place in EMRI and IMRI waveform modeling: it combines offline black hole perturbation theory or surrogate construction with online interpolated adiabatic evolution, mode selection, and hardware-accelerated summation, and it is used to test when reduced models remain adequate and when they induce detectable bias (Katz et al., 2021, Chapman-Bird et al., 11 Jun 2025).