# BSSN / Z4c Validation Results (vacuum, GMDBS-toroidal substrate)

Date: 2026-05-15
Solver: `wasm/src/bssn/` - 4th-order central finite differences, periodic
boundaries (T^3 topology), 4th-order RK4 time stepping, harmonic and 1+log
slicing, gamma-driver shift, Kreiss-Oliger 6th-order dissipation, algebraic
gamma-tilde determinant pullback. Two formulations select at run time:
`bssn` (pure vacuum BSSN, optional kappa_h Hamiltonian damping) and `z4c`
(Theta auxiliary plus Hilditch constraint damping).

Every number below is produced by `cargo test --release` with the named
test. Per-run JSON and CSV land in `_diag/bssn_*`. No number in this
document is hand-entered; they are copied from the test output.

## Summary

| Gate | Test | Observed | Reference |
|---|---|---|---|
| Convergence | factor-2 self-convergence on gamma_xx | 3.946 | 4.0 (4th-order FD) |
| Scaling | wall-time slope vs N | 3.048 | 3.0 (O(N) explicit RK4) |
| BSSN long-time | bounded gauge-wave evolution | 5 light crossings | 5-10 (classical BSSN) |
| Z4c long-time | bounded gauge-wave evolution | 7 light crossings | extends past BSSN |

Convergence and scaling match a 4th-order vacuum BSSN code on the
Apples-with-Apples gauge-wave testbed (Alcubierre 2004; Brown 2009;
Bernuzzi & Hilditch 2010). The BSSN bounded horizon and the Z4c extension
of that horizon are discussed in sections 4 and 5.

## 1. Convergence

Test: `bssn_convergence::gauge_wave_self_convergence_at_or_above_third_order`

Factor-2 self-convergence (Choptuik 1991), the community-standard NR
convergence measure. Apples-with-Apples gauge wave to t = 1.0 at four
resolutions; differences between successive doublings are restricted to
the coarsest grid by integer-stride sampling; observed order is
log2(d_N / d_2N).

```
d(N=12, 2N=24) = 3.558e-6
d(N=16, 2N=32) = 1.151e-6
d(N=24, 2N=48) = 2.310e-7
observed order = 3.946     (theoretical 4)
```

Self-convergence rather than analytic-solution comparison: the
gauge-wave-vs-exact error plateaus near 1e-3 across all resolutions
because the periodic-substrate algebraic-projection step has its own
resolution-independent floor. Self-convergence isolates the
finite-difference order. This is the same diagnostic Cactus, GRChombo,
and GRAthena++ report. Dissipation is set to zero for this test so the
measured order reflects the bare stencil.

Output: `_diag/bssn_convergence.json`

## 2. Scaling

Test: `bssn_scaling::wall_time_scales_as_n_cubed`

Wall time vs grid size N at fixed n_steps = 50, linear fit on
log(wall) vs log(N).

```
N=12   cells=  1,728   wall=  2.05s   23.73 us/cell/step
N=16   cells=  4,096   wall=  4.85s   23.66
N=24   cells= 13,824   wall= 16.54s   23.94
N=32   cells= 32,768   wall= 41.10s   25.09
N=48   cells=110,592   wall=138.68s   25.08
observed slope = 3.048     (theoretical 3 for O(N) explicit RK4)
```

CHE billing meter calibration: 23.94 us/cell/step (median). The BSSN
tab's CHE meter converts wall time and resolution into Computed-Core-Hour
Equivalents against this calibrated reference.

Output: `_diag/bssn_scaling.json`

## 3. Smoke gates

Test: `bssn_validation`

Three sanity checks the production defaults must hold:

- Minkowski is a fixed point. Hamiltonian L2 stays below 1e-6 and
  alpha, chi stay within 1% of unity for every snapshot.
- Small-amplitude gauge wave does not blow up over 100 steps; lapse and
  conformal factor stay positive.
- Default 16^3 x 100-step run finishes well under the native budget so
  the in-browser run stays under the 3 s target.

All three pass.

## 4. BSSN long-time bound and the constraint-mode gap

Pure vacuum BSSN with kappa_h Hamiltonian damping stays bounded for
about 5 light crossings of the gauge wave (N=16, harmonic slicing). Past
that the canonical formulation-level constraint-violating eigenmode
appears. This is not a bug and not under-resolution; it is the
limitation of free-evolution BSSN that motivated the Z4 family.

Diagnostic test: `bssn_spectral_diag::z_resolution_sensitivity`. Hold the
transverse resolution fixed (Nx=Ny=16) and vary nz over {8, 16, 32},
recording the blow-up time (first |H| > 1e-2). Pure BSSN, harmonic
slicing, no Theta:

```
(16,16, 8)   blow-up = 60.00   max_ham = 2.70e9
(16,16,16)   blow-up = 60.00   max_ham = 2.14e4
(16,16,32)   blow-up = 60.55   max_ham = 1.89e2
```

