Numerical Generalized-Harmonic Evolution

Generalized-harmonic evolution promotes coordinate choice to gauge-source functions and rewrites Einstein’s equations as coupled wave equations. Numerically, its appeal is a characteristic structure suited to constraint-preserving boundaries and black-hole excision. Its danger is equally clear: inconsistent gauge sources, reduction constraints, damping parameters, or incoming boundary fields can drive a stable-looking metric away from the intended Einstein solution. This page treats the discrete system and its tests; continuum hyperbolicity and gauge theory remain with Formulations of General Relativity.

Helpful background. Formulations of General Relativity develops harmonic formulations, Ordinary Differential Equations supplies time integration and stability, Computational Relativity supplies evidence and rerun conventions, and Black-Hole Dynamics and Collapse supplies the strong-field applications.

Gauge and first-order fields

Define the harmonic constraint

$$C_a=H_a+\Gamma_a,$$

where $H_a$ is the declared gauge source and $\Gamma_a=g^{bc}\Gamma_{abc}$ under a fixed sign convention. A common first-order reduction evolves $g_{ab}$ together with spatial derivatives $\Phi_{iab}=\partial_i g_{ab}$ and a normal-time derivative $\Pi_{ab}=-n^c\partial_cg_{ab}$. It introduces the reduction constraint $C_{iab}=\partial_i g_{ab}-\Phi_{iab}$ and the curl or integrability constraint

$$C_{ijab}=2\partial_{[i}\Phi_{j]ab}=2\partial_{[j}C_{i]ab},$$

in addition to $C_a$. First-order implementations may also evolve or diagnose derived constraints customarily denoted $F_a$ and $C_{ia}$; the exact complete constraint set must be stated rather than reduced to one representative norm.

The implementation record must state whether $H_a$ is analytic, frozen, damped-harmonic, or evolved; how its derivatives are obtained; and the signs and dimensions of damping parameters. Lindblom et al. 2006, §2.2 and §§3–5 constructs a first-order symmetric-hyperbolic system with explicit constraint damping and boundary fields. Translating it incompletely—for example damping $C_a$ but not controlling reduction constraints—defines a different algorithm.

Damping and discretization work together

Parameters commonly denoted $\gamma_0$, $\gamma_1$, and $\gamma_2$ control harmonic-constraint damping, characteristic linear degeneracy, and reduction-constraint damping. In coordinates whose time and spatial coordinates have length units, $\gamma_0$ and $\gamma_2$ have inverse-length units, whereas $\gamma_1$ is dimensionless. Their useful choices depend on local scales and the spatial discretization. An excessively large $\gamma_0$ can make explicit time integration stiff; too small a value allows long-wavelength contamination to persist. Damping does not remove every constraint mode, especially spatially constant modes outside the assumptions of a Fourier estimate.

For the first-order system of Lindblom et al. 2006, §3, Eqs. (32)–(34), define an outward spatial unit normal $s_i$ and $\beta^s=\beta^is_i$. The characteristic combinations and coordinate speeds are

$$\begin{array}{c|c} \text{field}&\text{normal speed}\\ \hline g_{ab}&-(1+\gamma_1)\beta^s\\ \Pi_{ab}\pm s^i\Phi_{iab}-\gamma_2g_{ab}&-\beta^s\pm\alpha\\ (\delta_i{}^j-s_i s^j)\Phi_{jab}&-\beta^s \end{array}$$

Negative speeds enter an outer boundary in this orientation. With $\alpha=1$, $\beta^s=0$, the minus field enters and the plus field leaves; zero-speed fields are handled by the evolution equations. A boundary moving normally at coordinate speed $v_b^s$ changes the classification to $v-v_b^s$. Check this sign with a translated flat-space wave before applying excision. The gauge identity $\Box x^a=-\Gamma^a=H^a$ also checks the plus sign in $C_a=H_a+\Gamma_a$.

