Cactus / Perspectives of Quantum Computing in Numerical GR

The Arnowitt-Deser-Misner (ADM) decomposition splits the four-dimensional spacetime metric into spatial slices evolved forward in a coordinate time t. The lapse function α controls the rate of proper time advance between slices, while the shift vector βⁱgoverns how spatial coordinates propagate from one slice to the next.

This 3+1 split is the foundation of all modern numerical relativity codes. By foliating spacetime into a family of spacelike hypersurfaces Σ_t, we convert Einstein’s covariant field equations into a Cauchy initial value problem: specify constraint-satisfying initial data on Σ₀, then evolve forward using hyperbolic PDEs. The gauge freedom of general relativity is encoded in the free choice of lapse and shift, which determine how the coordinate grid threads through spacetime.

Physically, the lapse α measures how much proper time elapses between adjacent slices per unit coordinate time. Setting α < 1 near a black hole (“maximal slicing” or “1+log” gauge) prevents the slice from hitting the singularity — the famous “singularity avoidance” property. The shift vector βⁱ prevents coordinate stretching by co-moving the grid with the black holes, which is essential for long-term stability of binary evolutions.

The ADM Line Element

Here γ_ij is the induced 3-metric on each spatial hypersurface Σ_t. The extrinsic curvature K_ij encodes how each slice is embedded in the four-dimensional manifold and satisfies the evolution equation relating it to α, βⁱ, and the spatial Ricci tensor.

The ADM evolution equations consist of 12 first-order-in-time PDEs for the 6 independent components of γ_ij and 6 independent components of K_ij. However, Einstein’s equations also impose 4 constraint equations (1 Hamiltonian + 3 momentum) that must be satisfied on every time slice. In the continuum these are preserved by the evolution; numerically, truncation errors cause them to drift, and the raw ADM system amplifies these violations exponentially — motivating the BSSN reformulation discussed in Chapter 2.

Video Lectures

Interactive Simulation

The Baumgarte-Shapiro-Shibata-Nakamura (BSSN) system reformulates the ADM equations by introducing a conformal decomposition of the 3-metric and trace-free extrinsic curvature. This cures the weak hyperbolicity of the raw ADM system and yields strongly/symmetric hyperbolic evolution equations that are numerically stable for long-term simulations.

The key idea is to factor out the determinant of the 3-metric by writing γ_ij = e^4φ ˜γ_ij with det(˜γ_ij) = 1. The trace-free part of the extrinsic curvature is similarly conformally rescaled: ˜A_ij = e^-4φ(K_ij - (1/3)γ_ijK). Additionally, the conformal connection functions ˜Γⁱ = ˜γ^jk˜Γⁱ_jkare promoted to independently evolved variables, which eliminates mixed second derivatives from the principal part and renders the system strongly hyperbolic.

BSSN Evolution Equations

The Z4c extension adds constraint-damping terms proportional to a new vector Θ, which drives violations of the Hamiltonian and momentum constraints exponentially toward zero. This is essential for stable binary black hole evolutions where truncation errors continuously generate small constraint violations.

In the Z4c system, the Einstein equations are extended to a larger system that includes constraint fields as dynamical variables. Constraint-violating modes are then exponentially damped with a user-specified timescale κ₁. The McLachlan thorn in the Einstein Toolkit implements both BSSN and CCZ4 (a covariant version of Z4c) and is auto-generated from tensor expressions using the Kranc code-generation package, ensuring that the complex derivative expressions are free of transcription errors.

Constraint Equations

Video Lectures

Interactive Simulation

Finite-difference schemes on multi-block and AMR grids inevitably produce high-frequency numerical noise at mesh-refinement boundaries and near coordinate singularities. Kreiss-Oliger (KO) dissipation adds a controlled amount of artificial viscosity that selectively damps these short-wavelength modes without degrading the accuracy of the physical solution at longer wavelengths.

The fundamental problem is that the semi-discrete system (continuous time, discrete space) supports modes at the Nyquist frequency π/Δx that are poorly resolved and can be driven unstable by interpolation at refinement boundaries. Without dissipation, these modes grow exponentially and eventually contaminate the physical solution. The KO operator acts as a high-pass filter: it has negligible effect on well-resolved modes but strongly damps modes near the grid frequency.

Dissipation Operator

The order p of the dissipation operator must be chosen to exceed the order of the finite-difference stencil, so that the added dissipation converges away faster than the truncation error. Typical production NR codes use 6th-order dissipation (p = 3) with 8th-order spatial derivatives.

In the Einstein Toolkit, the Dissipation thorn applies KO dissipation to all evolved BSSN variables. The dissipation order and strength ε are set in the parameter file and can be varied per refinement level.

