Estimating High-Order Time Derivatives of Kerr Orbital Functionals

TL;DR

The paper tackles the challenge of extracting high-order coordinate-time derivatives of Kerr orbital functionals by comparing two Fourier-based reconstruction approaches: (i) analytic extraction of coordinate-time coefficients via mapping from Mino-time expansions and (ii) fitting a time-series to a truncated coordinate-time Fourier series. It introduces a hybrid method that fits in Mino time, differentiates with respect to Mino time, and then converts to coordinate time, achieving superior accuracy for high-order derivatives with relatively few harmonics. A robust recursive procedure is also presented to compute high-order derivatives directly from the geodesic equations and translated to coordinate time, providing exact checks for functionals with explicit Boyer–Lindquist representations. Applications to a simple test function and to mass/current multipoles in harmonic coordinates demonstrate that the hybrid Mino-time fit yields ~10^{-6} fractional residuals for sixth derivatives and offers substantial computational efficiency, advancing accurate gravitational-wave predictions for EMRIs and kludge models. The work provides a general framework for high-order derivative evaluation along Kerr worldlines and suggests avenues for future improvements, including joint fits and automatic differentiation to further enhance robustness and speed.

Abstract

Functions of bound Kerr geodesic motion play a central role in many calculations in relativistic astrophysics, ranging from gravitational-wave generation to self-force and radiation-reaction modeling. Although these functions can be expressed as a Fourier series using the geodesic fundamental frequencies, reconstructing them in coordinate time is challenging due to the coupling of the radial and polar motions. In this paper, we compare two strategies for performing such reconstructions and their ability to estimate high-order coordinate-time derivatives of the orbital functional. The first method maps Fourier coefficients from Mino to coordinate time; the second method fits a sampled time series of the function to a truncated coordinate-time Fourier series. While the latter method is prone to overfitting, it yields more accurate reconstructions and derivatives than the mapping, but completely misrepresents the harmonic content of the orbital functional. For the purpose of accurate coordinate-time derivative estimation, we propose a hybrid method: fit for the Mino-time coefficients, differentiate with respect to Mino time, then convert to coordinate time. Applied to the mass quadrupole of a generic Kerr geodesic, this hybrid method recovers the sixth derivative with a fractional residual $\sim10^{-6}$ using only two harmonics. For orbital functionals that depend explicitly on the geodesic orbit expressed in Boyer--Lindquist coordinates, we also provide a recursive procedure for computing coordinate-time derivatives using exact analytic expressions. These results offer a general framework for accurately evaluating high-order time derivatives along Kerr geodesic worldlines, with direct relevance to applications such as extreme-mass-ratio inspiral kludge waveform modeling, where such derivatives are key ingredients for precise gravitational-wave predictions.

Figures (3)

• Figure 1: Reconstructed power spectrum of the orbital functional $f(t)=r\cos{\theta}$(left) and fractional residual $\epsilon_{\mathrm{frac}}=|1-\tilde{f}(t)/f(t)|$(right), where $\tilde{f}(t)$ is the time series of $f$ reconstructed using the analytic (dashed) or Fourier fit (solid) methods. For clarity, we show only one dashed curve for $N_\mathrm{harm}=20$ in the left panel since the power spectra obtained from the analytic formula for the Fourier coefficients does not change pointwise as $N_\mathrm{harm}$ is increased, so all the curves for lower $N_\mathrm{harm}$ overlap. The underlying geodesic was sampled at a rate greater than twice the maximum frequency shown in the left panel, i.e., the Nyquist limit is greater than the upper $x$ limit. The fractional residual in the right panel spikes at some points as the orbital functional crosses zero and changes sign. As seen in the left panel, the power spectrum obtained using the fitting method changes with the number of terms included in the truncated expansion but does not converge to that obtained from the analytic method with increasing $N_\mathrm{harm}$. Nevertheless, the right panel shows that the fitting method provides time series reconstructions which are orders of magnitude more accurate than the analytic method for a given truncated expansion.
• Figure 2: Derivative estimation fractional residual $\epsilon_{\mathrm{frac}}=|1-\tilde{f}^{(N)}(t)/f^{(N)}(t)|$, where $\tilde{f}^{(N)}(t)$ is the estimated $N$th (second, fourth, and sixth) time derivative obtained from the reconstructed time series of $f(t)=r\cos{\theta}$ from the analytic (dashed) or Fourier fit (solid) methods. The underlying Fourier coefficients used to evaluate Eq. \ref{['eq:Kerr:BLExpansionDeriv']} are the same as those depicted in Fig. \ref{['fig:SimpleTestFunction:PowerResidualSubplot']} for $N_\mathrm{harm}=20$. Although the fitting method fails to correctly capture the harmonic structure of $f$ due to overfitting, the resulting increased accuracy in the reconstructed $\tilde{f}$ relative to the analytic method yields derivative estimations which are more accurate than the analytic method by several orders of magnitude for each derivative order.
• Figure 3: Reconstruction and derivative estimation fractional residuals $\epsilon_{\mathrm{frac}}=|1-\tilde{f}^{(N)}(t)/f^{(N)}(t)|$, where $\tilde{f}^{(N)}(t)$ is the estimated $N$th (zeroth, second, fourth, and sixth) time derivative of $f(t)=\ddot{M}_{12}$, a mass quadrupole component. These residuals are obtained by fitting a truncated Mino time expansion of $f$ for its $\lambda$-coefficients, evaluating the Mino time version of Eq. \ref{['eq:Kerr:BLExpansionDeriv']}, and then converting these to derivatives with respect to coordinate time. Left: a fit with $N_\mathrm{harm}=2$ and $n_{p}=101$ which has $\epsilon_{\mathrm{frac}}\sim 10^{-6}$ near the center of the array for the estimated $\tilde{f}^{(6)}$ (i.e., the eighth coordinate time derivative of $M_{12}$). Right: a fit with $N_\mathrm{harm}=5$ and $n_{p}=501$ which has $\epsilon_{\mathrm{frac}}\sim10^{-7}$ near the center for the estimated $\tilde{f}^{(6)}$. In this example, the underlying geodesic was evolved for a time $\Delta\lambda=0.5\Lambda_{\mathrm{min}}$, where $\Lambda_{\mathrm{min}}=2\pi/\max\{\Upsilon_{r}, \Upsilon_{\theta}, \Upsilon_{\phi}\}$. Though the fit shown in the right panel has many more terms in the truncated Fourier expansion and a higher sampling rate of $f$, the fit in the left panel is less accurate by only roughly an order of magnitude, while being computationally faster by a factor about two orders of magnitude with respect to the case shown in the right panel.