Self-force calculations with numerical relativity methods

Abstract: To model gravitational waveforms from extreme mass-ratio inspirals (EMRIs) for the upcoming LISA space mission, gravitational self-force calculations are needed to second order in perturbation theory. However, to date these calculations have only been attempted for the simplest case of circular orbits in Schwarzschild spacetime. In this work, we present a new computational method aimed at performing generic second-order self-force calculations in Kerr spacetime using methods from the adjacent field of numerical relativity. We perform an $m$-mode separation of variables, add null ("$vtu$") slicing in horizon-penetrating coordinates, and solve the resulting elliptic PDEs using high-order discontinuous Galerkin discretization, adaptive mesh-refinement, and an iterative Krylov-type linear solver with parallelizable multigrid-Schwarz preconditioning. We find that our method achieves exponential convergence for the self-force on a scalar point charge in Kerr spacetime up to spins of $a=0.998$ (Thorne limit) on circular equatorial orbits as close as the ISCO (prograde and retrograde), despite the non-smooth puncture on the grid. We solve for 20 $m$-modes in parallel in a few seconds and retain the flexibility to extend the method to gravitational self-force and more generic orbits in the future. The code to perform these calculations is publicly available in the open-source numerical relativity code SpECTRE.

• The paper introduces a novel m‑mode elliptic formulation using DG hp‑AMR methods for computing self‑force in Kerr spacetimes with exponential error reduction.
• It employs advanced null slicing, horizon‑penetrating coordinates, and regularization techniques alongside GMRES solvers with multigrid preconditioning.
• The framework, implemented in the SpECTRE codebase, demonstrates rapid convergence and scalability, paving the way for high‑fidelity EMRI waveform generation.

Self-force Calculations with Numerical Relativity Methods: An Expert Overview

Motivation and Context

Accurate modeling of gravitational waveforms from extreme mass-ratio inspirals (EMRIs) is essential for forthcoming missions such as LISA. The generation of these waveforms requires precise evaluation of the gravitational self-force through second-order perturbation theory in curved spacetimes, particularly in the Kerr geometry. Historically, second-order self-force computations have only been realized for circular orbits in Schwarzschild spacetimes, primarily leveraging $l m$-mode decompositions and frequency-domain methods. However, the intrinsic lack of symmetry and the increased mode coupling in Kerr spacetime have made full separability and computational tractability substantially more challenging.

Approach: $m$-mode Elliptic Formulation with NR Techniques

The central contribution is a new approach to self-force calculations in Kerr spacetime that eschews $l m$-mode decompositions in favor of an $m$-mode Fourier basis, enabling adaptive discretization in the $(r,\theta)$ domain and simplified source construction. The key innovations are:

• Elliptic PDE formulation for each $m$-mode: The scalar-field problem sourced by a particle on a circular equatorial Kerr orbit is reduced via separation of variables to a set of 2D elliptic PDEs in $r$ and $\theta$ for each azimuthal mode $m$.
• Null slicing and horizon-penetrating coordinates: Introduction of a piecewise-constant boost ("$vtu$ slicing") enables horizon-penetrating coordinates and compactification at both the event horizon and spatial infinity, yielding improved numerical behavior and obviating the need for explicit, restrictive boundary conditions.
• Discontinuous Galerkin (DG) methods and AMR: The PDEs are discretized using high-order discontinuous Galerkin finite elements with both $m$0-refinement (splitting elements near non-smooth features) and $m$1-refinement (increasing polynomial order elsewhere), combined in an $m$2-AMR paradigm that enhances convergence even when the solution is non-smooth due to the presence of a worldline singularity.
• Advanced iterative solvers: Each discretized PDE is solved using a complex-valued GMRES Krylov method, preconditioned by a multigrid V-cycle and additive Schwarz smoothing, with subdomains solved either directly or iteratively, exploiting inherent parallelism across $m$3.
• Worldtube and regularization techniques: The domain is divided into regions inside and outside a worldtube surrounding the particle. Within the worldtube, regularization is effected through subtraction of an analytically-constructed puncture field (truncated local expansion of the Detweiler-Whiting singular field), rendering the residual field smooth enough for accurate numerical evaluation of the self-force.

Figure 1: Computational DG domain for $m$4, $m$5, $m$6; regularization is applied within the red rectangle surrounding the particle.

Numerical Implementation and Solver Infrastructure