Interactive Simulation

The Cactus Computational Toolkit provides the software infrastructure (the “flesh”) upon which physics modules (“thorns”) are built. For numerical relativity the Einstein Toolkit assembles thorns for initial data (TwoPunctures, Lorene), evolution (McLachlan, Lean), analysis (QuasiLocalMeasures, WeylScal4), and I/O (CarpetIOHDF5).

The flesh provides a scheduler that orchestrates thorn execution in a well-defined order: initial data → analysis → evolution → boundary conditions → output, repeated each timestep. Thorns communicate through grid functions (distributed arrays) and parameters (runtime configuration). This modular design allows researchers to swap physics components without modifying the infrastructure code.

Lichnerowicz Equation

The TwoPunctures thorn solves this nonlinear elliptic equation for the conformal factor ψ using a spectral method on a compactified domain. The “puncture” approach represents each black hole as a coordinate singularity that is analytically subtracted, leaving a regular correction u solved numerically.

The conformal factor ψ encodes the gravitational field strength: it diverges as 1/r near each puncture (representing the Schwarzschild-like singularity) and approaches 1 at spatial infinity. The contour plot below shows how ψ varies on the orbital plane for different mass ratios and separations. Unequal masses produce asymmetric contour patterns, with the heavier puncture having a stronger peak.

For spinning black holes, the Bowen-York extrinsic curvature ansatz introduces angular momentum parameters S_i for each puncture. The TwoPunctures solver then finds the conformal factor consistent with both the constraint equations and the specified spins and linear momenta.

Video Lectures

Interactive Simulation

Adaptive mesh refinement (AMR) is indispensable in numerical relativity. The strong-field region near each black hole demands resolution orders of magnitude finer than the wave zone where gravitational waves are extracted. The Berger-Oliger algorithm provides a recursive time-stepping scheme where each refinement level is evolved with its own CFL-limited timestep.

A typical binary black hole simulation uses 8–10 refinement levels with a 2:1 refinement ratio. The finest grid covering each black hole has Δx ≈ M/100 (where M is the total mass), while the coarsest grid extends to r ≈ 1000M for accurate wave extraction. Without AMR, covering the entire domain at the finest resolution would require ~10¹² grid points — far beyond any supercomputer.

CFL Condition

The Carpet driver (and its GPU-accelerated successor CarpetX built on AMReX) implements vertex-centered Berger-Oliger AMR with sub-cycling in time, buffer zones for interpolation at refinement boundaries, and regridding triggered by truncation error estimates.

The recursive Berger-Oliger algorithm works as follows: advance the coarsest level by one timestep Δt₀, then recursively advance finer levels. Level l takes 2^l steps of size Δt_l = Δt₀/2^lfor each coarse step. At each sub-step, boundary data for the fine level is interpolated in both space and time from the coarse level (prolongation). After the fine level completes its sub-cycle, the coarse level receives corrections from the fine level (restriction).

The interactive simulation below shows how nested refinement regions track two punctures as they orbit. Drag the punctures to see how the AMR grid structure adapts. In a real simulation, the regridding algorithm (e.g., based on truncation error estimates or fixed boxes centered on apparent horizons) automatically determines the box positions.

Video Lectures

Interactive Simulation

The primary observable output of a numerical relativity simulation is the gravitational waveform. It is typically extracted on coordinate spheres at large radius by decomposing the Weyl scalar Ψ₄ into spin-weighted spherical harmonics. The dominant (&ell; = 2, m = 2) mode encodes most of the radiated energy, but higher harmonics become important for asymmetric systems.

The waveform encodes three distinct phases of the binary coalescence:

• Inspiral: Post-Newtonian regime with slowly increasing frequency and amplitude (“chirp”)
• Merger: The plunge and coalescence, requiring full NR; peak amplitude and frequency
• Ringdown: The remnant Kerr black hole relaxes via quasi-normal mode (QNM) oscillations, exponentially damped

For LIGO/Virgo/KAGRA data analysis, NR waveforms are interpolated into surrogate models (NRSur7dq4, SEOBNRv5) that cover the full parameter space. The simulation viewer below shows the dominant (2,2) mode with optional higher harmonics that become visible for unequal-mass or precessing systems.

Spin-Weighted Decomposition

After merger the remnant black hole rings down with quasi-normal mode (QNM) frequencies determined solely by its mass and spin. The dominant QNM frequency for the (2,2) mode of a Kerr black hole provides a direct test of the no-hair theorem.

Interactive Simulation