The blow-up time is invariant under z-refinement; only the amplitude at
blow-up shrinks. That is the signature of a formulation-level
constraint-violating mode, not numerical aliasing (aliasing would push
the blow-up time out as nz grows). Adding kappa_h damping suppresses the
dominant mode by roughly 6x but a residual subspace remains. The fix is
Z4c, which promotes the Hamiltonian and momentum constraints to
dynamical fields with their own propagation and damping.

Output: `_diag/bssn_z_resolution.json`

## 5. Z4c long-time evolution

Z4c is now implemented in the CPU solver and on the GPU path. It adds one
scalar field Theta (in `state.rs`, alongside chi and K), the
kappa_1(1-kappa_2)Theta term in d_t K, the -2 alpha kappa_1 Z-hat^i term
in d_t Gamma-tilde^i, and the d_t Theta evolution equation in `rhs.rs`.
It is selected with `formulation = "z4c"`; kappa_1 and kappa_2 are the
Bernuzzi-Hilditch damping parameters. When Theta is zero and the kappa
terms are off, Z4c reduces to BSSN and reproduces sections 1-3 exactly.

Test: `bssn_long_time::long_time_gauge_wave_stays_bounded`
Config: N=16, gauge wave, harmonic slicing, KO=0.3, algebraic
pullback=0.05, formulation=z4c, kappa_1=0.02, kappa_2=0. Box L=10,
wavelength=L, amplitude 0.01. Evolved to t=70, which is 7 light crossings
(one crossing = L/c = 10 time units).

```
initial ham_l2     = 2.594e-5
max ham_l2         = 3.726e-5      (1.44x initial)
late/early ratio   = stable equilibrium, not exponential
max mom_l2         = 3.300e-4
min alpha          = 0.9945
min chi            = 0.9967
max |det gamma-tilde - 1| = 7.77e-16   (machine precision)
n_steps = 448      wall = 45.9 s
```

The constraint norm grows by 1.44x over the full run and then sits at a
bounded equilibrium rather than growing exponentially. Lapse and
conformal factor stay healthy, and the algebraic determinant constraint
holds at floating-point precision. Z4c extends the bounded horizon from
BSSN's roughly 5 crossings to 7 crossings at these parameters.

Scope of this result: this is bounded vacuum evolution to 7 crossings,
not the full 1000-crossing Apples-with-Apples horizon. Reaching that
horizon needs a kappa_1/kappa_2 parameter sweep (the Adv-tier batch
driver) and completing the Z4c equation set; a few terms in d_t
A-tilde and d_t chi are simplified in the current implementation. That
is scoped next-phase work, tracked in
[bssn_z4c_pathway.md](bssn_z4c_pathway.md).

Output: `_diag/bssn_long_time.json`, `_diag/bssn_long_time.csv`

## 6. What this solver is for

Short-to-medium-time vacuum GR evolution on T^3-periodic substrates:
about 5 crossings in pure BSSN, about 7 with Z4c at the validated
parameters. Convergence studies of new BSSN-related operators (4th order
verified). Resolution-to-runtime characterization for capacity planning
and CHE calibration. Educational and demo runs of the gauge-wave
testbed, including the in-browser run at the demo defaults (12^3, 60
steps, about 1.3 s native).

## 7. What it is not yet for

Full 1000-crossing long-time evolutions and binary-merger inspiral:
these need the completed Z4c equation set plus a damping-parameter
sweep. Strong-field collapse and singularity formation: need puncture
initial data, an apparent-horizon finder, and Z4c. Waveform extraction
at infinity: needs psi4 extraction and outer-boundary handling beyond
pure periodic. Matter coupling: this is a vacuum-only solver; no T_munu
source terms.

## 8. Reproducibility

```
cd wasm
cargo test --release --test bssn_convergence
cargo test --release --test bssn_scaling
cargo test --release --test bssn_validation
cargo test --release --test bssn_long_time -- --ignored
cargo test --release --test bssn_spectral_diag z_resolution_sensitivity -- --ignored
```

All numbers above reproduce verbatim. JSON outputs in `_diag/bssn_*.json`.

## 9. References

- Alcubierre, M. et al. Towards a stable numerical evolution of strongly
  gravitating systems in general relativity: the conformal treatments.
  Phys. Rev. D 62, 044034 (2000).
- Brown, J. D. Covariant formulations of BSSN and the standard gauge.
  Phys. Rev. D 79, 104029 (2009).
- Bernuzzi, S. & Hilditch, D. Constraint violation in free evolution
  schemes: comparing BSSNOK with a conformal decomposition of Z4.
  Phys. Rev. D 81, 084003 (2010).
- Hilditch, D. et al. Compact binary evolutions with the Z4c
  formulation. Phys. Rev. D 88, 084057 (2013).
- Choptuik, M. W. Consistency of finite-difference solutions of
  Einstein's equations. Phys. Rev. D 44, 3124 (1991).
- Alcubierre, M. et al. Towards standard testbeds for numerical
  relativity. Class. Quantum Grav. 21, 589 (2004). The
  Apples-with-Apples testbed paper.