matt-pharr · Nov 4, 2026 · Nov 6, 2026 · Nov 7, 2026
diff --git a/docs/development/architecture.md b/docs/development/architecture.md
@@ -210,7 +210,7 @@ The complete GPEC analysis pipeline:
 ### Stability
 - `SingType` - Singular surface data including:
   - Rational surface location (ψ, ρ, q = m/n, dq/dψ)
-  - Δ' (tearing stability parameter) — **stub**; the valid Δ' is `ForceFreeStatesResult.delta_prime.matrix`
+  - Δ' (tearing stability parameter) — **stub**; the valid Δ' is `ForceFreeStatesResult.delta_prime` (`DeltaPrimeData`)
   - Asymptotic solution bases at the inner-layer boundaries
 - `ForceFreeStatesResult` - Published product of a solve: mode space, metric/matrix fits, singular
   surfaces, and the per-integrator products (ξ solution and its basis, free-boundary energies,

diff --git a/docs/development/hdf5-conventions.md b/docs/development/hdf5-conventions.md
@@ -6,7 +6,7 @@ Conventions for the structure and naming of the `gpec.h5` output file. Every wri
 
 **The schema must be intuitive to a plasma physicist who is not a developer of this code.** The top level largely mirrors the TOML sections / major `src/` modules — which are themselves organized along physics lines — but **physics intuition wins whenever the two diverge**. Worked examples of that rule:
 
-- All per-rational-surface stability results consolidate under `SingularSurfaces/`, regardless of which algorithm produced them: the GGJ coefficients sit beside the Δ′/PEST-3 matching results, and the Riccati BVP and the Galerkin outer-region solve write the *same* `Delta_prime_matrix`/`Delta_prime_raw`/`Delta_coil`/`pest3_*` paths rather than each owning a producer-named subgroup.
+- All per-rational-surface stability results consolidate under `SingularSurfaces/`, regardless of which algorithm produced them: the GGJ coefficients sit beside the Δ′/PEST-3 matching results, and the Riccati BVP and the Galerkin outer-region solve write the *same* `Delta_prime_*`/`Delta_EF`/`pest3_*` paths rather than each owning a producer-named subgroup.
 - `EulerLagrangeMatrices` over a bare `Matrices` — group names must say what the data *is*, not which array it came from. Same reasoning renamed `matrices_<method>/` → `KineticMatrices/`.
 
 Physics-topic groups elevated to top level (rather than nested under their producer): `Info/`, `Input/`, `SingularSurfaces/`, `LocalStability/`, `SurfaceGeometries/`.
@@ -42,7 +42,7 @@ Top level (11 groups):
 | `Equilibrium/` | Scalars (`beta_N`, `q_axis`, `q_95`, `I_p`, …) plus `Profiles/` (1-D on `psi`: 2piF, mu0p, dVdpsi, q) and `Geometry/` (2-D on `psi`×`theta`: rcoords, offset, nu, jac) |
 | `ForceFreeStates/` | `Solutions/ForwardIntegration/` (u-solutions), `Solutions/GalerkinIntegration/` (closed ξ profiles in the shared layout, the gal surface list, `gal_basis_output`-gated `Basis/`), `write_el_matrices`-gated `EulerLagrangeMatrices/{Ideal,Kinetic}`, `FreeBoundaryStability/`, `EdgeScan/` |
 | `LocalStability/` | Mercier `D_I`, resistive interchange `D_R`, `ballooning_Delta_prime` on `psi`; the ballooning α boundary on `ballooning_psi` |
-| `SingularSurfaces/` | Per-rational-surface data: `rational_psi`/`rational_q`/`rational_m`/`rational_n`, GGJ coefficients, `Delta_prime_matrix`/`Delta_prime_raw`/`Delta_coil`/`pest3_A`/`pest3_B`/`pest3_Gamma` (Riccati or Galerkin alike), the inner-layer `Match/` results, `Kinetic/` |
+| `SingularSurfaces/` | Per-rational-surface data: `rational_psi`/`rational_q`/`rational_m`/`rational_n`, GGJ coefficients, `Delta_prime_raw_freeboundary`/`Delta_prime_raw_fixedboundary`/`Delta_EF`/`edge_response`, the free-boundary `Delta_prime_matrix`/`pest3_A`/`pest3_B`/`pest3_Gamma` and their `_fixedboundary` counterparts (Riccati or Galerkin alike), the inner-layer `Match/` results, `Kinetic/` |
 | `PerturbedEquilibrium/` | `ForcingModes/`, `Response/`, `ResponseMatrices/`, `SingularCoupling/`, `Energies/`, `TorqueResponse/` (kinetic forward solves: ψ-resolved `T_xe`, `T_source`), control-surface spectra |
 | `KineticForces/` | `<method>/` (torque/energy profiles, `dTdpsi_ell` resolved by bounce harmonic on the `ell` scale, `KineticMatrices/`); multi-ion runs add `PerSpecies/<species>/<method>/` with the same per-method layout, summing to the top-level total |
 | `ErrorFields/` | `CoilSensitivities/` (per-coil-set control-surface spectra and their rigid shift/tilt derivatives, `DominantMode/` full-window projection); `MonteCarlo/` (intrinsic and corrected `\|δ\|` histograms over the sampled tolerances, per batch and averaged); `Risk/` (threshold density, `P(lock\|δ)`, locking probabilities, `ToleranceScan/`); `NTV/` (correction-coil overlap and NTV torque couplings per kAt) |

diff --git a/docs/src/api.md b/docs/src/api.md
@@ -129,21 +129,21 @@ ffs = solve(eq, Riccati(; nchunks=40); nn=1, vac_flag=true)
 ffs = solve(eq, Galerkin(; nx=512); nn=1)                                        # finite-element outer-region Δ′
 ```
 
-A result carries `nothing` in the fields its integrator does not produce, and consumers warn and skip rather than erroring; `plasma_response` in particular needs a `Forward()` solve with `wall=nothing` (it errors on a wall or a Galerkin result and skips with a warning on a Riccati one), see [Stability](stability.md).
+A result carries `nothing` in the fields its integrator does not produce, and consumers warn and skip rather than erroring; `plasma_response` in particular needs a `Forward()` solve, or a `Galerkin(; error_field_response=true)` solve with `vac_flag=true`, with `wall=nothing` (it errors on a wall, any other Galerkin result or a matched one, and skips with a warning on a Riccati one), see [Stability](stability.md).
 A kinetic solve takes `Forward()` with the problem's `kinetic` keyword: `kinetic=KineticForcesControl(; ...)` evaluates the drift-kinetic matrices on the equilibrium's kinetic profiles, `kinetic=TestKineticMatrices(; scale)` uses synthetic matrices, `kinetic_factor` scales either, and `km = kinetic_matrices(prob)` evaluates the source once for reuse in later solves ([`KineticMatrices`](@ref GeneralizedPerturbedEquilibrium.ForceFreeStates.KineticMatrices)).
 
 ## Inner-layer matching
 
 A `MatchProblem` holds the outer Δ′ a solve published plus per-surface layer parameters, and the inner-layer model passed to `solve` computes the layer response at the prescribed rotation; the result is a new `ForceFreeStatesResult` with the closure changed from `:ideal` to `:matched` and the eigenfunctions replaced, so layer-parameter scans cost one cheap match per point:
 
 ```julia
-ffs = solve(eq, Galerkin(; rpec_flag=true, cut_solution=true); nn=1)
+ffs = solve(eq, Galerkin(; error_field_response=true, cut_solution=true); nn=1)
 matched = solve(MatchProblem(ffs; eta=1e-6, rho=1e-7, rotation=0.0), GGJModel())
 bpens = [solve(MatchProblem(ffs; eta=1e-6, rho=1e-7, rotation=f), GGJModel()).bpen for f in 0.0:50.0:500.0]
 ```
 
 A scalar applies to every surface and a vector gives one value per entry of `ffs.surfaces`, core to edge.
-The problem needs a Δ′ payload with coil-response columns, so it accepts Galerkin (`rpec_flag=true`) and Riccati (`vac_flag=true`) results; a Riccati-fed match fills `bpen` and the resonant data but keeps `solution === nothing`.
+The problem needs a Δ′ payload with error-field (EF) columns, so it accepts Galerkin and Riccati results with `error_field_response=true` (the default; Riccati also needs `vac_flag=true`); a Riccati-fed match fills `bpen`, the resonant data and the matched `wp` but keeps `solution === nothing`.
 `GGJModel()` is the closure-capable model; `SLAYERModel()` is slab-only and drives the tearing solve instead.
 
 ## Tearing stability

diff --git a/docs/src/galerkin.md b/docs/src/galerkin.md
@@ -23,24 +23,24 @@ integrator = "Galerkin"
 [Galerkin]
 nx = 256                # Hermite-cubic elements per interval between singular surfaces
 pfac = 0.03             # grid packing ratio toward the singular surfaces
-rpec_flag = true        # edge coil-response columns (Delta_coil), needed by [Match] and [Tearing]
+error_field_response = true        # error-field (EF) columns (Delta_EF): the ideal ξ solution and W_plasma, needed by [Match]
 ```
 
 `cut_solution = true` (default `false`) adds the cut solution that composite inner-region profiles need.
 
-It replaces the radial integration rather than supplementing it: the run computes its own vacuum response at the control surface (with `vac_flag`) but produces no free-boundary energies, no ODE trace and no `crit` scan, it does not support the kinetic solve, and its result cannot drive `[PerturbedEquilibrium]`, which needs the free-boundary δW of a `"Forward"` solve.
-A `[Match]` section matches the inner layer to its Δ′ and gives a driven ``\xi`` solution; the composite inner-region profiles need `cut_solution = true` (`Galerkin(; cut_solution=true)` followed by `solve(MatchProblem(ffs))` in a script).
+It replaces the radial integration rather than supplementing it: the run computes its own vacuum response at the control surface (with `vac_flag`), produces no ODE trace and no `crit` scan, and does not support the kinetic solve. With `error_field_response = true` and `vac_flag = true` it also produces the free-boundary energies of its ideal solution, so it can drive `[PerturbedEquilibrium]` like a `"Forward"` solve; a `[Match]`ed result cannot yet.
+With `error_field_response = true` the error-field (EF) columns, one per edge mode with unit edge displacement and no inner-layer drive, are solved with the edge held fixed; their small-solution coefficients are `Delta_EF`, they are the ideal (perfectly shielded) ``\xi`` solution, and the force they exert on the edge is the fixed-boundary plasma energy matrix ``W_P`` of Eq. 66 of [Glasser (2016)](citations.md#ForceFreeStates).
+A `[Match]` section matches the inner layer to its Δ′ and replaces them with the driven ``\xi`` solution, whose edge response ``W_P`` is complex: its anti-Hermitian part is the torque absorbed in the layers. The composite inner-region profiles need `cut_solution = true` (`Galerkin(; cut_solution=true)` followed by `solve(MatchProblem(ffs))` in a script).
 
 ## Outputs
 
 | Dataset | Meaning |
 |---|---|
-| `SingularSurfaces/Delta_prime_matrix` | the Δ′ matrix in the PEST-3 convention, on the same path as the Riccati result, so downstream stages never branch on which formalism ran |
-| `SingularSurfaces/Delta_prime_raw`, `Delta_coil` | the side-major matching matrix and the coil-response columns (`rpec_flag = true`) |
-| `SingularSurfaces/pest3_A`, `pest3_B`, `pest3_Gamma` | the PEST-3 interchange and mixed blocks ``A'``, ``B'``, ``\Gamma'`` |
+| `SingularSurfaces/Delta_prime_matrix` | the Δ′ matrix in the PEST-3 convention, on the same path as the Riccati result, so downstream stages never branch on which formalism ran; free-boundary with `vac_flag = true` (with `error_field_response = true` a second factorization re-solves the big-solution columns with the vacuum edge, so the EF columns never enter it; `free_boundary_delta_prime = false` skips it for matching-only scans, leaving no tearing Δ′) and fixed-boundary otherwise |
+| `SingularSurfaces/Delta_prime_raw_freeboundary`, `Delta_prime_raw_fixedboundary`, `Delta_EF`, `edge_response` | the side-major matching matrices with the vacuum and the fixed edge, the error-field (EF) response (`error_field_response = true`), and the edge momentum of each fixed-edge big-solution column |
+| `SingularSurfaces/pest3_A`, `pest3_B`, `pest3_Gamma` (and `_fixedboundary`) | the PEST-3 interchange and mixed blocks ``A'``, ``B'``, ``\Gamma'`` |
 | `ForceFreeStates/Solutions/GalerkinIntegration/` | the closed ``\xi`` profiles on the Galerkin grid, the surface list, the factorization used (`solver`) and, with `basis_output = true`, the `Basis/` group |
-
-Nothing is written under `ForceFreeStates/FreeBoundaryStability/`.
+| `ForceFreeStates/FreeBoundaryStability/` | `W_plasma`, the plasma energy matrix ``W_P`` of the closed solution (`error_field_response = true`), and with `vac_flag = true` the vacuum and free-boundary energies and eigenmodes as for a `"Forward"` solve |
 
 ## Numerical notes
 

diff --git a/docs/src/index.md b/docs/src/index.md
@@ -24,7 +24,7 @@ The [Workflow](workflow.md) page maps every stage to its physics, deck sections
 | `ForceFreeStates` — ideal (or drift-kinetic) MHD stability by the Newcomb criterion, with the `Vacuum` response of the `[Wall]` (a wall case is a separate deck from a response case) | `[ForceFreeStates]` | `ForceFreeStates/FreeBoundaryStability/`, `SingularSurfaces/Delta_prime_matrix` |
 | `LocalStability` — Mercier, resistive-interchange and ballooning criteria | `[LocalStability]` | `LocalStability/` |
 | `Match` — the resistive inner layer matched to the outer-region Δ′ | `[Match]` | `SingularSurfaces/Match/` |
-| `PerturbedEquilibrium` — the plasma response to the field of `[ForcingTerms]`: resonant fields, islands, overlap; needs a `"Forward"` solve without a wall | `[PerturbedEquilibrium]` | `PerturbedEquilibrium/SingularCoupling/` |
+| `PerturbedEquilibrium` — the plasma response to the field of `[ForcingTerms]`: resonant fields, islands, overlap; needs a `"Forward"` (or `error_field_response` `"Galerkin"`) solve without a wall | `[PerturbedEquilibrium]` | `PerturbedEquilibrium/SingularCoupling/` |
 | `KineticForces` — the neoclassical toroidal viscosity torque of that response | `[KineticForces]` | `KineticForces/` |
 | `ErrorFields` — coil tolerances, locking risk and correction limits | `[ErrorFields]` | `ErrorFields/Risk/` |
 | `Tearing` — tearing growth rates and critical resonant fields from an inner-layer model | `[Tearing]` | `Tearing/Roots/` |

diff --git a/docs/src/inner_layer.md b/docs/src/inner_layer.md
@@ -50,7 +50,7 @@ Penetration is the first resolved saddle-node bifurcation on the path from ``Q_0
 
 ## Running it
 
-From `examples/DIIID-like_SLAYER_example/gpec.toml`, after a `[ForceFreeStates]` solve with `integrator = "Riccati"` and `vac_flag = true` (or `"Galerkin"` with `rpec_flag = true`), which provides the Δ′ matrix, and with a kinetic file on the equilibrium:
+From `examples/DIIID-like_SLAYER_example/gpec.toml`, after a `[ForceFreeStates]` solve with `integrator = "Riccati"` or `"Galerkin"` and `vac_flag = true`, which provides the Δ′ matrix, and with a kinetic file on the equilibrium:
 
 ```toml
 [Tearing]

diff --git a/docs/src/perturbed_equilibrium.md b/docs/src/perturbed_equilibrium.md
@@ -92,7 +92,7 @@ The coupling matrices `C_*` map the applied ``\tilde b`` to each of these quanti
 ## Running it
 
 The response needs a force-free-states solve with `integrator = "Forward"` and no `[Wall]`.
-`"Forward"` is the only integrator that retains the dense ``\xi`` solution the operators are built from: a `"Riccati"` solve (the deck default) makes `plasma_response` warn that the response is skipped and return empty operators, so the failure surfaces only when a field is applied, and a `"Galerkin"` or `[Match]`ed solve, which publishes no free-boundary δW, is an error.
+The operators are built from a dense ``\xi`` solution with its free-boundary δW, which a `"Forward"` solve and a `"Galerkin"` solve with `error_field_response = true` and `vac_flag = true` both provide: a `"Riccati"` solve (the deck default) makes `plasma_response` warn that the response is skipped and return empty operators, so the failure surfaces only when a field is applied, while a `"Galerkin"` solve without those flags or a `[Match]`ed solve, whose singular coupling is not defined yet, is an error. The ψ-resolved torque response (`compute_torque_response`) needs a `"Forward"` kinetic solve.
 A solve with a wall is an error too, because the surface inductance is the free-space one and ``\Lambda`` from a wall run would mix two vacuum models; the wall case of a stability study is therefore a separate deck.
 
 From `examples/DIIID-like_ideal_example/gpec.toml`, after the `[ForceFreeStates]` solve and a `[ForcingTerms]` drive:

diff --git a/docs/src/stability.md b/docs/src/stability.md
@@ -59,7 +59,7 @@ qlow = 1.02              # inner edge of the integration domain
 qhigh = 1e3              # outer bound of qlim (see Integration domain)
 delta_mlow = 8           # extra poloidal modes below the resonant range
 delta_mhigh = 8          # extra poloidal modes above it
-integrator = "Forward"   # "Riccati" is the default; "Forward" (and no [Wall]) is required for [PerturbedEquilibrium]
+integrator = "Forward"   # "Riccati" is the default; "Forward" or an error_field_response "Galerkin" (and no [Wall]) is required for [PerturbedEquilibrium]
 truncate_at_rational_offset = true
 rational_offset_fraction = 0.2
 
@@ -81,7 +81,7 @@ ffs = solve(eq, Riccati(; eulerlagrange_tolerance=1e-10); nn=1, delta_mlow=8, de
 
 real(ffs.free_boundary.et[1])                 # least free-boundary eigenmode energy: negative is unstable
 for (s, sing) in enumerate(ffs.surfaces)
-    println("q = $(sing.m[1])/$(sing.n[1]) at ψ_N = $(sing.psifac): Δ′ = $(real(ffs.delta_prime.matrix[s, s]))")
+    println("q = $(sing.m[1])/$(sing.n[1]) at ψ_N = $(sing.psifac): Δ′ = $(real(ffs.delta_prime.freeboundary.Delta[s, s]))")
 end
 ```
 
@@ -103,7 +103,7 @@ The number of Julia threads is set at startup (`-t N` or `JULIA_NUM_THREADS`); t
 | `SingularSurfaces/rational_psi`, `rational_q`, `rational_m`, `rational_n`, `dqdpsi` | the rational surfaces of the domain, core to edge |
 | `SingularSurfaces/D_I`, `E`, `F`, `G`, `H`, `K`, `M` | the Mercier index and the Glasser–Greene–Johnson layer coefficients at each surface |
 | `SingularSurfaces/ca_left`, `ca_right` | the asymptotic (large, small) coefficients on either side of each surface, from ideal crossings |
-| `SingularSurfaces/Delta_prime_matrix` | the Δ′ matrix in the PEST-3 convention, from `"Riccati"` or `"Galerkin"`; `Delta_prime_raw` is the side-major (left, right) matrix it is folded from, `Delta_coil` the edge coil-response columns |
+| `SingularSurfaces/Delta_prime_matrix` | the free-boundary Δ′ matrix in the PEST-3 convention, from `"Riccati"` or `"Galerkin"`, with `pest3_A`, `pest3_B`, `pest3_Gamma`; `Delta_prime_raw_freeboundary` is the side-major (left, right) matrix they are folded from. The same set with a `_fixedboundary` suffix is the Δ′ with a conformal wall on the plasma, `Delta_EF` the error-field response, and `edge_response` the edge momentum of each fixed-boundary big-solution column (what matching needs for the matched `W_plasma`) |
 | `SingularSurfaces/Kinetic/` | the surfaces of a kinetic solve |
 
 The energy matrices are power-normalized: the eigenmode energies are the generalized eigenvalues of the pencil ``W v = \lambda N v``, with ``N`` the surface-norm matrix ``\xi^\dagger N \xi = \langle|\xi|^2\rangle`` built from the Jacobian Fourier coefficients, so each is an energy per unit surface-averaged squared displacement rather than joules, and the spectrum is invariant to the working coordinate because ``W`` and ``N`` transform by the same congruence.