PE - BUG - Im(py) is ~900x larger than Fortran GPEC, producing a spurious toroidal torque #401

Split out from #370, where the symptom (a toroidal torque that should be zero) was traced to its
source.

The measurement

toroidal_torque = -2·n·Im(py) where py = ⟨Φ_resp, Λ⁻¹·Φ_resp⟩/4 — the imaginary part of a
quadratic form, identically zero for Hermitian Λ. Julia gives ~0.05 on the DIII-D-like ideal deck.

Fortran GPEC does not reproduce it. Rebuilt from the exact checkout (v1.5.7-611-ge68d7ac2,
gfortran 15.2), running docs/examples/run_ideal_example with inputs unmodified. The decks match:
byte-identical TkMkr_D3Dlike_Hmode.geqdsk, same d3d C-coil currents
[1000, 500, −500, −1000, −500, 500], hamada, nn=1, singfac_min=1e-4, and resp_index=0 —
the energy-identity method, which is the one Julia implements.

The real part of the same complex number agrees to 8 parts in 100,000, so the chain feeding it
(equilibrium, DCON solve, coil forcing, inductances, response flux) is faithfully ported. The
imaginary part is off by ~900×. An older Fortran build (v1.5.5-323, ~290 commits back) gives
8.1e-5, so the Fortran value is stable and ~1e-5 relative — zero to numerical precision, as the
Hermitian argument requires.

What is already ruled out

Not a discretization error. |T| is refinement-insensitive across five axes: mpsi (128→2048),
mtheta (256→512), mthvac (512→2048), eulerlagrange_tolerance (1e-10→1e-12), and
singfac_min (1e-3→1e-6). The last of these converges — first-order, each decade cutting the
change ~11× — onto a finite nonzero 5.0758e-2 rather than vanishing. Details in #370.

Not a response-method mismatch: Fortran runs resp_index=0, the method Julia ports.

Not the plasma-inductance assembly formula. Fortran gpresp.f:201-213 and Julia
ResponseMatrices.jl:256-264 are algebraically identical — Fortran's
t1 = -1/(chi1·s_i·2π·i) equals Julia's t1 = im/(chi1·s[i]·2π), likewise t2, and both form
2·t1·wt0[i,j]·t2 and invert it.

Lead: the surface inductance is not symmetrized

Fortran symmetrizes it explicitly before use (gpec/gpresp.f:284):

temp1=0.5*(temp1+CONJG(TRANSPOSE(temp1)))
surf_indmats=temp1

Julia's calc_surface_inductance has no equivalent step — it returns
inv(I_v) .* (μ₀ * (2π)^2) (ResponseMatrices.jl:306).

py is built from Λ (plasma), not the surface inductance directly, but response_flux is obtained
through the permeability and reluctance, both of which are built from surface_inductance
(Response.jl:42,47). An anti-Hermitian component surviving there propagates into Φ_resp and so
into Im(py), while leaving Re(py) essentially untouched — which is the signature observed.

Worth checking alongside: Response.jl:45 asserts "Λ (plasma inductance) is not Hermitian — its
anti-Hermitian part is the dissipative/torque response". For an ideal run with no dissipation that
part should vanish, so it is worth confirming that comment describes the intended physics here
rather than licensing the residual.

I have not verified that symmetrizing fixes the number — flagging it as the most specific divergence
found, not a diagnosed cause.

Reproducing

Julia: examples/DIIID-like_ideal_example as shipped; read
PerturbedEquilibrium/Energies/{toroidal_torque, plasma_energy} from gpec.h5.

Fortran: docs/examples/run_ideal_example with dcon then gpec; the values are printed as
"Required energy/torque to perturb plasma" and written as toroidal_torque (gpout.f:1353).

Consequence for regression testing

The torque is pinned in the golden-value work (#397) as a bound, |T| < atol with target 0, rather
than as a value — pinning a number that should be zero would encode this as truth. Fortran's ~1e-5
relative is the scale that bound should eventually reflect.