ForceFreeStates - BUGFIX! - Use trace-aware Frobenius exponents at resonant surfaces

[jhalpern30]

Release note

• Audience: users
• Numerical impact: imaginary Δ′ at resonant surfaces drops from the spurious large value (up to 1e4) to roundoff-level on all but the near-axis surface of high-n decks. Real Δ′ and the free-boundary energies move ≤ 3×10⁻⁵ relative on the verified decks. Galerkin Δ′ real parts move ≤ 6×10⁻⁴ (q=5), and its imaginary parts shrink but are not at roundoff.
• Migration: none. A leftover extended_precision_bvp key is dropped with a warning.
• Regression harness: run at be852d436 against its merge-base eb8e5f1bc (14 cases, report below). The later commit only removes dead code left by the flag removal (no numerical change; Riccati tests pass).

Resonant Frobenius powers now come from the eigenvalues of each side's M₀ block. A traceless block is bitwise the same as before. extended_precision_bvp is removed: with the exponent fix the Float64 BVP solve is accurate, and it is 3.5–5.4× faster on large matrices.

Why this is a bugfix, and how it relates to Fortran

What Fortran does. Every Fortran singular-surface routine takes the resonant 2×2 block of M₀ and sets di = m0(1,1)*m0(2,2) − m0(2,1)*m0(1,2), then alpha = SQRT(−di) (dcon/sing.f:305, rdcon/sing.f:297, rdcon/sing1.f:491, stride/sing.F:395). It uses ±alpha as the two Frobenius powers (stride/sing.F:410), builds the zeroth-order eigenvectors as −(m0(1,1) ± sig·alpha)/m0(1,2) (stride/sing.F:447), and rescales the columns by dpsi**alpha (stride/sing.F:880). It computes alpha once, on the right side, and reuses it for the left. The Julia code previously reproduced this exactly, including the shared right-side alpha (alpha_override).

Why it is wrong. The solutions near a rational surface behave as |ψ − ψₛ|^λ, where λ solves det(M₀ − λI) = 0, so λ = tr/2 ∓ √(tr²/4 − det). The Fortran value ±√(−det) equals those roots only when tr = 0. It is a shortcut that holds for a traceless block and is silently assumed otherwise. When tr ≠ 0 both columns carry the wrong power (off by tr/2) and the zeroth-order eigenvectors are wrong at O(tr). The error is amplified by the matching distance, because Δ′ is read off as a coefficient that is |dpsi|^(2α) smaller than the dominant solution.

Why this is physical, not a tuning choice. For an ideal outer region with a real equilibrium, the Δ′ matrix must be real. With Fortran's exponents the old code returned a large imaginary part (DIII-D n=1, q=5: 67; Solov'ev q₀=1.85, q=2: −1.8×10⁹), and with the eigenvalue exponents it falls to roundoff (1×10⁻⁵ and −27, the latter against Re Δ′ of 3×10⁸). Real Δ′ barely moves (≤ 3×10⁻⁵ relative). The imaginary part is the symptom, and the mismatch with the block's eigenvalues is the cause. Mercier D_I comes from the local calculation and is unchanged; the asymptotics α changes by only O(tr²) (≤ 1.1×10⁻⁷ relative in the Galerkin case).

Divergence from Fortran. This is a deliberate divergence on non-traceless blocks. Where the block is traceless the output is bitwise identical to Fortran's formula (resonant_block_exponents keeps the old √(−det) operations on that path, and a unit test checks it). The same shortcut exists in the upstream Fortran. I have not yet checked the choice against the Glasser (2016) asymptotic-series derivation line by line; a physics reviewer should. Left unchanged: the shared-alpha reuse is removed only because each side now computes its own block (√(tr²/4 − det) is the same on both sides, so Mercier α is not affected).

Bugfix, not a feature. There is no new option or capability: the same quantity (Frobenius powers at a resonant surface) is computed as the math defines it, and a leftover extended_precision_bvp key is only dropped. It carries ! because results move (imaginary Δ′ and the 1e-9 to 1e-5 shifts in energies below), not because behaviour is added.

Regression report

regress --refs eb8e5f1bc,be852d436 on Feynman, 14 cases (the two ggj_* cases are inner-layer only and skipped). All ran; nothing failed. Compared against the merge-base rather than current develop, which has since gained #514 (axis start) and #496; rerun the harness on the merged head (origin has since merged develop) before merging.

Case	What moved
diiid_n1, diiid_multi_n, solovev_n1, solovev_multi_n	Energies 1e-9 to 1e-5 relative; forward ODE step count shifts by 1–7 (≤ 1%); island widths and Chirikov ≤ 0.03%
diiid_n1_riccati	Re Δ′ ≤ 2×10⁻⁷ relative. The 4.1% flagged is the imaginary part (92 → −4×10⁻⁶ at q=5)
gal_resistive_diiid	Δ′ matrix norm 12437 → 2184, from the spurious imaginary part. Re diagonal ≤ 6×10⁻⁴. Residual Im is 0.30 at q=4 (was 0.013) and 23.5 at q=5 (was −15173), so Galerkin Im is not at roundoff
diiid_slayer_n1, solovev_kinetic_ntv, solovev_kinetic_multiion, diiid_error_field	≤ 0.06%
gal_resistive_pe, efit_fixedbdy_separatrix, solovev_kinetic_calculated, solovev_kinetic_nuzero	unchanged

