ForceFreeStates - BUG FIX - Power-normalized Jacobian-invariant eigenmode energies via generalized eigenproblem

[logan-nc]

Closes #304.

Problem

FreeBoundaryStability/eigenmode_energies (and the EdgeScan defaults) wrote the raw root-area-weighted quadratic-form spectrum. The congruence jarea·M†·W·M preserves eigenvalue signs (instability counts were right), but the values are energies per unit ‖b̃‖₂ — an unphysical metric under which the low end of the spectrum collapses to the numerical noise floor (DIII-D: 1.5e-15, 13 decades of dynamic range). The regression harness re-pinned the junk as expected values, and tests/benchmarks had been quietly repointed to the XiNorm/ escape hatch.

Fix

Solve the Hermitian-definite pencil W·v = λ·N·v instead of eigen(W) + post-hoc division, where N = Toeplitz(jmat)/(dV/dψ) is the power-normalization surface norm both codes already computed as a loop (ξ†Nξ = ⟨|ξ|²⟩; Fortran free.f:193-206). The eigenvalues are stationary values of the power quotient:

• power-normalized — same convention/magnitudes as the DCON-comparable values (Solovev et[1] = 0.5749, DIII-D = 1.4541), never noise;
• exactly representation-invariant — W and N transform by the same congruence under a working-coordinate change, so the b̃-space and ξ-space spectra are now the same spectrum (verified to ~1e-8); no 1/(m−nq) singularity at rational crossings.

Consistency upgrades: the Jacobian band is now a ψ-spline (ffit.jmats), so N lives on the same surface as W (ψ_lim in free_run! — a 4% fix on truncated equilibria — and the local surface at each edge-scan step); vacuum_eigenvalue uses the pencil (wv, N). All raw rootA outputs, RootAreaWeighted.jl, and the RootAreaWeighted//XiNorm/ h5 groups are removed — there is one spectrum now.

A separate commit fixes the mode-converged kinetic example deck, broken on develop since 8bf5d01 (sol.toml was referenced but never loaded): [SOL_INPUT] is embedded in gpec.toml, restoring both kinetic regression cases.

Validation

Full report with equations and 15 figures (before/after spectra, cross-Jacobian invariance of eigenvalues and physical mode shapes, mpert-convergence study, energy bookkeeping, edge-scan trajectories, regression tables): https://claude.ai/code/artifact/581b0ab5-2850-460e-bd65-87d99e27e0a2

• h5 energies reproduced from independently reconstructed N: ≤1e-9
• Pencil congruence identity (b̃ ≡ ξ spectrum): ~1e-8
• Full-run tests: 17/17
• Regression: only energy rows change (intended re-pins). DIII-D: 38/46 quantities bit-identical, including ODE steps (truncation choice unmoved), all checksums, Δ′, islands, Chirikov, and all PerturbedEquilibrium outputs. solovev_multi_n re-pins the n=2 kink from the sign-correct-but-squashed −9.0e-3 to the physical −1.906.

Scope note: cross-run agreement between different jac_type choices remains at the discretization floor of W itself (~4% bulk modes; the near-marginal Solovev leading kink is hypersensitive and not mpert-converged) — identical for the old convention, not addressable by normalization. Details and figures in report §3; suggest a follow-up issue for the near-marginal edge-kink convergence.

🤖 Generated with Claude Code

[logan-nc]

Independent assessment from the overview-paper benchmarks (Sec. 3.2.1)

