3d Summation-by-Parts scheme for Linear Wave Equations on Hyperboloidal Slices

Abstract: We derive a fully 3-dimensional Summation-By-Parts scheme for a class of linear wave equations on hyperboloidal slices that meet future null infinity on a Minkowski background. The scheme is derived in spherical polar coordinates, with a major strength being that it is provably stable and allows having grid points at the origin and on the $z$-axis, despite coordinate singularities, and at infinity, by introducing compactification followed by rescaling. Reducing it to the standard Cauchy problem, or on finite spacelike slices with an outer boundary, will follow a similar procedure. Interesting relations are obtained between the rescaling and compactification factors that simplify the equations, and the conditions on constraint addition terms are discovered to maintain symmetric hyperbolicity. Numerical implementation is achieved using finite-difference methods at second-order accuracy, which can be generalized to higher-order or spectral accuracies as well. Dissipation operators are given a more abstract treatment, which makes it possible to define them everywhere in the domain, including at the boundary points, in curvilinear coordinates, such that they satisfy the dissipative property (DP) in our energy norms. These generalizations reduce to the well-known Kreiss-Oliger dissipation operators whenever defined on a Cartesian grid in the bulk and satisfy the DP in the standard $L<sup>2$-norms. We also propose new norm convergence tests that produce more accurate outputs. Promising results are obtained, giving hope for application to fully nonlinear systems, like the Einstein Field Equations, and extracting the resulting gravitational waves free of systematic errors or gauge ambiguities.

• The paper presents a novel 3D SBP finite-difference scheme for linear wave equations on hyperboloidal slices that achieves energy stability and regularizes coordinate singularities.
• It employs a first-order reduction with characteristic variables and incorporates constraint damping along with tailored dissipation to ensure second-order convergence.
• The approach enables direct gravitational wave extraction at null infinity, laying a rigorous foundation for future numerical relativity and waveform modeling.

3D Summation-By-Parts Scheme for Linear Wave Equations on Hyperboloidal Slices

Introduction and Context

Accurate extraction of gravitational wave (GW) signals from numerical relativity (NR) simulations necessitates computation in domains extending to future null infinity, $\mathscr{I}^+$, where GWs become unambiguously defined and free of gauge ambiguities. Traditional approaches such as Cauchy-characteristic extraction (CCE) and Cauchy-characteristic matching (CCM) typically struggle with systematic errors, boundary reflections, or gauge dependencies. Hyperboloidal slicing — spacelike foliations that intersect $\mathscr{I}^+$ — sidesteps these issues, enabling direct evaluation of GW signals at null infinity [Zen07, Zen08]. However, formulating fully regular, provably stable numerical schemes on such foliations, especially in three spatial dimensions and in coordinates with intrinsic singularities (e.g., spherical polar), presents significant mathematical and computational challenges.

This work establishes a robust 3D finite-difference summation-by-parts (SBP) discretization for linear wave equations (LWE) posed on hyperboloidal slices with compactification at $\mathscr{I}^+$. The construction is intimately connected with subsequent generalizations to the Einstein field equations (EFEs) in generalized harmonic gauge (GHG), which share the same principal part as the LWE. The approach maintains exact energy stability, regularizes all coordinate singularities (including at the origin, the $\pm z$-axis, and at infinity), and incorporates constraint damping and dissipation in a covariant manner.

Hyperboloidal Slicing, Regularization, and Compactification

The continuum setup recasts Minkowski spacetime in spherical coordinates and introduces hyperboloidal slices through a time coordinate $t$ (the 'hyperboloidal time') and a compactified radial coordinate $r$. The transformation $T = t + H(R)$, where $H(R)$ is a 'height function', ensures that each $t = \text{const}$ slice meets $\mathscr{I}^+$ at a finite $\mathscr{I}^+$0. Radial compactification is realized by $\mathscr{I}^+$1, with $\mathscr{I}^+$2 vanishing smoothly at $\mathscr{I}^+$3.

A critical insight is the need for appropriate variable rescalings, denoted collectively as $\mathscr{I}^+$4, that both regularize the equations at $\mathscr{I}^+$5 and preserve the parity required for handling the coordinate origin and axes. The natural choice $\mathscr{I}^+$6 ensures that the energy density, divergence operators, and SBP property remain well defined and asymptotically simple.

First-Order Reduction and Symmetric Hyperbolicity

Direct discretization of the second-order LWE is suboptimal. Instead, a first-order reduction (FOR) is introduced, in both standard ($\mathscr{I}^+$7, $\mathscr{I}^+$8, $\mathscr{I}^+$9, $\mathscr{I}^+$0) and characteristic ($\mathscr{I}^+$1) variables. The constraint subsystem resulting from this reduction is damped using parameters $\mathscr{I}^+$2 in each spatial direction, ensuring exponential decay of violations in the discrete setting.

On both the $\mathscr{I}^+$3-axis and the origin — locations where spherical coordinates become degenerate — the equations are modified to be compatible with symmetry and regularity. The scheme extracts and evolves axisymmetric ($\mathscr{I}^+$4) or spherically symmetric ($\mathscr{I}^+$5) modes as appropriate, ensuring consistent evolution across the entire domain.

