A 3D Summation-by-Parts scheme on a Hyperboloidal Foliation of Minkowski

Shalabh Gautam

Source-linked AI summary

A 3D Summation-by-Parts scheme on a Hyperboloidal Foliation of Minkowski

Shalabh Gautam

arXiv:2608.25363v1Aug 26, 2026gr-qcmath-phmath.APmath.NA

TL;DR

The paper addresses accurate gravitational-wave computation by presenting a full 3D Summation-by-Parts framework on hyperboloidal slices that includes future null infinity in the computational grid. The framework handles coordinate singularities, introduces energy-compatible dissipation and convergence tests, and yields two schemes with complementary accuracy and stability properties.

Problem
Accurate numerical-relativity calculations require including future null infinity in a well-posed computational formulation to capture spacetime backreaction and gravitational waves.
Method
The paper develops a full 3D Summation-by-Parts discretization in compactified hyperboloidal spherical-polar coordinates, with rescaling, coordinate-singularity treatment, energy-conserving operators, boundary-defined dissipation, and new all-grid-point convergence tests.
Results
The SBP-TEM and SBP-Stable schemes manage coordinate singularities while retaining convergence properties similar to completely regular systems; SBP-TEM favors accuracy, whereas SBP-Stable guarantees stability and negative-definite energy flux at future null infinity.
Takeaways & Limitations
The framework provides a promising basis for numerical evolution on hyperboloidal slices and possible extensions toward fully nonlinear systems such as the Einstein Field Equations.
Takeaways & Limitations
SBP-TEM can produce positive energy growth for highly noisy data near future null infinity, while formulations using I+-fixing coordinates do not yet give a completely regular system of equations.

Abstract

from arXiv · show

This paper summarises our previous work on a fully $3$D Summation-by-Parts scheme, derived for a class of linear wave equations on hyperboloidal slices on a fixed Minkowski background. The scheme is derived in spherical polar coordinates, 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, and is proved to be stable. Reducing it to the standard Cauchy problem, or to finite spacelike slices with an outer boundary, will follow a similarly. Second-order accurate finite-difference methods are used to implement this scheme numerically, but higher-order finite-difference or spectral methods could also be used. Kreiss-Oliger dissipation operators are generalized to curvilinear coordinates and are defined everywhere in the domain, including at the boundary points, such that they satisfy the dissipative property in the energy norms. We also propose new norm convergence tests that include all the grid points at all resolutions and produce more accurate results. 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.

1 Introduction

The paper targets numerical-relativity formulations that include future null infinity while avoiding outer-boundary, gauge, and extrapolation problems. It develops a full 3D SBP discretization for regularized wave systems on hyperboloidal slices, including coordinate singularities and stable boundary treatment.

Including future null infinity in the computational domain is required for highly accurate gravitational-wave calculations with spacetime backreaction.
Existing CCM, CCE, and extrapolation approaches retain well-posedness, feedback, gauge-dependence, outer-boundary, or systematic-error limitations.
Hyperboloidal foliations with conformal compactification address these limitations, although I +-fixing formulations remain analytically non-regular.
The proposed SBP scheme handles grid points at the origin and z-axis without intricate coordinate transformations or manual evolution-variable regularization.
The scheme combines energy-stable I + treatment, geometric constraint damping, and covariant artificial dissipation for hyperbolic systems.
A full 3D spatial discretization avoids spherical-harmonic decomposition and can generalize to d + 2 dimensional Minkowski spacetime.

2 Continuum Setup

The continuum setup introduces hyperboloidal coordinates that remain spacelike and reach future null infinity at finite compactified radius. Compactification and rescaling regularize the geometry and equations across the computational domain.

The outgoing and ingoing characteristic speeds are adjusted so the ingoing speed approaches zero at I +.
Hyperboloidal slices are spacelike everywhere and reach future null infinity, unlike standard Cauchy slices that approach spacelike infinity.
The height function transforms standard time as T = t + H(R), while R = R(r) maps the unbounded radial coordinate to a compactified coordinate.
The compactified radial coordinate r increases monotonically from 0 to a finite rI as R increases from 0 to infinity.
The transformed metric contains radial slicing and compactification factors through H′ and R′, while retaining the spherical angular terms.
The compactification satisfies Ω(rI) = 0, with origin conditions and asymptotic parameter n controlling the height-function behavior.

2.2 Linear Wave Equation

The paper studies a linear wave equation on Minkowski spacetime with a spatially defined potential, expressed in spherical polar coordinates. Regularized treatment is needed at the origin and for angular modes.

