diff --git a/benchmarks/benchmark_solovev_kinetic_stability.jl b/benchmarks/benchmark_solovev_kinetic_stability.jl
index ae98fe02..3ced63df 100644
--- a/benchmarks/benchmark_solovev_kinetic_stability.jl
+++ b/benchmarks/benchmark_solovev_kinetic_stability.jl
@@ -14,7 +14,7 @@ band (`delta_mlow=delta_mhigh=0`), same `nl`, collision operator, grid, and
 edge truncation (`sas_flag`/`dmlim`).
 
 The comparison is the least-stable total-energy eigenvalue:
-  - Julia:   `FreeBoundaryStability/XiNorm/eigenmode_energies[1]` in `gpec.h5`
+  - Julia:   `FreeBoundaryStability/eigenmode_energies[1]` in `gpec.h5`
   - Fortran: the "Energies: ... real = …, imaginary = …" line printed by `dcon`
              (== `W_t_eigenvalue[1]`; parsing the log avoids a NetCDF dependency).
 
@@ -157,7 +157,7 @@ function run_julia_reference()
     GPE.main([rundir])
     wall = time() - t0
     et = h5open(joinpath(rundir, "gpec.h5"), "r") do h5
-        read(h5["FreeBoundaryStability/XiNorm/eigenmode_energies"])
+        read(h5["FreeBoundaryStability/eigenmode_energies"])
     end
     return real(et[1]), imag(et[1]), wall
 end
diff --git a/benchmarks/compare_gal_vs_el.jl b/benchmarks/compare_gal_vs_el.jl
index f3facb53..cfb1d527 100644
--- a/benchmarks/compare_gal_vs_el.jl
+++ b/benchmarks/compare_gal_vs_el.jl
@@ -3,7 +3,7 @@
 #
 # EL :  ξ_EL(ψ)  = U_EL(ψ) · (U_EL_edge \ w)          (fundamental matrix, integration/xi_psi)
 # gal:  ξ_gal(ψ) = U_gal(ψ) · w                       (identity-at-edge ⇒ coefficient is w itself)
-# w = eigenvector of the total energy operator W = W_plasma + W_vacuum (FreeBoundaryStability/XiNorm).
+# w = eigenvector of the total energy operator W = W_plasma + W_vacuum (FreeBoundaryStability).
 #
 # Needs ONE gpec.h5 from a run with populate_dense_xi=true (EL dense u_store), gal_match_flag=true,
 # gal_ideal_flag=true (gal ideal matched), vac_flag=true (energy operator).
@@ -19,8 +19,8 @@ ksel = length(ARGS) >= 3 ? ARGS[3] : "highest"
 to_c(a) = eltype(a) <: Complex ? ComplexF64.(a) : map(x -> ComplexF64(x.re, x.im), a)
 
 et, wt, u1, psiE, gxi, psiG, issing, mlow, sing_psi = h5open(h5path) do f
