Systematic errors in fast relativistic waveforms for Extreme Mass Ratio Inspirals

Accurate modeling of \gls{EMRIs} is essential for extracting reliable information from future space-based gravitational wave observatories. Fast waveform generation frameworks adopt an offline/online architecture, where expensive relativistic computations (e.g. self-force and black hole perturbation theory) are performed offline, and waveforms are generated rapidly online via interpolation across a multidimensional parameter space. In this work, we investigate potential sources of error that result in systematic bias in these relativistic waveform models, focusing on radiation-reaction fluxes. Two key sources of systematics are identified: (i) the intrinsic inaccuracy of the flux data, for which we focus on the truncation of the multipolar mode sum, and (ii) interpolation errors from transitioning to the online stage. We quantify the impact of mode-sum truncation and analyze interpolation errors by using various grid structures and interpolation schemes. For circular orbits in Kerr spacetime with spins larger than $a \geq 0.9$, we find that $\ell_{\text{max}} \geq 30$ is required for the necessary accuracy. We also develop an efficient Chebyshev interpolation scheme, achieving the desired accuracy level with significantly fewer grid points compared to spline-based methods. For circular orbits in Kerr spacetimes, we demonstrate via Bayesian studies that interpolating the flux to a maximum global relative error that is equal to the small mass ratio is sufficient for parameter estimation purposes. For 4-year long quasi-circular EMRI signals with SNRs $= \mathcal{O}(100)$ and mass-ratios $10^{-4}-10^{-6}$, a global relative error of $10^{-6}$ yields mismatches $<10^{-3}$ and negligible parameter estimation biases.