The entire framework is implemented in the SpECTRE codebase, originally developed for large-scale numerical relativity (NR) simulations of binary black hole mergers. The code is extended to support:

• Complex-valued field arithmetic needed for $m$7-mode separation
• Dynamic insertion of analytic jump conditions at worldtube and $m$8 slicing interfaces using the DG flux infrastructure
• Efficient, lazy evaluation and caching of matrix structures for domain decomposition and AMR

The solver achieves substantial parallel efficiency, with $m$9-mode decoupling enabling trivial parallelization. The matrix structure, preconditioner design, and AMR interaction are exemplified in:

Figure 2: Matrix representation of the DG operator, showing block structure corresponding to domain decomposition; jump conditions at worldtube and $l m$0 interfaces are imposed via DG fluxes.

Figure 3: Preconditioning stack for GMRES iterations, showing multigrid hierarchy built via AMR and additive Schwarz smoothing with direct or iterative subdomain solutions.

Numerical Results and Empirical Performance

The method is validated in the context of scalar self-force calculations for a scalar charge on circular equatorial orbits in Kerr spacetimes up to extremal spin ($l m$1). Key findings include:

• Exponential convergence: Despite non-smoothness at the grid location of the worldline, the $l m$2-AMR DG scheme achieves exponential reduction of the error with increasing grid resolution, as long as refinement boundaries coincide with the singular region.

Figure 4: Exponential convergence of DG error for $l m$3-modes up to $l m$4 in the $l m$5, $l m$6 case; each mode is solved in a few seconds on a commodity processor.

• Efficient evaluation of large numbers of $l m$7-modes: Parallelized solution of 20 $l m$8-modes requires only a few seconds per mode.
• Quality of field solutions: The magnitude and structure of each $l m$9-mode field, including accurate resolution of regularized fields inside the worldtube, is preserved across spin and orbital radius configurations.

Figure 5: $m$0-mode solutions for scalar charge at ISCO for three Kerr spins; regularized field plotted inside worldtube for maximal accuracy and stability.

• Rapid solver convergence: GMRES iterations converge in fewer than 10 steps per AMR level, with full solution times dominated by one-time subdomain matrix construction.

Figure 6: GMRES convergence and walltime per AMR level for $m$1, $m$2, $m$3; demonstration of scale independence and effectiveness of multigrid preconditioning.

• Dominance of truncation error: For the total self-force composed via regularized sum over $m$4-modes, numerical resolution errors are well beneath truncation errors in the finite $m$5-mode sum.

Figure 7: Error budget in self-force extraction; numerical resolution error is subdominant to truncation error from finite $m$6-mode summation.

Implications, Applications, and Outlook

This work provides a viable, highly efficient framework for performing self-force calculations in Kerr spacetimes without recourse to $m$7-mode decompositions. The following aspects are significant:

• Practical self-force data generation: The method enables rapid, direct computation of self-forces for circular equatorial orbits in Kerr.
• Framework extensibility: The decoupling of modes, generic domain decomposition, and robust AMR open the door to extension to eccentric/inclined orbits and, crucially, the gravitational self-force with its larger, more complex PDE structure.
• Public availability: All code and implementation details are part of the open-source SpECTRE framework.

Theoretically, the use of $m$8-mode decomposition tailored with regularization permits direct evaluation of non-smooth solutions and may be pivotal for the future computation of second-order gravitational self-forces and the construction of high-fidelity EMRI waveforms for LISA and related observatories.

Future improvements could include higher-order or more sophisticated puncture fields (further accelerating $m$9-mode convergence), direct compactification to null infinity (for waveform extraction), and more advanced Helmholtz-type preconditioners or deflation strategies to further optimize solver convergence. There is also clear synergy with independent efforts in adapting coordinate systems and domain decompositions.

Conclusion

This work establishes a methodologically novel and computationally efficient paradigm for self-force computations in black hole perturbation theory by importing, adapting, and extending numerical relativity technologies. The use of $(r,\theta)$0-mode elliptic decomposition, DG $(r,\theta)$1-AMR discretizations, and multilevel iterative solvers facilitates practical, scalable, and accurate calculation of self-forces in Kerr spacetime. The approach is poised for straightforward generalization to gravitational perturbations and more general orbital configurations and positions itself as a foundational tool as the community advances toward EMRI waveform modeling for next-generation gravitational wave observatories.

(2606.04998)