While drafting the normalization/truncation section of the Julia GPEC 2.0 overview paper, we ran this PR through an independent truncation-scan study on the DIII-D-like reference case (PR #295 equilibrium; local merges of develop+#295 and develop+#295+#314; hamada/pest/boozer working coordinates; edge scan psiedge=0.85 crossing q=4, 5, 6). Summary: the pencil fix is validated and needed, but removing the b-tilde spectrum discards a property worth keeping — recommend reinstating it with null-space regularization.

What we confirmed

• Pencil math verified independently. We rebuilt N from scratch (Toeplitz Jacobian band / dV/dpsi, per-n blocks) in an external script and solved eigen(W, N) ourselves: agreement with this branch's et to all printed digits, on all three Jacobians.
• Stationarity fix is real. Cross-Jacobian spread of the least-stable eigenvalue at the final surface: 7.8% for the post-hoc convention on develop vs 4.0% for the pencil — the residual being W's own discretization floor, exactly as the PR's scope note says.
• The FreeBoundaryStability eigenmode energies written as ~0 in gpec.h5 #304 mechanism is what this PR says it is. The raw rootA spectrum's two lowest eigenvalues are metric null directions (~1e-13, eleven decades below the first physical eigenvalue at ~3.5e-4). The sqrt-area-weight operator is nearly rank deficient in its high-|m| tail. A global rescale cannot fix an in-spectrum gap; along an edge scan the raw rootA "least-stable eigenvalue" is pinned to these nulls in every stable window.

What we found that argues against full removal

Projecting out the metric null directions (threshold 1e-6 of the spectral radius) and tracking the least physical b-tilde eigenvalue along the same scans:

• The regularized field-norm margin is nearly flat (~10%) across each valid truncation window (the stable regions just outside rational crossings, where production truncations live via dmlim), falling only as the next rational is approached.
• The displacement-norm (pencil) margin swings by more than an order of magnitude across the same windows (0.04 -> 1.49 -> 0.63 in the q=4..5 window), because the near-resonant displacement content it weights varies strongly between rationals.
• Cross-Jacobian agreement of the regularized b-tilde margin is ~10% (small-eigenvalue discretization amplification); signs agree with the pencil everywhere (Sylvester).

So the two normalizations answer different questions: the pencil is the correct variational form of the DCON convention (comparability, tight invariance), while the regularized b-tilde margin is truncation-robust — a stability margin quoted at the standard truncation is insensitive to the exact truncation choice. That robustness is a core motivation of the normalization work and is lost if the b-tilde spectrum is deleted rather than repaired. NaN at exact rational crossings is fine (the metric is singular there and those truncations are invalid anyway).

[figure: least-stable eigenvalue vs truncation surface, panel (a) pencil / panel (b) regularized b-tilde — attached below]

Actionable suggestions (scoped for a follow-up session)

• Reinstate the b-tilde spectrum as a secondary output with null-space regularization: SVD (or eigendecomposition) of the sqrt-area-weight operator, project out directions below rtol * sigma_max (default rtol = 1e-6, exposed as a control), report the spectrum of the projected congruence.
• Name the h5 groups unambiguously, e.g. FreeBoundaryStability/eigenmode_energies (pencil, as in this PR) and FreeBoundaryStability/btilde_energies (regularized), so no "default vs escape hatch" ambiguity recurs.
• Add the regularized b-tilde least eigenvalue to the EdgeScan outputs alongside the pencil trajectory (both are cheap; the b-tilde trajectory is the truncation-robustness diagnostic).
• Keep NaN at exact rational crossings but state it in the EdgeScanState docstring/contract (valid truncations sit dmlim outside rationals, as in the Fortran lineage).
• Pin a regression quantity on the regularized b-tilde lambda_1 at psilim for the DIII-D case, so a null-direction bug of the FreeBoundaryStability eigenmode energies written as ~0 in gpec.h5 #304 class cannot be re-pinned silently again.
• Add a cross-Jacobian invariance test (hamada/pest/boozer) asserting the pencil spectrum to solver tolerance and the regularized b-tilde spectrum to the discretization floor, on the Solovev case for speed.
• Reproduction: figures/make_truncation_scan.jl, figures/make_rootA_flatness_scan.jl + committed CSVs in the overview-paper repo. The flatness scan needs this small collector hook in free_compute_total (Free.jl, on develop-era code), since transform_u! rewrites u_store after the edge scan:

# after: sqrtamat_local/jarea_local evaluated from es.sqrtamat_spline
rootA_result = compute_rootarea_eigenvalues(wt_saved, wp, wv, sqrtamat_local,
    jarea_local, equil, odet.psifac, intr; all_eigenvalues=true)
if isdefined(Main, :ROOTA_COLLECTOR)
    push!(Main.ROOTA_COLLECTOR, (odet.psifac, sort(real.(rootA_result.rootA_et_all))))
end

Numbers above were measured on an Apple M-series laptop with julia -t 8; eigenvalue ratios and spreads are hardware independent.