-    (to_c(read(f["FreeBoundaryStability/XiNorm/eigenmode_energies"])),
-        to_c(read(f["FreeBoundaryStability/XiNorm/W_freeboundary_eigenmodes"])),
+    (to_c(read(f["FreeBoundaryStability/eigenmode_energies"])),
+        to_c(read(f["FreeBoundaryStability/W_freeboundary_eigenmodes"])),
         to_c(read(f["integration/xi_psi"])), read(f["integration/psi"]),
         to_c(read(f["galerkin/match/xi"])), read(f["galerkin/solution/psi"]),
         Bool.(read(f["galerkin/solution/issing"])), read(f["info/mlow"]),
diff --git a/benchmarks/plot_xi_eigenmode.jl b/benchmarks/plot_xi_eigenmode.jl
index 842e0cc5..e7c67520 100644
--- a/benchmarks/plot_xi_eigenmode.jl
+++ b/benchmarks/plot_xi_eigenmode.jl
@@ -17,8 +17,8 @@ outpng = length(ARGS) >= 2 ? ARGS[2] : joinpath(@__DIR__, "xi_eigenmode.png")
 to_c(a) = eltype(a) <: Complex ? ComplexF64.(a) : map(x -> ComplexF64(x.re, x.im), a)
 
 et, wt, u1, psi, mlow, sing_psi = h5open(h5path) do f
-    (to_c(read(f["FreeBoundaryStability/XiNorm/eigenmode_energies"])),
-        to_c(read(f["FreeBoundaryStability/XiNorm/W_freeboundary_eigenmodes"])),
+    (to_c(read(f["FreeBoundaryStability/eigenmode_energies"])),
+        to_c(read(f["FreeBoundaryStability/W_freeboundary_eigenmodes"])),
         to_c(read(f["integration/xi_psi"])),
         read(f["integration/psi"]),
         read(f["info/mlow"]),
diff --git a/docs/src/conventions.md b/docs/src/conventions.md
index 268fde01..c9f6838c 100644
--- a/docs/src/conventions.md
+++ b/docs/src/conventions.md
@@ -317,9 +317,16 @@ The internal flux-conform operator is just ``R = \Sigma\sqrt{A} = `` `rootarea_t
   `reluctance`, `plasma_inductance`, `surface_inductance`) are stored in this coordinate-invariant space
   under `perturbed_equilibrium/response_matrices/`. The stored `rootarea_to_area_weight_operator` ``S``
   recovers the area-weighted field forms (``L_{\bar b} = S\,\tilde L\,S^\dagger``) and the scalar
-  `surface_area` ``A`` recovers flux (``\Phi = A\,\bar b``). The coordinate-invariant ideal-MHD energies
-  written by the stability stage (`FreeBoundaryStability/eigenmode_energies`) are the ``\tilde b``
-  quadratic form scaled by the scalar ``c = A``: ``\mathrm{d}W = c\,\tilde b^\dagger W_t\,\tilde b``.
+  `surface_area` ``A`` recovers flux (``\Phi = A\,\bar b``).
+- **Eigenmode energies** (`FreeBoundaryStability/eigenmode_energies`) — the generalized eigenvalues of the
+  pencil ``W\,v = \lambda\,N\,v``, where ``N`` is the power-normalization (surface-norm) matrix built from
+  the Jacobian Fourier coefficients: ``\xi^\dagger N\,\xi = \oint J\,|\xi(\theta)|^2\,d\theta\,/\,V' =
+  \langle|\xi|^2\rangle``. The eigenvalues are stationary values of the power quotient
+  ``\xi^\dagger W \xi / \xi^\dagger N \xi`` — energies per unit surface-averaged squared displacement
+  (a normalization convention, not raw joules). Because ``W`` and ``N`` transform by the same congruence
+  under a change of working (Jacobian) coordinate, the spectrum is invariant to that coordinate choice.
+  ``N`` is not an inertia matrix (the stability stage carries no velocities); it is the DCON power
+  normalization promoted from a post-hoc per-eigenvector division into the metric of the eigenproblem.
 
 
 ## Rotation Velocity Conventions (KineticForces)
diff --git a/examples/Solovev_kinetic_calculated_example/gpec.toml b/examples/Solovev_kinetic_calculated_example/gpec.toml
index ea3913fe..5614dbca 100644
--- a/examples/Solovev_kinetic_calculated_example/gpec.toml
+++ b/examples/Solovev_kinetic_calculated_example/gpec.toml
@@ -1,6 +1,5 @@
 [Equilibrium]
 eq_type = "sol"                         # Type of the input 2D equilibrium file
-eq_filename = "sol.toml"                # Path to equilibrium file
 jac_type = "pest"                       # Coordinate system (hamada, pest, boozer, equal_arc)
 power_bp = 0                            # Poloidal field power exponent for Jacobian
 power_b = 0                             # Toroidal field power exponent for Jacobian
@@ -48,3 +47,15 @@ singfac_min = 1e-4            # Fractional distance from rational q at which ide
 ucrit = 1e3                   # Maximum fraction of solutions allowed before re-normalized
 write_outputs_to_HDF5 = true
 verbose = false
+
+[SOL_INPUT]
+mr = 128       # Number of radial grid zones
+mz = 128       # Number of axial grid zones
+ma = 128       # Number of flux grid zones
+e = 1.6        # Elongation
+a = 0.33       # Minor radius
+r0 = 1.0       # Major radius
+q0 = 1.9       # Safety factor at the magnetic axis (O-point)
+p0fac = 1      # Scale on-axis pressure (β changes; Φ, q fixed)
+b0fac = 1      # Scale toroidal field at constant β (Bt changes; shape, β fixed)
+f0fac = 1      # Scale toroidal field at constant pressure (β, q change; Φ, p, Bp fixed)
diff --git a/examples/Solovev_kinetic_calculated_example/sol.toml b/examples/Solovev_kinetic_calculated_example/sol.toml
deleted file mode 100644
index e54fdeac..00000000
--- a/examples/Solovev_kinetic_calculated_example/sol.toml
+++ /dev/null
@@ -1,11 +0,0 @@
-[SOL_INPUT]
-mr = 128      # number of radial grid zones
-mz = 128      # number of axial grid zones
-ma = 128      # number of flux grid zones
-e = 1.6       # elongation
-a = 0.33      # minor radius
-r0 = 1.0      # major radius
-q0 = 1.9      # safety factor at the o-point
-p0fac=1       # scale on-axis pressure (P-> P+P0*p0fac. beta changes. Phi,q constant)
-b0fac=1       # scale toroidal field at constant beta (s*Phi,s*f,s^2*P. bt changes. Shape,beta constant)
-f0fac=1       # scale toroidal field at constant pressure (s*f. beta,q changes. Phi,p,bp constant)
diff --git a/regression-harness/cases/diiid_n1.toml b/regression-harness/cases/diiid_n1.toml
index 1a69b727..b8090835 100644
--- a/regression-harness/cases/diiid_n1.toml
+++ b/regression-harness/cases/diiid_n1.toml
@@ -8,12 +8,12 @@ name = "diiid_n1"
 description = "DIII-D-like equilibrium, n=1, ideal + perturbed equilibrium"
 example_dir = "examples/DIIID-like_ideal_example"
 
-# Energies — leading eigenvalues at the final truncation (psilim). eigenmode_* are now root-area-weighted (Φ-space, Jacobian-invariant); ξ-space counterparts live under FreeBoundaryStability/XiNorm/ and are not tracked.
+# Energies — leading generalized (W,N) pencil eigenvalues at the final truncation (psilim): power-normalized (⟨|ξ|²⟩ = 1 metric) and invariant to the working-coordinate Jacobian.
 [quantities.et_real]
 h5path = "FreeBoundaryStability/eigenmode_energies"
 type = "complex_vector"
 extract = "real_first"
-label = "root-area-weighted total energy Re(et[1])"
+label = "total energy Re(et[1])"
 noise_threshold = 1e-10
 order = 10
 
@@ -21,7 +21,7 @@ order = 10
 h5path = "FreeBoundaryStability/eigenmode_energies"
 type = "complex_vector"
 extract = "imag_first"
-label = "root-area-weighted total energy Im(et[1])"
+label = "total energy Im(et[1])"
 noise_threshold = 1e-10
 order = 11
 
@@ -29,7 +29,7 @@ order = 11
 h5path = "FreeBoundaryStability/eigenmode_plasma_energies"
 type = "complex_vector"
 extract = "real_first"
-label = "root-area-weighted plasma energy Re(ep[1])"
+label = "plasma energy Re(ep[1])"
 noise_threshold = 1e-10
 order = 12
 
@@ -37,15 +37,15 @@ order = 12
 h5path = "FreeBoundaryStability/eigenmode_vacuum_energies"
 type = "complex_vector"
 extract = "real_first"
-label = "root-area-weighted vacuum energy Re(ev[1])"
+label = "vacuum energy Re(ev[1])"
 noise_threshold = 1e-10
 order = 13
 
 [quantities.vacuum_eigenvalue]
-h5path = "FreeBoundaryStability/XiNorm/vacuum_eigenvalue"
+h5path = "FreeBoundaryStability/vacuum_eigenvalue"
 type = "real_scalar"
 extract = "value"
-label = "vacuum matrix min eigenvalue (XiNorm)"
+label = "vacuum matrix min eigenvalue"
 noise_threshold = 1e-10
 order = 14
 
@@ -54,7 +54,7 @@ order = 14
 h5path = "FreeBoundaryStability/eigenmode_plasma_energies"
 type = "complex_vector"
 extract = "all_complex"
-label = "root-area-weighted plasma energy (all)"
+label = "plasma energy (all)"
 noise_threshold = 1e-10
 order = 20
 
@@ -62,7 +62,7 @@ order = 20
 h5path = "FreeBoundaryStability/eigenmode_vacuum_energies"
 type = "complex_vector"
 extract = "all_complex"
-label = "root-area-weighted vacuum energy (all)"
+label = "vacuum energy (all)"
 noise_threshold = 1e-10
 order = 21
 
@@ -70,7 +70,7 @@ order = 21
 h5path = "FreeBoundaryStability/eigenmode_energies"
 type = "complex_vector"
 extract = "all_complex"
-label = "root-area-weighted total energy (all)"
+label = "total energy (all)"
 noise_threshold = 1e-10
 order = 22
 
diff --git a/regression-harness/cases/solovev_kinetic_calculated.toml b/regression-harness/cases/solovev_kinetic_calculated.toml
index aa59b4db..e5ebb29e 100644
--- a/regression-harness/cases/solovev_kinetic_calculated.toml
+++ b/regression-harness/cases/solovev_kinetic_calculated.toml
@@ -9,28 +9,28 @@ example_dir = "examples/Solovev_kinetic_calculated_example"
 
 # Energies — leading eigenvalues
 [quantities.et_real]
-h5path = "FreeBoundaryStability/XiNorm/eigenmode_energies"
+h5path = "FreeBoundaryStability/eigenmode_energies"
 type = "complex_vector"
 extract = "real_first"
 label = "total energy Re(et[1])"
 noise_threshold = 1e-10
 
 [quantities.et_imag]
-h5path = "FreeBoundaryStability/XiNorm/eigenmode_energies"
+h5path = "FreeBoundaryStability/eigenmode_energies"
 type = "complex_vector"
 extract = "imag_first"
 label = "total energy Im(et[1])"
 noise_threshold = 1e-10
 
 [quantities.ep_real]
-h5path = "FreeBoundaryStability/XiNorm/eigenmode_plasma_energies"
+h5path = "FreeBoundaryStability/eigenmode_plasma_energies"
 type = "complex_vector"
 extract = "real_first"
 label = "plasma energy Re(ep[1])"
 noise_threshold = 1e-10
 
 [quantities.ev_real]
-h5path = "FreeBoundaryStability/XiNorm/eigenmode_vacuum_energies"
+h5path = "FreeBoundaryStability/eigenmode_vacuum_energies"
 type = "complex_vector"
 extract = "real_first"
 label = "vacuum energy Re(ev[1])"
@@ -38,7 +38,7 @@ noise_threshold = 1e-10
 
 # Energies — full eigenvalue arrays
 [quantities.et_all]
-h5path = "FreeBoundaryStability/XiNorm/eigenmode_energies"
+h5path = "FreeBoundaryStability/eigenmode_energies"
 type = "complex_vector"
 extract = "all_complex"
 label = "total energy (all)"
diff --git a/regression-harness/cases/solovev_kinetic_nuzero.toml b/regression-harness/cases/solovev_kinetic_nuzero.toml
index 89f9b5aa..e51efe87 100644
--- a/regression-harness/cases/solovev_kinetic_nuzero.toml
+++ b/regression-harness/cases/solovev_kinetic_nuzero.toml
@@ -9,28 +9,28 @@ example_dir = "examples/Solovev_kinetic_calculated_example"
 
 # Energies — leading eigenvalues
 [quantities.et_real]
-h5path = "FreeBoundaryStability/XiNorm/eigenmode_energies"
+h5path = "FreeBoundaryStability/eigenmode_energies"
 type = "complex_vector"
 extract = "real_first"
 label = "total energy Re(et[1])"
 noise_threshold = 1e-10
 
 [quantities.et_imag]
-h5path = "FreeBoundaryStability/XiNorm/eigenmode_energies"
+h5path = "FreeBoundaryStability/eigenmode_energies"
 type = "complex_vector"
 extract = "imag_first"
 label = "total energy Im(et[1])"
 noise_threshold = 1e-10
 
 [quantities.ep_real]
-h5path = "FreeBoundaryStability/XiNorm/eigenmode_plasma_energies"
+h5path = "FreeBoundaryStability/eigenmode_plasma_energies"
 type = "complex_vector"
 extract = "real_first"
 label = "plasma energy Re(ep[1])"
 noise_threshold = 1e-10
 
 [quantities.ev_real]
-h5path = "FreeBoundaryStability/XiNorm/eigenmode_vacuum_energies"
+h5path = "FreeBoundaryStability/eigenmode_vacuum_energies"
 type = "complex_vector"
 extract = "real_first"
 label = "vacuum energy Re(ev[1])"
@@ -38,7 +38,7 @@ noise_threshold = 1e-10
 
 # Energies — full eigenvalue arrays
 [quantities.et_all]
-h5path = "FreeBoundaryStability/XiNorm/eigenmode_energies"
+h5path = "FreeBoundaryStability/eigenmode_energies"
 type = "complex_vector"
 extract = "all_complex"
 label = "total energy (all)"
diff --git a/regression-harness/cases/solovev_multi_n.toml b/regression-harness/cases/solovev_multi_n.toml
index c67882ca..e8cbaa21 100644
--- a/regression-harness/cases/solovev_multi_n.toml
+++ b/regression-harness/cases/solovev_multi_n.toml
@@ -7,12 +7,12 @@ name = "solovev_multi_n"
 description = "Solovev analytical equilibrium, multi-n, ideal stability"
 example_dir = "examples/Solovev_ideal_example_multi_n"
 
-# Energies — leading eigenvalues at the final truncation (psilim). eigenmode_* are now root-area-weighted (Φ-space, Jacobian-invariant); ξ-space counterparts live under FreeBoundaryStability/XiNorm/ and are not tracked.
+# Energies — leading generalized (W,N) pencil eigenvalues at the final truncation (psilim): power-normalized (⟨|ξ|²⟩ = 1 metric) and invariant to the working-coordinate Jacobian.
 [quantities.et_real]
 h5path = "FreeBoundaryStability/eigenmode_energies"
 type = "complex_vector"
 extract = "real_first"
-label = "root-area-weighted total energy Re(et[1])"
+label = "total energy Re(et[1])"
 noise_threshold = 1e-10
 order = 10
 
@@ -20,15 +20,15 @@ order = 10
 h5path = "FreeBoundaryStability/eigenmode_energies"
 type = "complex_vector"
 extract = "imag_first"
-label = "root-area-weighted total energy Im(et[1])"
+label = "total energy Im(et[1])"
 noise_threshold = 1e-10
 order = 11
 
 [quantities.vacuum_eigenvalue]
-h5path = "FreeBoundaryStability/XiNorm/vacuum_eigenvalue"
+h5path = "FreeBoundaryStability/vacuum_eigenvalue"
 type = "real_scalar"
 extract = "value"
-label = "vacuum matrix min eigenvalue (XiNorm)"
+label = "vacuum matrix min eigenvalue"
 noise_threshold = 1e-10
 order = 14
 
@@ -37,7 +37,7 @@ order = 14
 h5path = "FreeBoundaryStability/eigenmode_plasma_energies"
 type = "complex_vector"
 extract = "all_complex"
-label = "root-area-weighted plasma energy (all)"
+label = "plasma energy (all)"
 noise_threshold = 1e-10
 order = 20
 
@@ -45,7 +45,7 @@ order = 20
 h5path = "FreeBoundaryStability/eigenmode_vacuum_energies"
 type = "complex_vector"
 extract = "all_complex"
-label = "root-area-weighted vacuum energy (all)"
+label = "vacuum energy (all)"
 noise_threshold = 1e-10
 order = 21
 
@@ -53,7 +53,7 @@ order = 21
 h5path = "FreeBoundaryStability/eigenmode_energies"
 type = "complex_vector"
 extract = "all_complex"
-label = "root-area-weighted total energy (all)"
+label = "total energy (all)"
 noise_threshold = 1e-10
 order = 22
 
diff --git a/regression-harness/cases/solovev_n1.toml b/regression-harness/cases/solovev_n1.toml
index eecbf8cc..d6c255ef 100644
--- a/regression-harness/cases/solovev_n1.toml
+++ b/regression-harness/cases/solovev_n1.toml
@@ -7,12 +7,12 @@ name = "solovev_n1"
 description = "Solovev analytical equilibrium, n=1, ideal stability"
 example_dir = "examples/Solovev_ideal_example"
 
-# Energies — leading eigenvalues at the final truncation (psilim). eigenmode_* are now root-area-weighted (Φ-space, Jacobian-invariant); ξ-space counterparts live under FreeBoundaryStability/XiNorm/ and are not tracked.
+# Energies — leading generalized (W,N) pencil eigenvalues at the final truncation (psilim): power-normalized (⟨|ξ|²⟩ = 1 metric) and invariant to the working-coordinate Jacobian.
 [quantities.et_real]
 h5path = "FreeBoundaryStability/eigenmode_energies"
 type = "complex_vector"
 extract = "real_first"
-label = "root-area-weighted total energy Re(et[1])"
+label = "total energy Re(et[1])"
 noise_threshold = 1e-10
 order = 10
 
@@ -20,7 +20,7 @@ order = 10
 h5path = "FreeBoundaryStability/eigenmode_energies"
 type = "complex_vector"
 extract = "imag_first"
-label = "root-area-weighted total energy Im(et[1])"
+label = "total energy Im(et[1])"
 noise_threshold = 1e-10
 order = 11
 
@@ -28,7 +28,7 @@ order = 11
 h5path = "FreeBoundaryStability/eigenmode_plasma_energies"
 type = "complex_vector"
 extract = "real_first"
-label = "root-area-weighted plasma energy Re(ep[1])"
+label = "plasma energy Re(ep[1])"
 noise_threshold = 1e-10
 order = 12
 
@@ -36,15 +36,15 @@ order = 12
 h5path = "FreeBoundaryStability/eigenmode_vacuum_energies"
 type = "complex_vector"
 extract = "real_first"
-label = "root-area-weighted vacuum energy Re(ev[1])"
+label = "vacuum energy Re(ev[1])"
 noise_threshold = 1e-10
 order = 13
 
 [quantities.vacuum_eigenvalue]
-h5path = "FreeBoundaryStability/XiNorm/vacuum_eigenvalue"
+h5path = "FreeBoundaryStability/vacuum_eigenvalue"
 type = "real_scalar"
 extract = "value"
-label = "vacuum matrix min eigenvalue (XiNorm)"
+label = "vacuum matrix min eigenvalue"
 noise_threshold = 1e-10
 order = 14
 
@@ -53,7 +53,7 @@ order = 14
 h5path = "FreeBoundaryStability/eigenmode_plasma_energies"
 type = "complex_vector"
 extract = "all_complex"
-label = "root-area-weighted plasma energy (all)"
+label = "plasma energy (all)"
 noise_threshold = 1e-10
 order = 20
 
@@ -61,7 +61,7 @@ order = 20
 h5path = "FreeBoundaryStability/eigenmode_vacuum_energies"
 type = "complex_vector"
 extract = "all_complex"
-label = "root-area-weighted vacuum energy (all)"
+label = "vacuum energy (all)"
 noise_threshold = 1e-10
 order = 21
 
@@ -69,7 +69,7 @@ order = 21
 h5path = "FreeBoundaryStability/eigenmode_energies"
 type = "complex_vector"
 extract = "all_complex"
-label = "root-area-weighted total energy (all)"
+label = "total energy (all)"
 noise_threshold = 1e-10
 order = 22
 
diff --git a/src/ForceFreeStates/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl
index daa37337..9f93e60a 100644
--- a/src/ForceFreeStates/EulerLagrange.jl
+++ b/src/ForceFreeStates/EulerLagrange.jl
@@ -939,9 +939,6 @@ function findmax_dW_edge!(odet::OdeState, ctrl::ForceFreeStatesControl, equil::E
     # Create a rough spline for wv matrix between psiedge -> psilim so we can approximate dW
     es.wvmat = free_compute_wv_spline(ctrl, equil, intr)
 
-    # Pre-compute sqrtamat + jarea spline for root-area-weighted eigenvalues
-    es.sqrtamat_spline = free_compute_sqrtamat_spline(ctrl, equil, intr)
-
     # Loop with compact index j into EdgeScanState; ODE index is edge_start + j - 1.
     # Steps where free_compute_total hits a singular wp solve are left as NaN per the EdgeScanState contract.
     for j in 1:N_edge
@@ -954,10 +951,6 @@ function findmax_dW_edge!(odet::OdeState, ctrl::ForceFreeStatesControl, equil::E
             es.plasma_energy[j] = result.plasma_energy
             es.vacuum_energy[j] = result.vacuum_energy
             es.vacuum_eigenvalue[j] = result.vacuum_eigenvalue
-            es.rootA_total_eigenvalue[j] = result.rootA_total_eigenvalue
-            es.rootA_plasma_energy[j] = result.rootA_plasma_energy
-            es.rootA_vacuum_energy[j] = result.rootA_vacuum_energy
-            es.rootA_vacuum_eigenvalue[j] = result.rootA_vacuum_eigenvalue
         catch e
             e isa LinearAlgebra.SingularException || rethrow()
         end
diff --git a/src/ForceFreeStates/ForceFreeStates.jl b/src/ForceFreeStates/ForceFreeStates.jl
index 5705c1a4..65d31ab6 100644
--- a/src/ForceFreeStates/ForceFreeStates.jl
+++ b/src/ForceFreeStates/ForceFreeStates.jl
@@ -33,7 +33,6 @@ include("FixedKineticMatrices.jl")
 include("Kinetic.jl")
 include("FixedBoundaryStability.jl")
 include("Utils.jl")
-include("RootAreaWeighted.jl")
 include("Free.jl")
 include("Riccati.jl")
 
diff --git a/src/ForceFreeStates/ForceFreeStatesStructs.jl b/src/ForceFreeStates/ForceFreeStatesStructs.jl
index c8702e00..6c7bccc3 100644
--- a/src/ForceFreeStates/ForceFreeStatesStructs.jl
+++ b/src/ForceFreeStates/ForceFreeStatesStructs.jl
@@ -371,8 +371,9 @@ end
     # Shared hint for sequential evaluation (all splines evaluated at same psi)
     _hint::Base.RefValue{Int} = Ref(1)
 
-    # Used in Free.jl
-    jmat::Vector{ComplexF64} = Vector{ComplexF64}(undef, 2 * mpert - 1)
+    # Jacobian Fourier band ψ-spline (2·mpert−1 conjugate-symmetric coefficients per surface),
+    # used to assemble the power-normalization matrix N in Free.jl
+    jmats::S = _empty_series_interp_complex(2 * mpert - 1, itp_opts)
 end
 
 # Helper to create empty complex series interpolant for default initialization
@@ -402,17 +403,15 @@ Populated in `Free.jl`.
   - `numpoints::Int` - Total number of points in the vacuum calculation (mthvac * nzvac)
   - `numpert_total::Int` - Total number of modes (mpert × npert)
   - `mthvac::Int` - Number of vacuum poloidal grid points (corresponds to `mtheta` in VacuumInput) - only needed for GPEC functionality currently
-  - `wt::Array{ComplexF64, 2}` - Free-boundary eigenvector matrix after diagonalising W (numpert_total × numpert_total). Columns are the eigenmodes sorted most-unstable first. **ξ-space.**
-  - `wt0::Array{ComplexF64, 2}` - Free-boundary total-energy matrix W = wp + wv before diagonalisation (numpert_total × numpert_total). **ξ-space.**
-  - `wp::Array{ComplexF64, 2}` - Plasma energy matrix (numpert_total × numpert_total). **ξ-space.**
-  - `wv::Array{ComplexF64, 2}` - Vacuum energy matrix (numpert_total × numpert_total). **ξ-space.**
-  - `rootA_wt0, rootA_wp, rootA_wv::Array{ComplexF64,2}` - Root-area-weighted flux (Φ-space) counterparts of `wt0`, `wp`, `wv` at the plasma edge (see RootAreaWeighted.jl). Jacobian-invariant up to the M†·W·M transform.
-  - `rootA_wt::Array{ComplexF64,2}` - Φ-space eigenvector matrix of `rootA_wt0` (columns sorted most-unstable first, phase-normalized so each column's largest-magnitude entry is real-positive).
-  - `ep::Vector{ComplexF64}` - Plasma eigenvalues
-  - `ev::Vector{ComplexF64}` - Vacuum eigenvalues
-  - `et::Vector{ComplexF64}` - Total eigenvalues of plasma + vacuum
+  - `wt::Array{ComplexF64, 2}` - Free-boundary eigenvector matrix of the generalized eigenproblem W·v = λ·N·v (numpert_total × numpert_total). Columns are the eigenmodes sorted most-unstable first, normalized to unit power norm v†·N·v = 1.
+  - `wt0::Array{ComplexF64, 2}` - Free-boundary total-energy matrix W = wp + wv before diagonalisation (numpert_total × numpert_total). ξ Fourier basis.
+  - `wp::Array{ComplexF64, 2}` - Plasma energy matrix (numpert_total × numpert_total). ξ Fourier basis.
+  - `wv::Array{ComplexF64, 2}` - Vacuum energy matrix (numpert_total × numpert_total). ξ Fourier basis.
+  - `ep::Vector{ComplexF64}` - Plasma energy per eigenmode (power quotient v†·wp·v with v†·N·v = 1)
+  - `ev::Vector{ComplexF64}` - Vacuum energy per eigenmode (power quotient v†·wv·v with v†·N·v = 1)
+  - `et::Vector{ComplexF64}` - Total energy eigenvalues of the pencil (W, N): power-normalized and invariant to the working (Jacobian) coordinate; et = ep + ev per mode
   - `n_tor_idx::Vector{Int}` -  0-based toroidal mode number index of each sorted eigenvalue (numpert_total). Needed in `write_imas`
-  - `vacuum_eigenvalue::Float64` - Least stable (minimum) eigenvalue of the vacuum matrix wv, clamped to zero
+  - `vacuum_eigenvalue::Float64` - Least stable (minimum) eigenvalue of the pencil (wv, N), clamped to zero
   - `grri::Array{Float64, 2}` - Interior Green's function matrices (2 * mthvac * nzvac × 2 * numpert_total)
   - `grre::Array{Float64, 2}` - Exterior Green's function matrices (2 * mthvac * nzvac × 2 * numpert_total)
   - `plasma_pts::Array{Float64, 3}` - Cartesian coordinates of plasma points, shape (mthvac * nzvac) × 3 for (x, y, z)
@@ -433,17 +432,6 @@ Populated in `Free.jl`.
     n_tor_idx::Vector{Int} = zeros(Int, numpert_total)
     vacuum_eigenvalue::Float64 = NaN
 
-    # Root-area-weighted flux eigenvalues (Jacobian-invariant)
-    rootA_et::Vector{ComplexF64} = fill(complex(NaN), numpert_total)
-    rootA_ep::Vector{ComplexF64} = fill(complex(NaN), numpert_total)
-    rootA_ev::Vector{ComplexF64} = fill(complex(NaN), numpert_total)
-
-    # Root-area-weighted flux matrices at the plasma edge (Jacobian-invariant).
-    rootA_wt0::Array{ComplexF64,2} = fill(complex(NaN), numpert_total, numpert_total)
-    rootA_wp::Array{ComplexF64,2} = fill(complex(NaN), numpert_total, numpert_total)
-    rootA_wv::Array{ComplexF64,2} = fill(complex(NaN), numpert_total, numpert_total)
-    rootA_wt::Array{ComplexF64,2} = fill(complex(NaN), numpert_total, numpert_total)
-
     grri::Array{Float64,2} = Array{Float64}(undef, 2 * numpoints, 2 * numpert_total)
     grre::Array{Float64,2} = Array{Float64}(undef, 2 * numpoints, 2 * numpert_total)
     plasma_pts::Array{Float64,2} = Array{Float64}(undef, numpoints, 3)
@@ -456,18 +444,16 @@ VacuumData(numpoints::Int, numpert_total::Int, mthvac::Int) = VacuumData(; numpo
 EdgeScanState
 
 Holds the state and results for the edge dW stability scan over ψ ∈ [psiedge, psilim].
-Initialized and populated by `findmax_dW_edge!`; results written to HDF5 under `EdgeScan/`
-(root-area-weighted values) and `EdgeScan/XiNorm/` (ξ-space values, for benchmarking).
+Initialized and populated by `findmax_dW_edge!`; results written to HDF5 under `EdgeScan/`.
+The energies are generalized (W, N) pencil values: power-normalized and invariant to the
+working (Jacobian) coordinate (see `power_norm_matrix!`).
 
 ## Fields
 
   - `wvmat` - Precomputed wv matrix spline (raw, no singfac); singfac applied analytically in `free_compute_total`.
   - `wv_hint::Base.RefValue{Int}` - Search hint for wvmat spline (different grid from equilibrium profiles).
-  - `sqrtamat_spline` - Precomputed √A convolution matrix + surface area spline for root-area-weighted eigenvalues.
-  - `sqrtamat_hint::Base.RefValue{Int}` - Search hint for sqrtamat spline.
   - `psi, q` - ψ and q values at each edge scan step.
-  - `total_eigenvalue, plasma_energy, vacuum_energy, vacuum_eigenvalue` - ξ-space energy components at each step (NaN for failed steps). Kept for Fortran benchmarking.
-  - `rootA_total_eigenvalue, rootA_plasma_energy, rootA_vacuum_energy, rootA_vacuum_eigenvalue` - Root-area-weighted flux (Φ-space) energies (Jacobian-invariant; NaN for failed/singular steps). These are the values written to `EdgeScan/` by default.
+  - `total_eigenvalue, plasma_energy, vacuum_energy, vacuum_eigenvalue` - Power-normalized energy components at each step (NaN for steps where the wp solve was singular). These drive the truncation choice and are written to `EdgeScan/`.
 """
 @kwdef mutable struct EdgeScanState
     numpert_total::Int
@@ -477,10 +463,6 @@ Initialized and populated by `findmax_dW_edge!`; results written to HDF5 under `
     wvmat::CubicSeriesInterpolant{Float64,ComplexF64} = _empty_series_interp_complex(numpert_total^2)
     wv_hint::Base.RefValue{Int} = Ref(1)
 
-    # Root-area-weighted sqrtamat + jarea spline (mpert^2 + 1 series: flattened sqrtamat then jarea)
-    sqrtamat_spline::CubicSeriesInterpolant{Float64,ComplexF64} = _empty_series_interp_complex(numpert_total^2 + 1)
-    sqrtamat_hint::Base.RefValue{Int} = Ref(1)
-
     # Scan results (written to HDF5 under EdgeScan/; NaN where free_compute_total raised SingularException)
     psi::Vector{Float64} = Vector{Float64}(undef, N_edge)
     q::Vector{Float64} = Vector{Float64}(undef, N_edge)
@@ -488,12 +470,6 @@ Initialized and populated by `findmax_dW_edge!`; results written to HDF5 under `
     plasma_energy::Vector{ComplexF64} = fill(complex(NaN), N_edge)
     vacuum_energy::Vector{ComplexF64} = fill(complex(NaN), N_edge)
     vacuum_eigenvalue::Vector{Float64} = fill(NaN, N_edge)
-
-    # Root-area-weighted flux eigenvalues (Jacobian-invariant)
-    rootA_total_eigenvalue::Vector{ComplexF64} = fill(complex(NaN), N_edge)
-    rootA_plasma_energy::Vector{ComplexF64} = fill(complex(NaN), N_edge)
-    rootA_vacuum_energy::Vector{ComplexF64} = fill(complex(NaN), N_edge)
-    rootA_vacuum_eigenvalue::Vector{Float64} = fill(NaN, N_edge)
 end
 
 EdgeScanState(numpert_total::Int, N_edge::Int) = EdgeScanState(; numpert_total, N_edge)
diff --git a/src/ForceFreeStates/Fourfit.jl b/src/ForceFreeStates/Fourfit.jl
index 2e9d7109..1ce7fde5 100644
--- a/src/ForceFreeStates/Fourfit.jl
+++ b/src/ForceFreeStates/Fourfit.jl
@@ -361,6 +361,7 @@ function make_matrix(equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStates
     g12 = Vector{ComplexF64}(undef, 2 * intr.mpert - 1)
     jmat = Vector{ComplexF64}(undef, 2 * intr.mpert - 1)
     jmat1 = Vector{ComplexF64}(undef, 2 * intr.mpert - 1)
+    jmats_flat = zeros(ComplexF64, mpsi, 2 * intr.mpert - 1)
     a_inv_dmat_temp = Matrix{ComplexF64}(undef, intr.numpert_total, intr.numpert_total)
     a_inv_cmat_temp = Matrix{ComplexF64}(undef, intr.numpert_total, intr.numpert_total)
 
@@ -409,6 +410,8 @@ function make_matrix(equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStates
             jmat[k+intr.mpert] = conj(jmat[intr.mpert-k])
             jmat1[k+intr.mpert] = conj(jmat1[intr.mpert-k])
         end
+        # Save the Jacobian Fourier band at every surface for the power-normalization spline
+        @views jmats_flat[ipsi, :] .= jmat
 
         # TODO: for 3D, would need an additional nlow:nhigh loop here for n/n' coupling
         for n in intr.nlow:intr.nhigh
@@ -508,8 +511,8 @@ function make_matrix(equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStates
     # TODO: set powers
     # Do we need this yet? Only called if power_flag = true
 
-    # This is used in free_run
-    ffit.jmat = jmat
+    # Jacobian Fourier band ψ-spline, used for the power normalization in Free.jl
+    ffit.jmats = cubic_interp(metric.xs, Series(jmats_flat); ffit.itp_opts...)
 
     return ffit
 end
diff --git a/src/ForceFreeStates/Free.jl b/src/ForceFreeStates/Free.jl
index a5184af4..12491107 100644
--- a/src/ForceFreeStates/Free.jl
+++ b/src/ForceFreeStates/Free.jl
@@ -1,3 +1,29 @@
+"""
+    power_norm_matrix!(Nmat, jmat, mpert, npert, dV_dpsi) -> Nmat
+
+Assemble the power-normalization (surface-norm) matrix N from the conjugate-symmetric Jacobian
+Fourier band `jmat` (length 2·mpert−1, evaluated from the `ffit.jmats` spline), such that
+
+    ξ†·N·ξ = ∮ J |ξ(θ)|² dθ / (dV/dψ) = ⟨|ξ|²⟩
+
+is the flux-surface average of the squared boundary displacement — the DCON power normalization
+as a quadratic form. N is Hermitian Toeplitz within each n-block (block-diagonal over n since the
+Jacobian is axisymmetric) and positive definite (J > 0). It is the metric of the generalized
+eigenproblem W·v = λ·N·v solved in `free_run!` and `free_compute_total`: because W and N
+transform by the same congruence under a change of working (Jacobian) coordinate, the pencil
+eigenvalues are power-normalized mode energies that are invariant to that coordinate choice.
+"""
+function power_norm_matrix!(Nmat::AbstractMatrix{ComplexF64}, jmat::AbstractVector{ComplexF64}, mpert::Int, npert::Int, dV_dpsi::Float64)
+    fill!(Nmat, 0.0 + 0.0im)
+    for ipert_n in 1:npert
+        off = (ipert_n - 1) * mpert
+        for ipert_m in 1:mpert, jpert_m in 1:mpert
+            Nmat[off+jpert_m, off+ipert_m] = jmat[jpert_m-ipert_m+mpert] / dV_dpsi
+        end
+    end
+    return Nmat
+end
+
 """
     free_run!(odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal) -> VacuumData
 
@@ -34,16 +60,27 @@ and data dumping.
     singfac = vec((mlow:mhigh) .- qlim .* (nlow:nhigh)')
     vac_data.wv .*= singfac .* singfac'
 
-    # Least stable eigenvalue of the vacuum matrix alone (should be PSD; clamp numerical noise to zero)
-    vac_data.vacuum_eigenvalue = max(0.0, minimum(real.(eigvals(Hermitian(vac_data.wv)))))
+    # Power-normalization matrix N at the plasma edge: ξ†·N·ξ = ⟨|ξ|²⟩ (see power_norm_matrix!).
+    # The Jacobian band is evaluated at psilim (same surface as W), not at the last grid surface.
+    Nmat = zeros!(pool, ComplexF64, numpert_total, numpert_total)
+    jmat_edge = zeros!(pool, ComplexF64, 2 * mpert - 1)
+    ffit.jmats(jmat_edge, psilim; hint=ffit._hint)
+    power_norm_matrix!(Nmat, jmat_edge, mpert, npert, dV_dpsi)
+
+    # Least stable eigenvalue of the vacuum matrix alone, power-normalized via the pencil
+    # (wv, N) so it shares the units of the mode energies (should be PSD; clamp noise to zero)
+    vac_data.vacuum_eigenvalue = max(0.0, minimum(real.(eigvals(Hermitian(vac_data.wv), Hermitian(Nmat)))))
 
     # Preserve wp in vac_data so it can be written to HDF5 as W_plasma
     vac_data.wp .= wp
 
-    # Compute complex energy eigenvalues and vectors
+    # Complex energy eigenvalues and vectors of the generalized eigenproblem W·v = λ·N·v.
+    # The eigenvalues are stationary values of the power quotient ξ†Wξ/ξ†Nξ — power-normalized
+    # mode energies invariant to the working (Jacobian) coordinate, since W and N transform by
+    # the same congruence under a poloidal coordinate change.
     vac_data.wt .= wp .+ vac_data.wv
     vac_data.wt0 .= vac_data.wt
-    Ev = eigen(vac_data.wt)
+    Ev = eigen(vac_data.wt, Nmat)
     vac_data.et .= Ev.values
     eindex = sortperm(real.(vac_data.et); rev=true)
 
@@ -59,17 +96,11 @@ and data dumping.
         vac_data.n_tor_idx[ipert] = (imax - 1) ÷ mpert
     end
 
-    # Normalize eigenfunction and energy.
+    # Normalize eigenvectors to unit power norm v†·N·v = 1. The generalized eigenvalues are
+    # already the power-normalized energies, so et is not rescaled here.
     for isol in 1:numpert_total
-        norm = 0.0 + 0.0im
-        for ipert_n in 1:npert, ipert_m in 1:mpert, jpert_m in 1:mpert
-            ipert = (ipert_n - 1) * mpert + ipert_m
-            jpert = (ipert_n - 1) * mpert + jpert_m
-            norm += ffit.jmat[jpert_m-ipert_m+mpert] * vac_data.wt[ipert, isol] * conj(vac_data.wt[jpert, isol])
-        end
-        norm /= dV_dpsi
-        vac_data.wt[:, isol] ./= sqrt(norm)
-        vac_data.et[isol] /= norm
+        v = @view vac_data.wt[:, isol]
+        v ./= sqrt(real(dot(v, Nmat, v)))
     end
 
     # Normalize phase
@@ -89,21 +120,6 @@ and data dumping.
     vac_data.ep .= diag(wpt)
     vac_data.ev .= diag(wvt)
 
-    # Eigenspectrum of W_Φ at psilim — Jacobian-invariant energy values; see RootAreaWeighted.jl.
-    # Computed directly here (no spline); spline is used by the edge scan.
-    mtheta_eq = length(equil.rzphi_ys)
-    ft_rootA = Utilities.FourierTransforms.FourierTransform(mtheta_eq, mpert, mlow)
-    sqrtamat_rootA = Equilibrium.compute_sqrtamat(equil, psilim, ft_rootA)
-    jarea_rootA = Equilibrium.flux_surface_area(equil, psilim, mtheta_eq)
-    rootA_result = compute_rootarea_eigenvalues(vac_data.wt0, wp, vac_data.wv, sqrtamat_rootA, jarea_rootA, equil, psilim, intr; all_eigenvalues=true)
-    vac_data.rootA_et .= rootA_result.rootA_et_all
-    vac_data.rootA_ep .= rootA_result.rootA_ep_all
-    vac_data.rootA_ev .= rootA_result.rootA_ev_all
-    vac_data.rootA_wt0 .= rootA_result.wt_rootA
-    vac_data.rootA_wp .= rootA_result.wp_rootA
-    vac_data.rootA_wv .= rootA_result.wv_rootA
-    vac_data.rootA_wt .= rootA_result.rootA_eigenvectors
-
     # Normalize eigenvectors based on scaled wt
     coeffs = odet.u[:, :, 1, end] \ (vac_data.wt .* (2π * equil.psio * 1e-3))
     @views for istep in 1:odet.step
@@ -190,6 +206,8 @@ wv matrix spline to `free_compute_wv_spline` and pass it in `odet.edge_scan.wvma
     eigenvalues = zeros!(pool, ComplexF64, Npert)
     wt = zeros!(pool, ComplexF64, Npert, Npert)
     wv = zeros!(pool, ComplexF64, Npert, Npert)
+    Nmat = zeros!(pool, ComplexF64, Npert, Npert)
+    jmat_local = zeros!(pool, ComplexF64, 2 * intr.mpert - 1)
     eindex = zeros!(pool, Int, Npert)
     evals_real = zeros!(pool, Float64, Npert)
     tmp_v = zeros!(pool, ComplexF64, Npert)
@@ -209,13 +227,15 @@ wv matrix spline to `free_compute_wv_spline` and pass it in `odet.edge_scan.wvma
     singfac = vec((intr.mlow:intr.mhigh) .- q_at_psifac .* (intr.nlow:intr.nhigh)')
     wv .*= singfac .* singfac'
 
-    # Compute total energy matrix and eigen-decomposition
-    wt .= wp .+ wv
+    # Local power-normalization matrix N(ψ) from the Jacobian Fourier band spline, so the
+    # power quotient uses the same surface as W (see power_norm_matrix!)
+    ffit.jmats(jmat_local, odet.psifac; hint=ffit._hint)
+    power_norm_matrix!(Nmat, jmat_local, intr.mpert, intr.npert, dV_dpsi)
 
-    # Save wt before eigen (which overwrites the input) for root-area-weighted computation
-    wt_saved = copy(wt)
-
-    Ev = eigen(wt)
+    # Total energy matrix and generalized eigen-decomposition of the pencil (W, N) — the
+    # eigenvalues are power-normalized, Jacobian-invariant mode energies (see free_run!)
+    wt .= wp .+ wv
+    Ev = eigen(wt, Nmat)
 
     # Sort eigenvalues by descending real part
     @inbounds for i in 1:Npert
@@ -227,43 +247,19 @@ wv matrix spline to `free_compute_wv_spline` and pass it in `odet.edge_scan.wvma
         eigenvalues[ipert] = Ev.values[eindex[Npert+1-ipert]]
     end
 
-    # Compute kinetic norm xi'*J(psi)*xi / dV_dpsi for the leading eigenvector.
-    # This normalizes total_eigenvalue, plasma_energy, vacuum_energy to be dimensionally consistent with free_run!.
-    isol = 1
-    v = @view wt[:, isol]
-    norm = 0.0 + 0.0im
-    for ipert_n in 1:intr.npert, ipert_m in 1:intr.mpert, jpert_m in 1:intr.mpert
-        ipert = (ipert_n - 1) * intr.mpert + ipert_m
-        jpert = (ipert_n - 1) * intr.mpert + jpert_m
-        norm += ffit.jmat[jpert_m-ipert_m+intr.mpert] * v[ipert] * conj(v[jpert])
-    end
-    norm /= dV_dpsi
-    eigenvalues[isol] /= norm
-
-    # Plasma and vacuum energy components for the leading eigenvector, normalized by the same norm.
-    # plasma_energy + vacuum_energy = total_eigenvalue by construction (wt = wp + wv; eigenvalue = v'*wt*v / norm).
+    # Plasma and vacuum energy components for the leading eigenvector via power quotients:
+    # plasma_energy + vacuum_energy = total_eigenvalue (wt = wp + wv; λ = v'*wt*v / v'*N*v).
+    v = @view wt[:, 1]
+    vnorm = real(dot(v, Nmat, v))
     mul!(tmp_v, wp, v)
-    plasma_energy = ComplexF64(dot(v, tmp_v)) / norm
+    plasma_energy = ComplexF64(dot(v, tmp_v)) / vnorm
     mul!(tmp_v, wv, v)
-    vacuum_energy = ComplexF64(dot(v, tmp_v)) / norm
-
-    # Smallest eigenvalue of the vacuum matrix alone, normalized by the same kinetic norm as the other energies
-    # so all four outputs are directly comparable. The singfac-scaled wv should be PSD by
-    # construction (congruence of PSD wv_raw), but numerical noise can make eigenvalues slightly
-    # negative. Clamp to zero to enforce the physical constraint.
-    vacuum_eigenvalue = real(max(0.0, minimum(real.(eigvals(Hermitian(wv))))) / norm)
-
-    # Eigenspectrum of W_Φ — Jacobian-invariant energy values; see RootAreaWeighted.jl.
-    # Evaluate sqrtamat + jarea from the pre-computed spline.
-    mpert = intr.mpert
-    sqrtamat_flat = Vector{ComplexF64}(undef, mpert^2 + 1)
-    es.sqrtamat_spline(sqrtamat_flat, odet.psifac; hint=es.sqrtamat_hint)
-    sqrtamat_local = reshape(@view(sqrtamat_flat[1:(mpert^2)]), mpert, mpert)
-    jarea_local = real(sqrtamat_flat[end])
-
-    rootA_result = compute_rootarea_eigenvalues(wt_saved, wp, wv, sqrtamat_local, jarea_local, equil, odet.psifac, intr)
-
-    return (total_eigenvalue=eigenvalues[1], plasma_energy=plasma_energy, vacuum_energy=vacuum_energy, vacuum_eigenvalue=vacuum_eigenvalue,
-        rootA_total_eigenvalue=rootA_result.rootA_total_eigenvalue, rootA_plasma_energy=rootA_result.rootA_plasma_energy, rootA_vacuum_energy=rootA_result.rootA_vacuum_energy,
-        rootA_vacuum_eigenvalue=rootA_result.rootA_vacuum_eigenvalue)
+    vacuum_energy = ComplexF64(dot(v, tmp_v)) / vnorm
+
+    # Smallest eigenvalue of the vacuum matrix alone via the pencil (wv, N), so all four outputs
+    # share the power-normalized units. The singfac-scaled wv should be PSD by construction
+    # (congruence of PSD wv_raw); clamp numerical noise to zero.
+    vacuum_eigenvalue = max(0.0, minimum(real.(eigvals(Hermitian(wv), Hermitian(Nmat)))))
+
+    return (total_eigenvalue=eigenvalues[1], plasma_energy=plasma_energy, vacuum_energy=vacuum_energy, vacuum_eigenvalue=vacuum_eigenvalue)
 end
diff --git a/src/ForceFreeStates/RootAreaWeighted.jl b/src/ForceFreeStates/RootAreaWeighted.jl
deleted file mode 100644
index bbd45153..00000000
--- a/src/ForceFreeStates/RootAreaWeighted.jl
+++ /dev/null
@@ -1,221 +0,0 @@
-"""
-Root-area-weighted field eigenvalue computation for the FFS edge scan.
-
-Transforms the energy matrix W from ξ-space into the coordinate-invariant root-area-weighted
-field (b̃) space and returns the eigenvalues of the scaled quadratic form `c·W_t` (energies).
-The **eigenvalues** (energies) are coordinate-invariant: energy is a physically meaningful
-scalar and does not depend on the choice of straight-field-line coordinate.
-
-    dW = ξ†·W·ξ = c·b̃†·(M†·W·M)·b̃ = c·b̃†·W_t·b̃
-
-Transformation chain (no poloidal flux Φ appears — it would only be the scalar product `Φ = A·b̄`):
-  - T = diag(i·2π·χ₁·singfac) converts ξ → flux-like amplitude (χ₁ = 2π·ψ₀, singfac = m − n·q).
-  - Σ = sqrtamat is the √weight (the bare → root-area-weighted field map).
-  - M = T⁻¹·Σ maps the root-area-weighted field b̃ → ξ.
-  - c = jarea = ∫ J|∇ψ| dθ is the scalar surface area carrying the energy's area normalization.
-
-What is **NOT** invariant across Jacobian choices (see
-`scripts/test_power_norm_invariance.jl` for numerical proof):
-  - ξ itself (W is Jacobian-dependent).
-  - The θ-space field reconstructed from the eigenvector.
-
-What **IS** invariant (verified numerically in the test script, within the
-area-integral precision floor ~1e-6):
-  - Flux-surface area A = ∫ J|∇ψ| dθ.
-  - The √weight operator identity ‖sqrtamat·b_fft‖² = N²·∫|b|²·J|∇ψ| dθ.
-  - The angle-map convmat (b_jac2 = convmat·b_jac1) to machine precision.
-  - The eigenspectrum of `c·W_t`.
-"""
-
-# The √weight building blocks `compute_sqrt_jac_delpsi`, `compute_sqrtamat`, and the field-translation
-# operators now live in `Equilibrium/CoordinateInvariant.jl` so the ForceFreeStates and
-# PerturbedEquilibrium modules share one coordinate-invariant definition.
-
-"""
-    compute_rootarea_eigenvalues(wt, wp, wv, sqrtamat, jarea, equil, psi, intr; all_eigenvalues=false)
-
-Transform W from ξ-space to the root-area-weighted field (b̃) space and eigendecompose.
-
-The transformation chain:
-  1. T = diag(i·2π·χ₁·singfac) converts ξ → flux-like amplitude (singfac = m − n·q).
-  2. Σ = sqrtamat is the √weight (bare → root-area-weighted field).
-  3. M = T⁻¹·Σ maps the root-area-weighted field b̃ → ξ.
-  4. Energy matrix = c·W_t = c·M†·W·M with the scalar c = jarea (the surface area).
-
-Only the eigenspectrum of `c·W_t` is physically invariant across Jacobian choices —
-it is the coordinate-independent energy. The eigenvectors in this basis are
-normalized under Julia's `eigen` (‖v‖₂ = 1) rather than the ∫|b|²dA = 1
-convention, so `v` should not be reused as a physical b̃-space mode shape — here
-it is only used to split λ = v†Wt v into ⟨v,Wp v⟩ + ⟨v,Wv v⟩, both of which are
-spectrally invariant because Wp + Wv = Wt.
-
-Returns the standard eigenvalues of `c·W_t` sorted ascending (most negative/unstable
-first), matching the ξ-space convention in `free_run!` and `free_compute_total`.
-Also returns `rootA_vacuum_eigenvalue = max(0, λ_min(Hermitian(c·wv_t)))` — the b̃-space
-counterpart of the ξ-space `vacuum_eigenvalue` diagnostic.
-Returns NaN when any singfac ≈ 0 (rational surface crossing makes T singular).
-"""
-function compute_rootarea_eigenvalues(
-    wt::AbstractMatrix{ComplexF64},
-    wp::AbstractMatrix{ComplexF64},
-    wv::AbstractMatrix{ComplexF64},
-    sqrtamat::AbstractMatrix{ComplexF64},
-    jarea::Float64,
-    equil::Equilibrium.PlasmaEquilibrium,
-    psi::Float64,
-    intr::ForceFreeStatesInternal;
-    all_eigenvalues::Bool=false
-)
-    Npert = intr.numpert_total
-    mpert = intr.mpert
-    npert = intr.npert
-    chi1 = 2π * equil.psio
-
-    # Compute singfac = m - n·q for each mode
-    q_at_psi = equil.profiles.q_spline(psi)
-    singfac = vec((intr.mlow:intr.mhigh) .- q_at_psi .* (intr.nlow:intr.nhigh)')
-
-    # Check for rational surface proximity — T is singular there
-    if any(abs.(singfac) .< 1e-6)
-        if all_eigenvalues
-            nan_vec = fill(complex(NaN), Npert)
-            nan_mat = fill(complex(NaN), Npert, Npert)
-            return (rootA_total_eigenvalue=complex(NaN), rootA_plasma_energy=complex(NaN), rootA_vacuum_energy=complex(NaN),
-                rootA_vacuum_eigenvalue=NaN, rootA_et_all=nan_vec, rootA_ep_all=nan_vec, rootA_ev_all=nan_vec,
-                wt_rootA=nan_mat, wp_rootA=nan_mat, wv_rootA=nan_mat, rootA_eigenvectors=nan_mat)
-        else
-            return (rootA_total_eigenvalue=complex(NaN), rootA_plasma_energy=complex(NaN), rootA_vacuum_energy=complex(NaN),
-                rootA_vacuum_eigenvalue=NaN)
-        end
-    end
-
-    # T = diag(i·2π·chi1·singfac): ξ → flux-like amplitude
-    T_diag_inv = Vector{ComplexF64}(undef, Npert)
-    for i in 1:Npert
-        T_diag_inv[i] = 1.0 / (im * 2π * chi1 * singfac[i])
-    end
-
-    # √weight operator Σ = sqrtamat maps the bare field → root-area-weighted field b̃.
-    # For multi-n, expand the mpert×mpert block to Npert×Npert block-diagonal.
-    sqrtamat_full = sqrtamat
-    if npert != 1
-        sqrtamat_full = zeros(ComplexF64, Npert, Npert)
-        for in in 1:npert
-            r = ((in - 1) * mpert + 1):(in * mpert)
-            sqrtamat_full[r, r] .= sqrtamat
-        end
-    end
-
-    # M = diag(T_inv) · Σ: row-scale Σ by T_inv (maps the root-area-weighted field b̃ → ξ)
-    M = similar(sqrtamat_full)
-    for i in 1:Npert
-        M[i, :] .= T_diag_inv[i] .* sqrtamat_full[i, :]
-    end
-
-    # Energy is the b̃ quadratic form scaled by the scalar surface area c = jarea:
-    #   dW = c·b̃†·W_t·b̃,   W_t = M†·W·M,   c = jarea  [Pharr 2026, Eq. for coordinate-invariant energy]
-    # Folding the scalar c here keeps the eigenspectrum (the coordinate-invariant energy) physical [J].
-    wt_rootA = jarea .* (M' * wt * M)
-    wp_rootA = jarea .* (M' * wp * M)
-    wv_rootA = jarea .* (M' * wv * M)
-
-    # Smallest eigenvalue of the vacuum matrix alone in b̃-space, clamped to zero.
-    # wv_rootA should be PSD by congruence of the PSD ξ-space wv; numerical noise
-    # can make the smallest eigenvalue slightly negative.
-    rootA_vacuum_eigenvalue = real(max(0.0, minimum(real.(eigvals(Hermitian(wv_rootA))))))
-
-    # Eigendecompose the total energy in root-area-weighted space. Only the eigenspectrum
-    # is Jacobian-invariant; the eigenvectors are coordinate-dependent.
-    Ev = eigen(wt_rootA)
-    eigenvalues = Ev.values
-    vectors = Ev.vectors
-
-    # Sort ascending by real part (most negative/most unstable first) — matches ξ-space convention in free_run!/free_compute_total
-    eindex = sortperm(real.(eigenvalues); rev=true)
-
-    # λ = v†Wt v = v†Wp v + v†Wv v ; the plasma/vacuum split is the projection
-    # onto the same eigenvector and is therefore also Jacobian-invariant.
-    if all_eigenvalues
-        rootA_et_all = Vector{ComplexF64}(undef, Npert)
-        rootA_ep_all = Vector{ComplexF64}(undef, Npert)
-        rootA_ev_all = Vector{ComplexF64}(undef, Npert)
-        rootA_eigenvectors = Matrix{ComplexF64}(undef, Npert, Npert)
-
-        for i in 1:Npert
-            v = vectors[:, eindex[Npert+1-i]]
-            rootA_eigenvectors[:, i] .= v
-            rootA_et_all[i] = eigenvalues[eindex[Npert+1-i]]
-            rootA_ep_all[i] = dot(v, wp_rootA * v)
-            rootA_ev_all[i] = dot(v, wv_rootA * v)
-        end
-
-        # Phase convention: rotate each column so its largest-magnitude entry is real-positive
-        # (matches the ξ-space convention in free_run!). Magnitudes are preserved since the
-        # eigenvectors come from Julia's `eigen` with ‖v‖₂ = 1, the natural b̃-space norm.
-        for isol in 1:Npert
-            imax = argmax(abs.(@view rootA_eigenvectors[:, isol]))
-            phase = abs(rootA_eigenvectors[imax, isol]) / rootA_eigenvectors[imax, isol]
-            @view(rootA_eigenvectors[:, isol]) .*= phase
-        end
-
-        return (rootA_total_eigenvalue=rootA_et_all[1], rootA_plasma_energy=rootA_ep_all[1], rootA_vacuum_energy=rootA_ev_all[1],
-            rootA_vacuum_eigenvalue=rootA_vacuum_eigenvalue,
-            rootA_et_all=rootA_et_all, rootA_ep_all=rootA_ep_all, rootA_ev_all=rootA_ev_all,
-            wt_rootA=wt_rootA, wp_rootA=wp_rootA, wv_rootA=wv_rootA, rootA_eigenvectors=rootA_eigenvectors)
-    else
-        idx = eindex[Npert]
-        v = vectors[:, idx]
-        rootA_total = eigenvalues[idx]
-        rootA_plasma = ComplexF64(dot(v, wp_rootA * v))
-        rootA_vacuum = ComplexF64(dot(v, wv_rootA * v))
-
-        return (rootA_total_eigenvalue=rootA_total, rootA_plasma_energy=rootA_plasma, rootA_vacuum_energy=rootA_vacuum,
-            rootA_vacuum_eigenvalue=rootA_vacuum_eigenvalue)
-    end
-end
-
-"""
-    free_compute_sqrtamat_spline(ctrl, equil, intr) -> CubicSeriesInterpolant
-
-Pre-compute sqrtamat and jarea over the edge scan ψ range and return a cubic series
-interpolant. Uses the same q-evenly-spaced grid as `free_compute_wv_spline`.
-
-The interpolant stores `mpert^2 + 1` complex series per grid point:
-  - First mpert^2 values: flattened sqrtamat matrix
-  - Last value: jarea (stored as complex with zero imaginary part)
-"""
-function free_compute_sqrtamat_spline(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal)
-    profiles = equil.profiles
-    mtheta_eq = length(equil.rzphi_ys)
-    mpert = intr.mpert
-
-    # Same grid as wv spline: evenly spaced in q
-    qedge = profiles.q_spline(ctrl.psiedge)
-    npsi = max(4, ceil(Int, (intr.qlim - qedge) * intr.nhigh * 4))
-    psi_array = zeros(Float64, npsi + 1)
-
-    n_series = mpert^2 + 1  # sqrtamat (flattened) + jarea
-    data_array = zeros(ComplexF64, npsi + 1, n_series)
-
-    # Create FourierTransform once (same for all psi since mtheta/mpert/mlow don't change)
-    ft = Utilities.FourierTransforms.FourierTransform(mtheta_eq, mpert, intr.mlow)
-
-    for i in 1:(npsi+1)
-        qi = qedge + (intr.qlim - qedge) * ((i - 1) / npsi)
-        psii = ctrl.psiedge + (intr.psilim - ctrl.psiedge) * ((i - 1) / npsi)
-        psi_array[i] = find_zero(
-            (psi -> profiles.q_spline(psi) - qi,
-                psi -> profiles.q_deriv(psi)),
-            psii, Roots.Newton()
-        )
-
-        sqrtamat = Equilibrium.compute_sqrtamat(equil, psi_array[i], ft)
-        jarea = Equilibrium.flux_surface_area(equil, psi_array[i], mtheta_eq)
-
-        # Flatten sqrtamat into first mpert^2 entries, jarea as last entry
-        data_array[i, 1:mpert^2] .= vec(sqrtamat)
-        data_array[i, end] = complex(jarea, 0.0)
-    end
-
-    return cubic_interp(psi_array, Series(data_array); extrap=ExtendExtrap())
-end
diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl
index 51bd422b..298bf537 100755
--- a/src/GeneralizedPerturbedEquilibrium.jl
+++ b/src/GeneralizedPerturbedEquilibrium.jl
@@ -724,21 +724,16 @@ function write_outputs_to_HDF5(
         out_h5["integration/crit"] = odet.crit_store
 
         # Write edge stability scan data (only present when psiedge < psilim).
-        # Root-area-weighted (Φ-space) energies are the default — they are Jacobian-
-        # invariant. The ξ-space values sit under EdgeScan/XiNorm/ and are retained for
-        # benchmarking against the Fortran GPEC lineage.
+        # Generalized (W, N) pencil energies — power-normalized, Jacobian-invariant; these are
+        # the values findmax_dW_edge! uses to choose the truncation point.
         if !isempty(odet.edge_scan.psi)
             es = odet.edge_scan
             out_h5["EdgeScan/psi"] = es.psi
             out_h5["EdgeScan/q"] = es.q
-            out_h5["EdgeScan/total_energy"] = es.rootA_total_eigenvalue
-            out_h5["EdgeScan/plasma_energy"] = es.rootA_plasma_energy
-            out_h5["EdgeScan/vacuum_energy"] = es.rootA_vacuum_energy
-            out_h5["EdgeScan/vacuum_eigenvalue"] = es.rootA_vacuum_eigenvalue
-            out_h5["EdgeScan/XiNorm/total_energy"] = es.total_eigenvalue
-            out_h5["EdgeScan/XiNorm/plasma_energy"] = es.plasma_energy
-            out_h5["EdgeScan/XiNorm/vacuum_energy"] = es.vacuum_energy
-            out_h5["EdgeScan/XiNorm/vacuum_eigenvalue"] = es.vacuum_eigenvalue
+            out_h5["EdgeScan/total_energy"] = es.total_eigenvalue
+            out_h5["EdgeScan/plasma_energy"] = es.plasma_energy
+            out_h5["EdgeScan/vacuum_energy"] = es.vacuum_energy
+            out_h5["EdgeScan/vacuum_eigenvalue"] = es.vacuum_eigenvalue
         end
 
         # Write singular surface data
@@ -782,27 +777,20 @@ function write_outputs_to_HDF5(
         out_h5["singular/kinetic/scan_cond"] = intr.kinsing_scan_cond
         out_h5["singular/kinetic/scan_threshold"] = intr.kinsing_scan_threshold
 
-        # Write free-boundary stability data. Root-area-weighted (Φ-space) is the
-        # default — Jacobian-invariant. ξ-space counterparts sit under
-        # FreeBoundaryStability/XiNorm/ for Fortran benchmarking.
-        # W_freeboundary_eigenmodes holds the eigenvector matrix of W_freeboundary with
-        # columns sorted most-unstable first; the same phase normalization is applied in
-        # both spaces (largest-magnitude entry made real-positive).
-        out_h5["FreeBoundaryStability/W_freeboundary"] = ctrl.vac_flag ? vac_data.rootA_wt0 : ComplexF64[]
-        out_h5["FreeBoundaryStability/W_plasma"] = ctrl.vac_flag ? vac_data.rootA_wp : ComplexF64[]
-        out_h5["FreeBoundaryStability/W_vacuum"] = ctrl.vac_flag ? vac_data.rootA_wv : ComplexF64[]
-        out_h5["FreeBoundaryStability/W_freeboundary_eigenmodes"] = ctrl.vac_flag ? vac_data.rootA_wt : ComplexF64[]
-        out_h5["FreeBoundaryStability/eigenmode_energies"] = ctrl.vac_flag ? vac_data.rootA_et : ComplexF64[]
-        out_h5["FreeBoundaryStability/eigenmode_plasma_energies"] = ctrl.vac_flag ? vac_data.rootA_ep : ComplexF64[]
-        out_h5["FreeBoundaryStability/eigenmode_vacuum_energies"] = ctrl.vac_flag ? vac_data.rootA_ev : ComplexF64[]
-        out_h5["FreeBoundaryStability/XiNorm/W_freeboundary"] = ctrl.vac_flag ? vac_data.wt0 : ComplexF64[]
-        out_h5["FreeBoundaryStability/XiNorm/W_plasma"] = ctrl.vac_flag ? vac_data.wp : ComplexF64[]
-        out_h5["FreeBoundaryStability/XiNorm/W_vacuum"] = ctrl.vac_flag ? vac_data.wv : ComplexF64[]
-        out_h5["FreeBoundaryStability/XiNorm/W_freeboundary_eigenmodes"] = ctrl.vac_flag ? vac_data.wt : ComplexF64[]
-        out_h5["FreeBoundaryStability/XiNorm/eigenmode_energies"] = ctrl.vac_flag ? vac_data.et : ComplexF64[]
-        out_h5["FreeBoundaryStability/XiNorm/eigenmode_plasma_energies"] = ctrl.vac_flag ? vac_data.ep : ComplexF64[]
-        out_h5["FreeBoundaryStability/XiNorm/eigenmode_vacuum_energies"] = ctrl.vac_flag ? vac_data.ev : ComplexF64[]
-        out_h5["FreeBoundaryStability/XiNorm/vacuum_eigenvalue"] = ctrl.vac_flag ? vac_data.vacuum_eigenvalue : NaN
+        # Write free-boundary stability data. The eigenmode energies are the generalized
+        # eigenvalues of the pencil (W, N) with N the power-normalization (surface-norm) matrix:
+        # power-normalized mode energies (⟨|ξ|²⟩ = 1 metric) that are invariant to the choice of
+        # working (Jacobian) coordinate. W_freeboundary_eigenmodes holds the generalized
+        # eigenvectors, columns sorted most-unstable first, normalized to unit power norm with
+        # the largest-magnitude entry made real-positive.
+        out_h5["FreeBoundaryStability/W_freeboundary"] = ctrl.vac_flag ? vac_data.wt0 : ComplexF64[]
+        out_h5["FreeBoundaryStability/W_plasma"] = ctrl.vac_flag ? vac_data.wp : ComplexF64[]
+        out_h5["FreeBoundaryStability/W_vacuum"] = ctrl.vac_flag ? vac_data.wv : ComplexF64[]
+        out_h5["FreeBoundaryStability/W_freeboundary_eigenmodes"] = ctrl.vac_flag ? vac_data.wt : ComplexF64[]
+        out_h5["FreeBoundaryStability/eigenmode_energies"] = ctrl.vac_flag ? vac_data.et : ComplexF64[]
+        out_h5["FreeBoundaryStability/eigenmode_plasma_energies"] = ctrl.vac_flag ? vac_data.ep : ComplexF64[]
+        out_h5["FreeBoundaryStability/eigenmode_vacuum_energies"] = ctrl.vac_flag ? vac_data.ev : ComplexF64[]
+        out_h5["FreeBoundaryStability/vacuum_eigenvalue"] = ctrl.vac_flag ? vac_data.vacuum_eigenvalue : NaN
 
         # Cartesian surface point clouds used downstream for visualisation and
         # perturbed-equilibrium plotting.
diff --git a/test/runtests_fullruns.jl b/test/runtests_fullruns.jl
index ac11df52..b3f8daa8 100644
--- a/test/runtests_fullruns.jl
+++ b/test/runtests_fullruns.jl
@@ -25,9 +25,7 @@ using HDF5
             # not numeric regression — the mode-converged physical value is pinned in the
             # regression harness (examples/Solovev_kinetic_calculated_example). Assert only
             # nerfed-grid-robust facts: finite and positive (Solovev is stable).
-            # ξ-space (XiNorm) eigenvalues are read here because root-area weighting collapses
-            # these stable modes to ~0, making them unusable as a regression anchor.
-            et = read(h5["FreeBoundaryStability/XiNorm/eigenmode_energies"])
+            et = read(h5["FreeBoundaryStability/eigenmode_energies"])
             @test isfinite(real(et[1]))
             @test real(et[1]) > 0
         end
@@ -45,9 +43,7 @@ using HDF5
             # harness on the mode-converged deck, not here — et[1] is a near-marginal,
             # ill-conditioned, FP-reassociation-sensitive quantity on this grid (sign not even
             # robust across platforms), so no value is pinned.
-            # ξ-space (XiNorm) eigenvalues are read here because root-area weighting collapses
-            # these modes to ~0, making them unusable as a regression anchor.
-            et = read(h5["FreeBoundaryStability/XiNorm/eigenmode_energies"])
+            et = read(h5["FreeBoundaryStability/eigenmode_energies"])
             @test isfinite(real(et[1]))
         end
         rm(joinpath(ex4, "gpec.h5"); force=true)
@@ -59,7 +55,7 @@ using HDF5
     @test begin
         GeneralizedPerturbedEquilibrium.main([ex5])
         h5open(joinpath(ex5, "gpec.h5"), "r") do h5
-            et = read(h5["FreeBoundaryStability/XiNorm/eigenmode_energies"])
+            et = read(h5["FreeBoundaryStability/eigenmode_energies"])
             # Smoke test (nerfed mpsi=16, delta_m=0 deck): exercises the full self-consistent
             # KF→FFS kinetic-MHD path end-to-end. NO numeric value is pinned here — the prior
             # imag(et[1]) ≈ -0.711 rtol=0.08 pin was platform-fragile (failed on macOS aarch64
@@ -82,7 +78,7 @@ using HDF5
     @test begin
         GeneralizedPerturbedEquilibrium.main([ex6])
         h5open(joinpath(ex6, "gpec.h5"), "r") do h5
-            et = read(h5["FreeBoundaryStability/XiNorm/eigenmode_energies"])
+            et = read(h5["FreeBoundaryStability/eigenmode_energies"])
             # Smoke test (nerfed mpsi=16, delta_m=0 deck): exercises the collisionless
             # (nutype="zero") real-x-space energy-integral path end-to-end — the #281 fix —
             # without faulting/NaN (the bug this guards against). The precise collisionless