The model is a linear wave equation with the Minkowski d’Alembertian and a potential F defined everywhere as a function of spatial coordinates.
In spherical polar coordinates, the wave operator contains radial and angular divergence terms with explicit R and sin θ factors.
These modes can be extracted from data defined at all angular coordinates using relations provided by the formulation.

2.3 First-Order Reduction

The second-order wave equation is rewritten as a first-order reduction using Cauchy variables, constraint damping, and parity or periodicity conditions at coordinate singularities and ghost points.

The first-order reduction introduces ψT, ψR, ψθ, and ψϕ as derivatives of ψ with respect to time and spherical spatial coordinates.
The resulting system evolves the reduction variables while adding damping terms proportional to their derivative constraints.
Constraint violations decay exponentially as CR(T) = CR(0)e^-ζRT, with analogous expressions for the angular constraints.
The system is symmetric hyperbolic when the three damping parameters are equal: ζR = ζθ = ζϕ = ζ.
At the origin and z-axis, specialized reductions and parity conditions provide values needed for the singular-coordinate evolution.
Parity and periodicity rules populate ghost points beyond the θ, R, and ϕ coordinate domains during discretization.

2.4 Characteristic Variables

The formulation uses characteristic variables to obtain a system that is regular at future null infinity and to handle parity conditions at the origin.

The system is entirely regular at future null infinity when described using characteristic variables.
The characteristic formulation leads to a first-order system for the rescaled field and its radial and angular components.
At the origin, the m = 0 and l = 0 modes require parity conditions for the characteristic variables.
At the origin, the parity conditions imply ∂R(ψ+ − ψ−)(l=0) = 2∂R(ψ+)(l=0) = −2∂R(ψ−)(l=0), allowing the ψ+ and ψ− equations to be rewritten.

2.5 Introducing Hyperboloidal Coordinates: Dual Foliation approach

The dual-foliation formulation expresses characteristic equations in hyperboloidal coordinates while separating gauge choices from coordinate choices, then regularizes formally singular terms at future null infinity by rescaling.

The dual-foliation formulation uses a 3 + 1 split of the Jacobian to express equations in hyperboloidal coordinates.
The Jacobian relates (T, R, θ, ϕ) to (t, r, θ, ϕ) and decouples gauge selection from coordinate selection.
The time and radial derivatives transform as ∂T = ∂t and ∂R = (R′−1)∂t + ∂r, producing the hyperboloidal characteristic system.
Although terms involving R′ are formally singular at future null infinity, the wave-equation falloff motivates rescaling the characteristic variables.
Choosing χ monotonically with χ(0) = 1 and χ ∼ R, and taking the natural conformal rescaling, gives R/χ = r and simplifies the equations and volume element.
For n = 2, the resulting system is fully regular at future null infinity when ζ̃ = O(1) and F falls off like 1/R^2.
For F falling off like 1/R^(1+ϵ), regularity can be achieved with n ∈ (1, 1 + ϵ], whereas slower falloff makes the system singular and requires sufficient initial-data decay.

2.7 Introducing Regularized Covariant Divergence Operators

The scheme replaces a numerically infeasible radial divergence operator at infinity with regularized covariant divergence operators, while retaining endpoint equations on the origin and axis.

The spatial divergence operator is introduced componentwise along the radial, polar, and azimuthal directions.
The radial operator cannot be defined numerically at future null infinity because r → rI corresponds to R → ∞.
The regularized angular operators are ˜∂θf = ∂θf + (cot θ)f and ˜∂ϕ = ∂ϕ.
Substituting the regularized operators yields equations at the origin and, for n = 2, a corresponding system at future null infinity.
The final equations also cover the z-axis, with ζ̃ chosen as 0 or 1 for the stated purposes.

2.8 Conserved energy on hyperboloidal slices

The paper derives a conserved-energy relation for the characteristic formulation on truncated hyperboloidal slices and evaluates its limit at future null infinity.

Because the SBP scheme preserves conserved energy discretely, the formulation derives the corresponding energy on hyperboloidal slices using rescaled characteristic variables.
The slices are truncated at r0 < rI with a time-dependent outer boundary r0 = r0(t), defining the total energy on the truncated slices.
The energy change is obtained by substituting the equations of motion and applying Gauss’s theorem on the spatial slices.
The outer boundary is an incoming-null hypersurface when dr0/dt = cr.
Taking the limit r0 → rI determines the total energy change at future null infinity.
The limiting energy expression excludes angular-variable contributions ˜ψθ and ˜ψϕ because they fall off like 1/R^2.

3 Discretization

The paper constructs a semi-discrete SBP scheme for spherical-polar hyperboloidal coordinates, including coordinate-singular locations and future null infinity. The scheme enforces a discrete energy-flux condition, with alternative outer-boundary treatments trading accuracy for guaranteed stability.

