Efficient and Stable Computation of Gravitational-Wave Fluxes from Generic Kerr Orbits via a Unified HeunC Framework

Abstract: Modeling extreme-mass-ratio inspirals hinges on the accurate and efficient computation of gravitational-wave fluxes from generic Kerr orbits. Conventional frequency-domain techniques are often limited by costly auxiliary parameter searches and numerical instabilities in the strong-field or high-frequency regimes. We address these challenges by reformulating both the angular and radial Teukolsky equations in terms of confluent Heun functions. Employing a hybrid analytic continuation algorithm to compute the connection coefficients eliminates the dependence on auxiliary parameters, directly yielding globally convergent solutions and scattering amplitudes. To resolve the highly oscillatory source integrands for generic orbits, we implement an adaptive bi-power mapping quadrature. Comprehensive benchmarks under standard double-precision arithmetic demonstrate that, for the total radiative flux summed over 168 low-order modes, our method achieves relative errors of order $10^{-11}$, with computational costs typically reduced by factors of 3--13 compared to the state-of-the-art GeneralizedSasakiNakamura. jl and pybhpt packages. Notably, for highly oscillatory high-order modes, our framework achieves a speedup of up to 60 times compared to specialized oscillatory integrators like GeneralizedSasakiNakamura. jl. These demonstrated gains in precision and efficiency establish the framework as a robust tool for strong-field perturbation theory, providing the numerical foundation for high-order self-force calculations and rapid, high-precision waveform generation.

• The paper introduces a unified framework utilizing HeunC functions for efficient calculation of gravitational-wave fluxes from generic Kerr orbits.
• The approach eliminates auxiliary parameter searches, achieves numerical errors as low as 10⁻¹¹, and significantly reduces computation time.
• Adaptive bi-power mapping quadrature stabilizes oscillatory integrals, yielding precise waveforms for EMRI and IMRI modeling.

Efficient and Stable Gravitational-Wave Flux Computation from Generic Kerr Orbits via HeunC Functions

Background and Motivation

Accurate and efficient computation of gravitational-wave (GW) fluxes from extreme-mass-ratio inspirals (EMRIs) and intermediate-mass-ratio inspirals (IMRIs) in Kerr spacetime is of paramount importance for modeling GW signals in the millihertz band targeted by next-generation observatories such as LISA, TianQin, and Taiji. EMRI systems provide strong-field tests of general relativity, mapping of spacetime geometry, and constraints on astrophysical environments and cosmological parameters. High-fidelity theoretical waveform templates are required for parameter estimation, necessitating sub-radian phase coherence over thousands to millions of orbits.

The established framework is black hole perturbation theory, where flux-balance relations govern orbital evolution at leading order, evaluated via the separable Teukolsky equation for gravitational perturbations. Standard frequency-domain methods, notably the Mano–Suzuki–Takasugi (MST) formalism, allow convergent expansion of homogeneous solutions, but suffer from limitations: costly searches for auxiliary parameters (e.g., renormalized angular momentum $\nu$), numerical instability in the strong-field/high-frequency regime, and bottlenecks for high eccentricity or inclination.

Recent software tools (GeneralizedSasakiNakamura.jl, pybhpt) implement sophisticated solution strategies, combining spectral expansion, semi-analytic MST series, and numerical integration, but each incurs significant computational overhead for high-order modes and generic orbits.

Unified HeunC Reformulation

The paper introduces a unified framework wherein both the angular and radial Teukolsky equations are reformulated as confluent Heun equations (CHE). The central innovation is leveraging HeunC functions as the fundamental basis for both angular and radial sectors, thus eliminating the need for an auxiliary parameter search and enabling globally convergent solutions.

Motygin's numerical algorithm constructs connection coefficients without auxiliary parameters or series expansions. Hybrid analytic continuation is employed: local Frobenius power-series expansions are matched to asymptotic forms via a bidirectional chain at an intermediate matching point. Once connection coefficients are computed and cached, evaluation of far-field solutions reduces to efficient summation.

Physical boundary conditions (ingoing/outgoing waves at horizon/infinity) are enforced directly using the HeunC basis. SWSHs and separation eigenvalues for the angular equation are determined via Wronskian vanishing, with normalization via either power-series summation or Sturm–Liouville Wronskian derivative, the latter providing numerical stability even for high $\ell$ and near-extremal spin.

