Spectral Methods in Numerical Relativity

Spectral methods approximate a function by projecting it onto a set of global basis functions. Unlike finite-difference or finite-element methods, whose accuracy scales as a fixed power of the grid spacing, spectral methods achieve exponential convergence for smooth problems — a property that is indispensable in numerical relativity, where the dynamical fields must be resolved with extreme precision over long evolution times.

1.1 — Fourier / Modal Expansion

Given a smooth function on a computational domain, we expand it in a truncated series of orthogonal basis functions:

Here the φ_k are the basis functions (Fourier modes, Chebyshev polynomials, or spherical harmonics depending on the geometry), and the û_k are the spectral coefficients determined by the inner product of u with each basis function. The coefficients carry the full information content of the discretisation: once they are known, the solution is available everywhere in the domain without interpolation.

1.2 — Exponential vs Algebraic Convergence

The hallmark of spectral methods is the rate at which the truncation error decreases as the number of modes N grows. For a C^∞ (infinitely differentiable) function, the spectral coefficients decay faster than any polynomial:

Compare this with a standard finite-difference scheme of order p, whose error behaves as:

The exponential convergence means that doubling the number of spectral modes can gain several orders of magnitude in accuracy, whereas doubling the number of grid points in a finite-difference scheme only improves accuracy by a constant factor 2^p. This distinction is critical in numerical relativity: extracting gravitational waveforms at distances of hundreds of M from the source requires errors below 10⁻¹⁰, a regime where algebraic convergence demands prohibitively fine grids.

1.3 — Lab: Convergence Comparison

The code below compares the convergence rate of a Chebyshev spectral approximation against finite-difference schemes of orders 2 and 4 for the function f(x) = sin(πx) on [−1, 1].

For problems on a finite, non-periodic interval — which is the generic case in numerical relativity — the natural basis is the family of Chebyshev polynomials. They inherit the optimal approximation properties of Fourier series through the simple change of variable x = cos(θ), which maps the interval [−1, 1] to the periodic domain [0, π].

2.1 — Chebyshev Polynomials

The Chebyshev polynomial of degree n is defined by:

These polynomials satisfy the orthogonality relation with respect to the weight function w(x) = (1 − x²)^−1/2 and obey a three-term recurrence: T₀ = 1, T₁ = x, T_n+1 = 2x T_n − T_n−1. Their extremal property — among all monic polynomials of degree n, T_n/2ⁿ⁻¹deviates least from zero in the maximum norm — makes them the ideal interpolation basis.

2.2 — Gauss-Lobatto Collocation Points

In practice, spectral methods in NR almost always use the pseudospectral (collocation) approach: the unknown is represented by its values at specific nodes rather than by its modal coefficients. The optimal choice of nodes for Chebyshev methods is the Gauss-Lobatto grid:

These points cluster near the boundaries of the interval, which counteracts the Runge phenomenon and guarantees stable interpolation. Crucially, the endpoints x = ±1 are always included, which simplifies the imposition of boundary conditions.

2.3 — The Differentiation Matrix

On the collocation grid, differentiation reduces to multiplication by a dense matrix D. If u is the vector of nodal values [u(x₀), ..., u(x_N)]^T, then the vector of derivative values is u′ = D u. The entries of D are obtained from the analytical derivative of the Lagrange interpolant through the Gauss-Lobatto nodes. Second derivatives are computed as D² = D · D, giving us a direct way to discretise elliptic and parabolic operators.

2.4 — Lab: Spectral BVP Solver

We construct the Chebyshev differentiation matrix and use it to solve the boundary value problem u″ = −π² sin(πx) with u(−1) = u(1) = 0, whose exact solution is u(x) = sin(πx).

A single spectral domain cannot efficiently resolve a binary black hole spacetime: the geometry contains two excision regions (one per black hole), a near zone with steep gradients, and a wave zone extending to hundreds of M. The solution is to tile the computational domain with multiple subdomains, each carrying its own spectral expansion, and couple them through interface conditions.

3.1 — Multi-Domain Strategy

In the SpEC code (precursor to SpECTRE), a typical BBH domain contains O(60–80) subdomains: spherical shells around each black hole, cylindrical regions along the orbital axis, and a series of concentric spherical shells in the wave zone. At each interface between adjacent subdomains, the solution and its normal derivative must be continuous.

3.2 — Penalty Method for Interface Coupling

Rather than enforcing exact matching (which destroys the banded structure of the time-stepper), modern spectral NR codes use a penalty method. At each interface, a term proportional to the mismatch between the two sides is added to the evolution equation:

The penalty parameter τ must be chosen large enough to suppress interface errors without introducing numerical instability. For hyperbolic systems, the optimal choice is related to the characteristic speeds at the interface. This approach is the spectral analogue of the discontinuous Galerkin (DG) method and is sometimes called the spectral penalty method or multi-domain DG-spectral method.

3.3 — Interface Matching Conditions

For the first-order reduction of Einstein’s equations commonly used in spectral NR codes, the interface conditions are imposed on the characteristic variables: incoming modes take their values from the neighbouring domain, while outgoing modes retain their local values. This upwind-type flux ensures that information propagates correctly across domain boundaries and is essential for stability of the coupled multi-domain system.

3.4 — Lab: Multi-Domain Wave Equation

We solve the 1D wave equation u_tt = u_xx on two abutting spectral domains [−1, 0] and [0, 1], coupled via the penalty method. A Gaussian pulse propagates seamlessly across the interface.

SpECTRE is the next-generation spectral/DG code for numerical relativity, developed primarily at Caltech, Cornell, and CITA. It supersedes the earlier SpEC code and is designed from the ground up for modern supercomputers with hundreds of thousands of cores. SpECTRE combines a discontinuous Galerkin (DG) scheme with spectral accuracy, task-based parallelism, and adaptive mesh refinement.