Before mapping numerical relativity operators onto quantum circuits, we must establish the language of quantum computation. A qubit is a two-level quantum system whose state lives on the Bloch sphere. Multi-qubit entanglement, enabled by two-qubit gates like CNOT, is the resource that potentially enables exponential speedups for structured linear algebra problems.

The computational model of quantum computing differs fundamentally from classical computing. While a classical register of n bits stores exactly one of 2ⁿ states, a quantum register of n qubits exists in a superposition of all 2ⁿ basis states simultaneously. Quantum algorithms exploit constructive and destructive interference to amplify correct answers and suppress incorrect ones. However, the output is probabilistic: measurement collapses the superposition, and the algorithm must be designed so that the desired outcome has high probability.

Qubit State Space

CNOT Gate

The CNOT (controlled-NOT) gate is a universal two-qubit gate: together with single-qubit rotations it can approximate any unitary operation to arbitrary precision (Solovay-Kitaev theorem). In the context of HHL and variational algorithms, CNOT gates are the primary source of circuit depth and noise.

For numerical relativity applications, the key quantum primitives are:

• Quantum Phase Estimation (QPE): Extracts eigenvalues of a unitary operator to n-bit precision using O(2ⁿ) controlled-U operations. Central to HHL.
• Hamiltonian Simulation: Implements e^iAt via Trotter-Suzuki decomposition or quantum signal processing. Required for QPE of the sparse Laplacian.
• Quantum Fourier Transform (QFT): O(n²) gates on n qubits; used in QPE and Shor’s algorithm.
• Amplitude Estimation: Quadratic speedup for computing expectation values, useful for extracting observables from quantum linear solver output.

Video Lectures (Qiskit Series)

The Harrow-Hassidim-Lloyd (HHL) algorithm solves the linear system Ax = b with an exponential speedup in N (the matrix dimension) compared to classical direct solvers, provided the matrix is sparse and well-conditioned. For the elliptic equations arising in numerical relativity (e.g., the Lichnerowicz equation), the condition number κ grows polynomially with grid size, which determines whether the quantum speedup survives in practice.

The HHL algorithm proceeds in three stages: (1) quantum phase estimation (QPE) to extract the eigenvalues of A encoded in ancilla qubits, (2) controlled rotation to map eigenvalues λ_jto amplitudes proportional to 1/λ_j, and (3) inverse QPE to uncompute the ancilla. The result is a quantum state |x〉 proportional to A^-1|b〉. The crucial caveat is that reading out the full classical solution vector requires O(N) measurements, erasing the exponential advantage — so HHL is useful only when one needs expectation values 〈x|M|x〉 rather than the full vector.

Complexity Comparison

Variational Quantum Linear Solvers (VQLS) offer a near-term alternative that trades the deep circuits of HHL for a hybrid classical-quantum optimization loop. VQE-style ansatze can be hardware-efficient but may suffer from barren plateaus at large qubit counts.

For the specific case of the discrete Laplacian on a 3D grid (the operator appearing in the Lichnerowicz equation), the matrix A is sparse with O(1) nonzero entries per row. Its condition number scales as κ ~ N^2/3 where N is the total number of grid points. This means the HHL complexity becomes O(N^4/3 log N / ε), which must be compared against classical multigrid at O(N). The crossover point — where HHL becomes faster — depends sensitively on the constant prefactors and the overhead of quantum error correction.

The interactive calculator below lets you explore these scaling relationships and see where (if ever) quantum algorithms gain an advantage over classical solvers for NR-relevant problem sizes.

Video Lectures

Interactive Simulation

A realistic path to quantum-accelerated NR involves hybrid workflows: the bulk of the evolution (explicit time-stepping of BSSN) remains classical, while quantum processors handle the expensive elliptic solves (initial data, gauge conditions) and potentially the implicit steps in IMEX schemes. The key question is whether current and near-term hardware can provide any advantage, or whether fault-tolerant machines with millions of physical qubits are required.

The hybrid architecture envisions a classical HPC node running the BSSN evolution loop. At each elliptic solve step (e.g., computing initial data, enforcing the conformal thin-sandwich gauge, or solving the Poisson equation in constrained evolution), the classical node packages the sparse matrix A and right-hand-side b, transmits them to a quantum co-processor, receives the solution (or expectation values thereof), and continues the classical evolution. This model mirrors the GPU-offloading paradigm already used in CarpetX.

T-gate Resource Formula

The distillation overhead C_distill for surface-code magic state factories is currently estimated at ~10³ physical qubits per logical T-gate. Combined with the κ² scaling and the precision requirement ε^-1, this makes the total physical qubit count formidable for production-scale grids.

Interactive Simulation