Finite-difference, finite-volume, and discontinuous-Galerkin realizations differ in numerical flux and penalty choices. Record characteristic speeds, interface fluxes, dissipation or filtering, time integrator, CFL factor, and any projection. Re-evaluate constraints with an independent derivative operator and output the damping source contribution separately from the undamped evolution residual.

Boundaries and a benchmark ladder

At an outer boundary, decompose fields into incoming and outgoing characteristic combinations. Supply physical radiation data only for the physical incoming fields, enforce harmonic and reduction constraints on their incoming sectors, and choose gauge data consistently. Freezing every incoming component is generally reflective and need not preserve constraints; Rinne 2006, §§2.1–2.3 and §§3–5 derives and tests radiation-controlling, constraint-preserving conditions for this characteristic system.

Begin with robust stability, linear-wave, gauge-wave, and shifted-gauge-wave tests. Compare runs with $\gamma_0$ and $\gamma_2$ enabled and disabled, holding all other controls fixed. The damped run should suppress the targeted finite-wavelength constraint contamination without changing the convergent physical wave. Lindblom et al. 2006, Figures 1–3 demonstrates constraint behavior and black-hole evolutions under the stated formulation. Preserve the undamped control run even when it fails, because its growth rate tests whether the observed improvement follows the predicted damping mechanism rather than unrelated filtering.

Then evolve stationary Schwarzschild or Kerr data. Monitor $C_a$, $C_{iab}$, $C_{ijab}$, the formulation’s derived first-order constraints, horizon mass and spin, curvature invariants, and spurious radiation. Compare constraint-preserving and frozen outer boundaries at several radii. A reflected pulse should move in arrival time with boundary radius; physical stationary invariants should converge independently of that placement.

Use at least three resolutions and vary damping parameters around the production values. Convergence of the metric with fixed damping does not establish robustness to an unresolved damping length scale. The method’s domain is the tested gauge and formulation; success for stationary excision does not automatically transfer to matter shocks or moving-puncture interiors.

Environment and failure criteria

Pin the exact generalized-harmonic equations, $H_a$ prescription, damping profiles, first- or second-order reduction, complete constraint set, numerical flux or stencil, filters, time integrator, CFL factor, boundary characteristic treatment, mesh hierarchy, precision, and code revision. A clean interface could be evolve-gh --config stationary-kerr.toml --output run/, emitting all constraints, boundary fluxes, damping terms, invariants, resolution and parameter sweeps, and environment metadata. It specifies a rerun record, not an executed artifact on this page.

Reject an evolution if harmonic constraints grow with resolution, physical fields change materially when damping is disabled after convergence, incoming boundary contamination appears, or invariant horizon quantities drift behind stable coordinates. Failure modes include a sign error in $H_a$, inconsistent derivative reduction, damping-induced stiffness, an excision surface with incoming characteristics, and a boundary algorithm that controls $C_a$ while injecting $C_{iab}$ or $C_{ijab}$.

Diagnostic

Increasing $\gamma_0$ by a factor of ten lowers $\lVert C_a\rVert$ by two orders of magnitude, but forces the time step down by twenty and shifts the horizon mass by $10^{-3}$ without resolution convergence. Is stronger damping better?

Solution

No. Constraint damping is part of the numerical system, not an objective to maximize. The mass shift and stiffness show that the chosen scale is unresolved or is altering the computed solution. Perform a joint resolution, time-step, and damping sweep; require invariant quantities to converge while the intended finite-wavelength constraints decay, and check reduction and boundary constraints separately.

References

• Lindblom, L., Scheel, M. A., Kidder, L. E., Owen, R., and Rinne, O. (2006). “A new generalized harmonic evolution system.” Classical and Quantum Gravity 23, S447–S462. arXiv.
• Rinne, O. (2006). “Stable radiation-controlling boundary conditions for the generalized harmonic Einstein equations.” Classical and Quantum Gravity 23, 6275–6300. arXiv.