For the radial sector, the transformation maps horizon and infinity to canonical CHE singularities. Scattering amplitudes (transmission, reflection, incidence) follow analytically from computed connection coefficients, bypassing branch-cut instabilities and costs associated with $\nu$-searches. The framework supplies exact global mode functions for both homogeneous and inhomogeneous cases.

Adaptive Bi-Power Mapping Quadrature

The source integrand for generic orbits exhibits highly localized oscillatory structures (particularly near apastron), rendering uniform quadrature inefficient and prone to under-sampling. The paper introduces an adaptive bi-power mapping: a coordinate transformation concentrates points where rapid phase oscillations occur, while maintaining sparse coverage in smooth regions. The mapping function and analytical Jacobian ensure numerical stability and perfect domain coverage, with linear scaling in cost.

Coupled with the HeunC solution architecture, this integration strategy suppresses quadrature error and resolves the mode amplitudes for high eccentricity/inclination and high-order harmonics.

Numerical Benchmarks

Comprehensive benchmarks are provided under standard double-precision arithmetic (machine epsilon $\sim 10^{-16}$):

• SWSH Functions and Eigenvalues: HeunC method, via Wronskian derivative, achieves relative errors $\sim 10^{-14}$, outperforming spectral and continued-fraction methods especially near extremality ($a=0.9999$). Computational time is one order of magnitude lower.
• Schwarzschild Circular Orbit Fluxes: Double-precision HeunC achieves $10^{-12}$ error with $70$ ms runtime, vs. MST and Nekrasov–Shatashvili methods requiring higher precision and $10$–$57$ s.
• Generic Kerr Orbits: For modes $\ell$0 up to $\ell$1 and orbital parameters $\ell$2 up to extremality and high eccentricity, the HeunC method maintains errors $\ell$3–$\ell$4 with runtimes 2–3 times shorter than GSN and 3–10 times shorter than pybhpt. Total radiative flux summed over $\ell$5 low-order modes achieves $\ell$6 error, with runtime varying less than $\ell$7 across spin range.
• QNM Ringdown Tests: Asymptotic amplitude errors remain below $\ell$8 for $\ell$9, $\nu$0, $\nu$1; machine-precision is achieved without instability or overhead.

Implications and Future Directions

The unified HeunC framework advances numerical stability, computational efficiency, and analytic consistency across both angular and radial sectors, establishing a robust backend for strong-field perturbation theory. Efficiency gains become pronounced for large-scale mode summation and regimes with challenging orbital parameters.

Practically, this enables rapid generation of high-precision waveform templates for EMRIs, accurate modeling of near-horizon tidal response and absorption, and robust computation of ringdown excitation factors. The framework is inherently portable: for any type-D geometry or perturbation field, governing equations are mapped to CHE parameters, permitting immediate computation of scattering and flux quantities.

This architecture streamlines solution construction—eliminating basis-matching, redundant parameter searches, and numerical instability. It is especially suited to high-order self-force calculations, adiabatic EMRI evolution, and extended source modeling. The methodology addresses longstanding bottlenecks in special function evaluation and quadrature for black hole perturbation theory.

Theoretically, the explicit parameter correspondence between HeunC and Kerr spacetime global structure may facilitate further advances in analytic studies, QNM spectrum diagnostics, and strong-field tests of general relativity.

Conclusion

The paper delivers a unified, highly efficient computational framework for GW flux calculation in Kerr spacetime, based entirely on confluent Heun functions. Motygin's hybrid analytic continuation provides robust global solutions and scattering amplitudes without auxiliary parameter overhead. Adaptive bi-power mapping quadrature resolves oscillatory source integrals, achieving high numerical precision and efficiency.

Benchmarks demonstrate relative errors $\nu$2, 2–10× speedup versus prior tools, stable performance across spin/eccentricity/inclination, and machine-precision accuracy for QNM excitation factors. The approach offers a scalable, consistent, and practically universal solution for GW flux evaluation and waveform generation in strong-field relativistic astrophysics, with broad theoretical and computational implications for the next era of high-precision GW astronomy.