3 Discretization: Spatial coordinates are discretized while time remains continuous, using grid arrays for the radial, polar, and azimuthal directions.The grid includes radial, polar, and azimuthal indices with corresponding spacings.
3 Discretization: The quadrature matrix decomposes into radial, polar, and azimuthal factors encoding grid spacing at each point.The multiplicative operators are represented as diagonal matrices, preserving continuum-like products.
3 Discretization: The SBP relations are derived by requiring the discrete energy change to equal an outer-boundary flux approaching the continuum flux.The resulting relations define the final three-dimensional SBP scheme and constrain the angular derivative operators.
3 Discretization: The SBP-TEM boundary construction preserves outer-boundary accuracy but can produce positive energy growth for states associated with a negative boundary-matrix eigenvalue.The problematic state corresponds to highly noisy data near future null infinity and can be mitigated with dissipation.
3 Discretization: Virtual ghost points populated by parity conditions define second-order finite-difference operators at coordinate endpoints and preserve derivative approximations throughout the domain.The construction handles the origin, polar endpoints, and periodic azimuthal direction while retaining antisymmetry properties in the bulk.
3 Discretization: The SBP-Stable treatment reduces outer-boundary accuracy while ensuring that the discrete energy derivative remains negative semi-definite for all times.Constraint damping leaves the SBP scheme unchanged and yields exponentially decaying discrete constraints when enabled.

4 Numerical Implementation and Results

Numerical experiments show stable, energy-conserving evolution and approximately second-order convergence for both SBP discretizations, including at coordinate singularities and null infinity. The new all-grid-point norms improve convergence assessment where conventional coarse-grid tests lose correlation.

Convergence: Energy-norm convergence approaches second order for both schemes, while pointwise convergence is second order in radial, polar, and azimuthal directions despite coordinate singularities.Radial convergence shows origin spikes, whereas polar-axis and periodic azimuthal convergence remain perfect second order.
Stability and dissipation: All simulations remain stable across resolutions up to the empirically observed CFL = 2.6785, including runs without artificial dissipation.Minimal dissipation reduces high-frequency noise, using a = 0.002 for SBP-TEM and a = 0.008 for SBP-Stable.
Energy conservation: The numerical energy agrees excellently with the analytical energy for both SBP-TEM and SBP-Stable at (Nr, Nθ, Nϕ) = (200, 50, 50).The comparison quantifies the schemes' reliability over long evolutions.
Convergence: At I +, the outgoing mode converges at second order, with higher-order wiggles diminishing as resolution increases.The convergence is integrated over the 2-sphere.
Physical behavior: The scattering-potential case produces late-time tails at I + whose power-law slopes approach a constant limit with increasing resolution and are anticipated to asymptote to −2.This behavior is associated with Price's law.
Massive Klein-Gordon equation: For F = m2, stable evolution and total discrete-energy conservation persist despite potential singularities at I +, while dissipation reduces convergence order over time only in the radial direction.The E1 and E2 norms avoid the loss of correlation seen in conventional coarse-grid convergence estimates because they use all grid points across resolutions.

5 Conclusions

The framework combines covariant regularization, energy-based SBP discretization, constraint damping, and boundary-inclusive dissipation for 3D hyperboloidal evolution. Two boundary treatments trade accuracy against stability, while new convergence tests improve assessment across resolutions.

Framework: The framework redefines the linear wave equation at coordinate singularities and uses compactified hyperboloidal coordinates with rescaling to simplify regularized divergence operators.The rescaling also reduces to conformal rescaling under conformal compactification.
Framework: Constraint damping is added geometrically to the first-order reduction, while artificial dissipation modifies its temporal part.The combined modification is expressed in geometric form.
Framework: The SBP discretization enforces energy conservation and defines spherical-polar dissipation operators everywhere, including boundary points, while satisfying energy-norm dissipativity.These properties are built into the discretization of partial derivatives and regularized covariant divergence operators.
Schemes: SBP-TEM preserves accuracy across the domain, whereas SBP-Stable guarantees stability and negative-definite energy flux at future null infinity.SBP-TEM converges better, while SBP-Stable can evolve singular massive Klein-Gordon fields.
Results: Both schemes manage curvilinear-coordinate singularities while sustaining convergence properties similar to completely regular systems.This is identified as a notable strength of the approach.
Results: Second-order finite differences capture wave propagation to future null infinity, late-time tails, and energy conservation even at low resolutions.Higher-order finite differences and pseudo-spectral discretizations can also adapt the scheme.
Convergence tests: New norm convergence tests use all grid points at all resolutions, expose higher-order numerical-error contributions, and improve convergence results where traditional tests can fail.Their effectiveness is illustrated in Figs. 12 and 13.