Hyperboloidal slicing is handled via a dual-foliation approach [Hil15, HilHarBug16], enabling the principal part of the continuum system to match the standard LWE everywhere and facilitating separation between gauge and coordinate effects.

Discrete SBP Framework in Spherical Coordinates

The primary technical contribution is a finite-difference SBP discretization in three spatial dimensions compatible with energy stability. The spatial domain is discretized on a uniform, non-staggered grid in compactified spherical-polar coordinates; the discretization directly incorporates the origin, $\mathscr{I}^+$6-axis, and outer boundary at $\mathscr{I}^+$7.

Key features include:

• Discrete covariant divergence operators: Angular and radial operators ($\mathscr{I}^+$8, $\mathscr{I}^+$9, etc.) are constructed to ensure SBP property and vanish appropriately at symmetry points. Discrete coefficients are handled with proper weightings to reflect continuum measures (e.g., factors of $\pm z$0).
• Summation-by-Parts energy norm: Both the discrete energy and the boundary flux match their continuum analogues to high order, ensuring stability under time evolution even for formally singular systems at $\pm z$1.
• Dissipation and constraint damping: Dissipation is introduced covariantly using generalized Kreiss-Oliger operators tailored to act on the relevant variables and designed to satisfy the dissipative property (DP) in the chosen energy norm. Constraint damping is likewise implemented geometrically, acting only on the spatial reduction constraints, preserving symmetric hyperbolicity.
• Boundary and singularity treatment: Discrete stencils and parity conditions are rigorously derived to enable stable evolution at the origin and on the axes, obviating the need for excision or ad hoc regularization.

Two variants of the SBP scheme are provided: SBP-TEM (which maintains uniform truncation error order throughout the domain by using wider stencils at the outer boundary) and SBP-stable (which ensures strictly negative-definite boundary flux in the discrete energy, though with locally degraded accuracy at the outer boundary).

Numerical Results and Empirical Performance

Numerical experiments spanning massless and massive wave equations as well as cases with nontrivial scattering potentials are conducted. The key empirical findings include:

• Stability and minimal artificial dissipation: Both SBP-TEM and SBP-stable schemes exhibit stability without artificial dissipation, though minor dissipation (e.g., $\pm z$2) effectively suppresses high-frequency noise.
• Convergence: Second-order convergence is observed in energy and pointwise norms across all directions, including at $\pm z$3 and at the axis and origin, consistent with the finite-difference order. Spurious oscillations in convergence estimates at low base resolution are attributed to higher-order errors and are mitigated at higher resolutions.
• Regularity at $\pm z$4 and late-time tails: The constructed scheme robustly captures physically relevant features, including propagation to and extraction at null infinity and the correct late-time tail decay exponents for solutions with scattering potentials (e.g., $\pm z$5).
• Singular source terms: For the massive LWE ($\pm z$6), despite the singularity of the mass term at infinity, the SBP-stable scheme maintains exact energy conservation at the discrete level and converges in the angular directions; only the radial direction evidences loss of convergence at late times, commensurate with the development of unresolved short-wavelength modes.

Norms incorporating all resolutions (not just those restricted to coarse grids) — as constructed via the new $\pm z$7 and $\pm z$8 definitions — are also introduced, providing more informative convergence diagnostics in challenging regimes.

Theoretical and Practical Implications

By delivering an energy-stable, fully explicit, and coordinate-agnostic 3D SBP discretization on hyperboloidal slices, this work addresses a standing challenge in numerical relativity and hyperbolic PDEs. Importantly, the framework generalizes naturally to higher-order accurate and pseudo-spectral schemes, and its structural properties (constraint preservation, symmetry, SBP, and dissipation) are directly portable to quasilinear, nonlinear, and fully coupled systems — including the 3D Einstein field equations in GHG form on hyperboloidal slices [PetGauRai23, PetGauVan24].

This scheme thus provides a rigorous mathematical foundation for future developments in:

• High-fidelity simulation of compact object mergers and isolated systems with direct GW extraction at $\pm z$9,
• Construction of waveform catalogs free from gauge/systematic extrapolation artifacts for GW data analysis and theoretical modeling,
• Stable evolution of nonlinear and matter-coupled systems on unbounded domains,
• Adaptation to generalized coordinates and curved backgrounds, including dual-foliation formulations.

Conclusion

This work constitutes a comprehensive, technically rigorous advance in the numerical treatment of hyperbolic PDEs — specifically the linear wave equation — on 3D hyperboloidal slices compactified at null infinity. The SBP scheme regularizes all coordinate singularities, preserves a discrete energy norm, and remains empirically robust and convergent under practical grid resolutions. These properties make it a compelling foundation for next-generation numerical relativity codes targeting the direct extraction of gravitational radiation at null infinity and for general hyperbolic problems in curved spacetime backgrounds.

---

References

• Anuraag Reddy, Shalabh Gautam, Prayush Kumar, "3d Summation-by-Parts scheme for Linear Wave Equations on Hyperboloidal Slices" (2606.02051)
• Additional technical background: [GauVanHil21], [Hil15], [PetGauRai23], [Zen07], [GunGarGar10], [CalLehReu03], [Str94]