4.1 — Task-Based Parallelism with Charm++

Traditional MPI-based codes assign each subdomain to a fixed MPI rank. Load imbalances (e.g., when one black hole passes through a processor boundary) lead to idle cores and wasted wall-clock time. SpECTRE instead uses the Charm++ runtime system, which decomposes the computation into lightweight chares (tasks) that can be migrated between processors at runtime. The runtime scheduler automatically overlaps computation with communication, hiding latency and balancing load dynamically.

4.2 — Domain Geometry

SpECTRE supports a rich library of coordinate maps that compose to build complex domain geometries:

• Spherical shells: concentric shells around each black hole and in the wave zone. The inner boundary sits just inside the apparent horizon (excision surface).
• Cubes / Wedges: fill the regions between spherical shells and provide smooth transitions between topologically different domains.
• Cylinders: wrap the orbital axis and accommodate the elongated geometry of the binary’s near zone.
• Time-dependent maps: the domain deforms in time to track the orbiting, inspiralling black holes, keeping the excision surfaces co-moving with the horizons.

4.3 — Code Structure

SpECTRE is written in modern C++ (C++17/20) with heavy use of template metaprogramming for compile-time configuration. The key layers are:

Each element (subdomain) is a Charm++ chare that independently advances its local state, communicates fluxes with neighbours, and periodically writes observational data. The entire simulation is expressed as a directed acyclic graph (DAG) of tasks, and the Charm++ runtime schedules these tasks across available processors.

4.4 — Lab: Task Graph Visualization

We demonstrate the task-based parallelism concept by building and visualizing a simple dependency graph that models the operations on a single SpECTRE element during one time step.

Even with spectral accuracy, a fixed resolution may be wasteful in smooth regions and insufficient near sharp features (e.g., matter surfaces in neutron-star mergers, or gauge dynamics near punctures). hp-Adaptive mesh refinement (hp-AMR)combines two strategies to allocate resolution optimally.

5.1 — Error Estimation

On each element K, an a-posteriori error indicator is computed from the spectral coefficients of the local solution. A common choice is the truncation error estimator based on the highest-mode coefficients:

If η_K exceeds a tolerance, the element is flagged for refinement. If it is well below the tolerance, the element may be coarsened to save computational resources.

5.2 — h-Refinement vs p-Refinement

Two complementary refinement strategies are available:

• h-refinement: split the element into smaller sub-elements, keeping the polynomial degree fixed. This is effective when the solution has localised features — discontinuities, sharp gradients, or singularities — where increasing polynomial order would not help.
• p-refinement: increase the polynomial degree within the element, keeping the element boundaries fixed. This is optimal for smooth solutions, where the spectral coefficients decay rapidly and adding a few more modes yields exponential accuracy gains.

5.3 — The hp Decision Criterion

The key question is: when should we split (h) and when should we enrich (p)? The standard approach examines the decay rate of the spectral coefficients:

• If the coefficients decay exponentially (smooth solution): choose p-refinement.
• If the coefficients decay algebraically or plateau (non-smooth solution): choose h-refinement.

A practical diagnostic is to fit log|û_k| to both a linear model (exponential decay) and a power-law model, then pick the strategy whose model fits better.

5.4 — Lab: Automatic hp Decision

We approximate a function that is smooth in one region and has a sharp feature (a narrow Gaussian spike) in another. The code automatically decides whether to use h- or p-refinement in each element based on the spectral coefficient decay rate.

The Laser Interferometer Space Antenna (LISA), scheduled for launch in the mid-2030s, will observe gravitational waves in the millihertz band. Its primary sources — massive black hole binaries with total mass 10⁵–10⁷ M_⊙ — spend up to 10⁵ orbits in the LISA sensitivity band before merger. Producing accurate numerical waveform templates for these signals is one of the grand challenges motivating spectral methods in NR.

6.1 — The Long-Inspiral Challenge

Ground-based detectors (LIGO, Virgo) observe the final tens of orbits before merger, where finite-difference codes like those built on Cactus/Carpet can produce waveforms of sufficient accuracy in a few days of wall-clock time. LISA, by contrast, demands waveforms covering 10⁴–10⁵ orbits. The cumulative phase error over such a long evolution must remain below ∼0.1 radians for matched filtering to be effective.

Finite-difference methods with algebraic convergence accumulate truncation error proportional to T · h^p over evolution time T. For 10⁵ orbits this becomes unmanageably large unless the grid is refined to extreme levels. Spectral methods, with their exponential convergence, keep the per-step error at O(e^−σN), and the cumulative error remains controllable even over very long evolutions.

6.2 — Cactus (FD) vs SpECTRE (Spectral)

6.3 — Waveform Extraction

Gravitational waveforms are extracted at large radius by decomposing the Weyl scalar Ψ₄ into spin-weighted spherical harmonics₋₂Y_lm(θ, φ). The spectral representation is natural here: the angular decomposition is itself a spectral expansion, and the radial extraction at multiple radii enables Richardson extrapolation to spatial infinity. SpECTRE’s Cauchy-characteristic extraction (CCE) module propagates the waveform to future null infinity ℐ⁺, eliminating finite-radius gauge artifacts entirely.

6.4 — Lab: Quantum Encoding of Spectral Coefficients

As an exploratory exercise, we encode a vector of Chebyshev spectral coefficients as quantum amplitudes in a multi-qubit state. This demonstrates the basic idea behind quantum algorithms that could, in principle, manipulate spectral representations of gravitational-wave data in superposition.

Further Reading & References