Notes for reviewers

alpha_override is gone. Each side computes its own block. D_I still uses α, not the column power.

Removing extended_precision_bvp

Same-machine Double64 baselines against the final code, DIII-D n=1,2 and Solov'ev, at singfac_min = 1e-4: resonant diagonals agree to ≤ 2×10⁻¹³ and the edge energy is bitwise identical. The largest production discrepancy is one coupling of 56 on Solov'ev q₀ = 1.85 (relative 8×10⁻⁶). singfac_min = 1e-6 at the same q₀ misses that coupling by 2%, which is a matching distance outside the recommended range.

Larger matrices (DIII-D, Riccati, 4 threads, same deck, base uses Double64 and head Float64):

Case	Force-free stage	Total wall
n=2, 10×10	104 s → 68 s (1.5×)	135 s → 98 s
n=3, 15×15	272 s → 77 s (3.5×)	303 s → 108 s (2.8×)
n=4, 20×20	735 s → 136 s (5.4×)	765 s → 167 s (4.6×)

At n=4 the base with the flag off takes 119 s, so the whole speedup is the Double64 solve. Double64 and Float64 agree to 1.7×10⁻⁶ on the worst entry of the base 20×20 matrix, so precision does not drive the differences below.

Open question: near-axis surface at high n

On DIII-D n=4 the first surface (q=1.25, |Δ′| ≈ 2300) moves by 5×10⁻³ on its diagonal and by O(1) on its couplings to other surfaces (one entry 0.52−1.75i → −1.58+3.06i, about 1% of √(Δ′ᵢᵢΔ′ⱼⱼ)). This comes from the exponents, not the precision change. The other 19 diagonals move ≤ 8×10⁻⁶, and the head still shows Im Δ′ = −25 at that surface against ~10⁻⁴ elsewhere, so it is not clean there. I could not establish which value is closer to the truth.

DoubleFloats stays. The inner-layer GGJ march still uses it.

[jhalpern30]

@matt-pharr and @d-burg the start of my Delta prime investigation. This one seemed super promising and allowed removing of that extended_precision_bvp flag and large speedups now for higher n. The numbers seem to checkout, and the physics also seems legit, but I think it needs an expert eye to see if what I am doing here is garbage.

Not sure which one of you was most applicable to look at this so I requested both for review - maybe @d-burg since the git blame says you added this extended precision thing?

[matt-pharr]

Ah @jhalpern30 this is what I just messaged Daniel about this morning and which I mentioned in our DMs lol. I think I wanted to discuss if this makes Delta' more right or just more self-consistent. Glasser says the eigenvalues of $M_0$ should be $\pm\sqrt{-D_I}$, but that this depends on Fourier mode convergence. So should non-hermiticity in a resulting Delta' be a sign of non-convergence or no? I am leaning toward what you implemented here, as it seems important to get the series expansion right for the asymptotics. It also strongly affects other RDCON results because using Glasser's "eigenvalues" results in larger imaginary components on the diagonals of the Galerkin matrix which should be hermitian, and makes it so a Cholesky matrix solve is incorrect because it discards this quantity, so I was working on this in a different branch. I can move that over to here though.

[matt-pharr]

Here is what my claude has to say in addition:

I hit the same bug independently today while chasing a Galerkin LU/Cholesky gap on #513. The implementation and the numbers match mine. A few notes:

1. "Δ′ must be real" is too strong. That holds only for up-down symmetric equilibria; the general constraint is Hermitian with a real diagonal. The DIII-D-like case is asymmetric. With this fix its diagonal is real, but the off-diagonals are genuinely complex: the loop phases are −26° to −84°, and no choice of per-surface phases removes them. LAR comes out real. Suggest rewording to "Hermitian (real diagonal)".
2. Δ′ is still non-Hermitian by a per-surface scaling. The scaling is the ratio of each surface's big/small cross-Wronskian W. Normalizing each pair by √|W| on top of this fix makes Δ′ directly Hermitian: 5e-9 on DIII-D Riccati and about 1e-9 on LAR. That's what @d-burg's PrincetonUniversity/GPEC#191 was after, so it would make a good small follow-up.
3. The residual Galerkin Im is expected. q=5 is under-resolved at mpsi=128: mpsi=256 moves Delta_coil there by about 90%. There's also a small F̄ spline mismatch between the series and the assembly. It isn't an exponent issue, and with this fix Galerkin no longer depends on gal_dx0.
4. #513: with this fix, Galerkin LU and Cholesky agree to 2e-5 (was 6e-3), so a Cholesky-first default would be essentially results-neutral.

[jhalpern30]

@matt-pharr ohh I see, glad to see out Claude's both independently verified it. If you're ok with it, I'll assign you to this PR to merge if you want it since it sounds like it affects your work and I don't want to mess up what you're doing