Local-stability scan: numerical instability near core (ψ < 1e-2) when q′ → 0

[logan-nc]

Summary

The local-stability scan (local_stability_flag = true, src/ForceFreeStates/Bal.jl) shows
numerical instability near the magnetic axis — roughly ψ < 1e-2 — when the equilibrium is run
with psilow = 1e-4. The affected outputs are the per-surface profiles locstab/di (Mercier
D_I), locstab/dr (resistive interchange D_R), and locstab/ballooning_Delta_prime. Outside
the core (ψ ≳ 1e-2) the profiles are smooth and well-behaved.

This was observed by @jmlmir369 during review of #234 and is split out here per that discussion.

Root cause

The core instability is traced to the magnetic shear q′ → 0 in the core region. With small
q′, several coupled terms in the local-stability formulas become ill-conditioned:

• Resistive interchange D_R has q′ (q1) in a denominator: Bal.jl:122, factor
  1/(q1·chi1³).
• The ballooning sigma term (Bal.jl:489) and the m0_12 / d0bar → det Mercier cascade
  (Bal.jl:492, Bal.jl:515) inherit the small-shear conditioning.

This is not new to the rewrite — the same behavior existed in the original standalone
Mercier.jl, so it is a pre-existing property of the local-stability formulation, surfaced now
that the scan runs down to psilow = 1e-4 during testing.

Impact — benign for the global result

This affects only the local (high-n) diagnostic profiles near the axis. It does not affect
the global ideal/perturbed-equilibrium result, which is the primary purpose of the code: the
global path keeps psilow = 1e-4 (so the Frobenius initial conditions get close to the axis) and
remains bit-identical in regression. The core noise is confined to the locstab/* diagnostics and
does not propagate outward (confirmed in @jmlmir369's core-vs-outer comparison plots). The merged
examples therefore keep psilow = 1e-4.

Why not a quick fix

A blunt mitigation (flooring the shear in the denominators, or skipping all surfaces below a
shear/ψ threshold) is unsatisfactory: q′ → 0 near the core is physical (low central shear), so
a floor distorts the interchange criteria and a core-skip drops genuine interior surfaces. A proper
fix needs an analytic small-shear limit / regularization of the sigma–m0_12–det(d0bar)
chain, which is a deeper change than this PR's scope.

Suggested follow-up

• Derive the q′ → 0 asymptotic limit of D_I, D_R, and the ballooning Δ' and apply it as a
  regularized branch in Bal.jl for low-shear surfaces.
• Add a regression check on the core-region locstab/* profiles once a regularization exists.

Reference image

@jmlmir369's comparison of the determinant / resistive-interchange D_I with profile derivatives,
split into core (1e-4–1e-2) and outer (1e-2–0.99) regions:

Split out from #234.