matt-pharr · Oct 27, 2026 · Oct 23, 2026 · Oct 24, 2026 · Oct 24, 2026 · Oct 24, 2026
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
@@ -77,5 +77,5 @@ repos:
       - id: toml-no-deprecated-keys
         name: 'TOML conventions: no deprecated config keys'
         language: pygrep
-        entry: '^(mer_flag|force_wv_symmetry|ode_flag|cyl_flag|mat_flag|reform_eq_with_psilim|use_riccati|use_parallel|parallel_threads|populate_dense_xi|gal_flag|nstep|diagnose_ca|diagnose|ksing|power_bp|power_b|power_r|power_rc|use_galgrid|moment|eq_out|theta_out|xlmda_out|eqpsi_out|save_records|gpec_file|fnml_flag|ellip_flag|diag_flag|force_xialpha|fixed_boundary|output_eigenmodes|filter_modes|singular_point_method|local_stability_flag)\s*='
+        entry: '^(mer_flag|force_wv_symmetry|ode_flag|cyl_flag|mat_flag|reform_eq_with_psilim|use_riccati|use_parallel|parallel_threads|populate_dense_xi|gal_flag|nstep|diagnose_ca|diagnose|ksing|power_bp|power_b|power_r|power_rc|use_galgrid|moment|eq_out|theta_out|xlmda_out|eqpsi_out|save_records|gpec_file|fnml_flag|ellip_flag|diag_flag|force_xialpha|fixed_boundary|output_eigenmodes|filter_modes|singular_point_method|compute_response|local_stability_flag)\s*='
         files: ^(examples/.*\.toml|test/test_data/.*\.toml)$
diff --git a/benchmarks/benchmark_against_fortran_run.jl b/benchmarks/benchmark_against_fortran_run.jl
@@ -440,7 +440,6 @@ function write_gpec_toml_coil(
         end
         println(io)
         println(io, "[PerturbedEquilibrium]")
-        println(io, "compute_response          = true")
         println(io, "compute_singular_coupling = true")
         println(io, "verbose                   = true")
         println(io, "write_outputs_to_HDF5     = true")
@@ -492,7 +491,6 @@ function write_gpec_toml_file(
         println(io, "forcing_data_format = \"ascii\"")
         println(io)
         println(io, "[PerturbedEquilibrium]")
-        println(io, "compute_response          = true")
         println(io, "compute_singular_coupling = true")
         println(io, "verbose                   = true")
         println(io, "write_outputs_to_HDF5     = true")

diff --git a/benchmarks/benchmark_diiid_ideal_ntv_torque.jl b/benchmarks/benchmark_diiid_ideal_ntv_torque.jl
@@ -265,6 +265,7 @@ function run_benchmark(fortran_dir::String=default_fortran_dir())
     # Construct PerturbedEquilibriumState with Fortran Clebsch data.
     # Fortran gpec_xclebsch already writes ξ^α = xms/χ₁, exactly what the KF operator consumes — feed directly.
     pe_state = PE.PerturbedEquilibriumState(;
+        response=PE.PlasmaResponse(; ffs),
         psi_grid=psi_grid_f,
         xi_modes=(
             clebsch_psi=clebsch_psi,

diff --git a/docs/development/architecture.md b/docs/development/architecture.md
@@ -215,12 +215,12 @@ The complete GPEC analysis pipeline:
   is built
 
 ### Perturbed Equilibrium
-- `PerturbedEquilibriumControl` - User-facing TOML configuration parameters
-- `PerturbedEquilibriumInternal` - Internal state with mode arrays
-- `PerturbedEquilibriumState` - Results including:
+- `PlasmaResponse` - Source-independent operators of one FFS result (`plasma_response`): inductances,
+  permeability and reluctance, the singular coupling matrices [n_rational × numpert_total] and their
+  dominant modes; applied to any number of external fields
+- `PerturbedEquilibriumState` - The perturbed equilibrium driven by one field (`perturbed_equilibrium`), including:
   - Response fields (ξ_modes, b_modes) in mode space
-  - Singular coupling matrices [msing × numpert_total]
-  - Island diagnostics (half-widths, Chirikov parameters)
+  - Applied resonant vectors and island diagnostics (half-widths, Chirikov parameters)
 - `ForcingMode` - External forcing specification (m, n, amplitude, phase)
 
 ## Module Dependencies

diff --git a/docs/development/toml-conventions.md b/docs/development/toml-conventions.md
@@ -11,7 +11,7 @@ modify them. **Do not invent a new convention** — match what the existing exam
 2. **Inline annotation on every variable line.** Every `key = value` gets a trailing
    `# description`. Source the wording from the matching config struct's docstring
    (`EquilibriumConfig`, `WallShapeSettings`, `ForceFreeStatesControl`, `ForcingTermsControl`,
-   `PerturbedEquilibriumControl`, `SolovevConfig`, `TJAnalyticConfig`, `KineticForcesControl`)
+   the `plasma_response`/`perturbed_equilibrium` docstrings, `SolovevConfig`, `TJAnalyticConfig`, `KineticForcesControl`)
    and keep it to one terse line. **Use the same description for the same variable across all
    files** — `examples/DIIID-like_ideal_example/gpec.toml` is the canonical reference for the
    common sections. Exception: in `regression-harness/cases/*.toml` the repeated

diff --git a/docs/src/api.md b/docs/src/api.md
@@ -11,7 +11,8 @@ using GeneralizedPerturbedEquilibrium
 eq  = PlasmaEquilibrium("input.geqdsk"; jac_type="hamada")
 ffs = solve(eq, Riccati(); nn=1, delta_mlow=8, delta_mhigh=8, vac_flag=true)
 rmp = RMPField("coils.dat")
-pe  = perturbed_equilibrium(ffs, rmp)
+resp = plasma_response(ffs)
+pe   = perturbed_equilibrium(resp, rmp)
 ```
 
 `solve(eq, alg; kwargs...)` is sugar for the canonical problem form: a
@@ -25,7 +26,7 @@ prob = EulerLagrangeProblem(eq; nn=1, delta_mlow=8, delta_mhigh=8, vac_flag=true
 ffs  = solve(prob, Riccati())
 ```
 
-The four objects map one-to-one onto the pipeline stages:
+The objects map one-to-one onto the pipeline stages:
 
   - [`PlasmaEquilibrium`](@ref) reads and processes the equilibrium — the `[Equilibrium]`
     section. Analytic equilibria take their parameters from their config object, with the
@@ -42,8 +43,9 @@ The four objects map one-to-one onto the pipeline stages:
     a `ForceFreeStatesResult`, exactly the object the TOML driver publishes.
   - [`RMPField`](@ref) describes the external field — the `[ForcingTerms]` section — without
     reading anything from disk until it is applied.
-  - `perturbed_equilibrium` runs the plasma-response stage — the `[PerturbedEquilibrium]`
-    section.
+  - [`plasma_response`](@ref) builds the plasma-response and singular-coupling operators of a
+    `ForceFreeStatesResult` once, and `perturbed_equilibrium(resp, rmp)` drives them with one field —
+    together the `[PerturbedEquilibrium]` section.
 
 ## Combining forcing sources
 
@@ -57,13 +59,15 @@ phase-rotates a source.
 ```julia
 nominal    = RMPField("nominal_efc.dat")
 weld_field = RMPField("weld_fields.h5")
-pe = perturbed_equilibrium(ffs, 2.0 * nominal + 0.5 * weld_field)
+pe = perturbed_equilibrium(resp, 2.0 * nominal + 0.5 * weld_field)
 ```
 
 This is the substrate for error-field workflows: build per-unit sources independently
 (a coil set per ampere, a displacement field per millimeter), then assemble physical cases
 by weighting and summing — without recomputing anything per combination until the final
-`perturbed_equilibrium` call.
+`perturbed_equilibrium` call, and with one `plasma_response` serving every combination.
+`RMPField(modes)` wraps a control-surface spectrum computed elsewhere, already in the unit-norm
+`Φ_x` convention, as a leaf that combines like any other.
 
 The weights are linear-combination coefficients, not physical amplitudes: sources that are
 not scalar multiples of each other get their own description and the algebra. A coil set
@@ -154,6 +158,10 @@ GeneralizedPerturbedEquilibrium.solve
 GeneralizedPerturbedEquilibrium.perturbed_equilibrium
 ```
 
+[`plasma_response`](@ref GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.plasma_response), which builds the
+[`PlasmaResponse`](@ref GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.PlasmaResponse) that `perturbed_equilibrium`
+applies, is documented with the [perturbed equilibrium](perturbed_equilibrium.md) module.
+
 ## Integrator selectors
 
 ```@autodocs

diff --git a/docs/src/error_fields.md b/docs/src/error_fields.md
@@ -439,7 +439,7 @@ table.abs_delta_rim_per_mm         # the same for tilt, as rim displacement at t
 table.tilt_sensitivity_per_deg[1, :]               # ∂δ/∂θx per degree, per coil set
 
 # In memory, from the run's returned state (no I/O)
-rc  = PerturbedEquilibrium.ResonantCoupling(run.pe, run.ffs)
+rc  = PerturbedEquilibrium.ResonantCoupling(run.pe.response)
 dom = PerturbedEquilibrium.dominant_coupling(rc; psi_low=0.5)
 table = EF.sensitivity_table(run.coil_sensitivities, dom; mode=1)
 ```

diff --git a/docs/src/kinetic_forces.md b/docs/src/kinetic_forces.md
@@ -153,7 +153,7 @@ The torque is three nested adaptive Gauss–Kronrod integrals: over ψ (paneled
 
 QuadGK bounds the error of the whole integral, not of each step as the LSODE integrators of Fortran PENTRC did, so a level tighter than its margin only multiplies evaluations: the ψ quadrature cannot see λ-level noise below `rtol_psi/10`. The calculated kinetic matrix build (`kinetic_source = "calculated"`) has no ψ quadrature; it uses `rtol_pitch` when set and 1e-3 otherwise. Tightening `rtol_psi` tightens the whole ladder.
 
-## Regularization: why kinetic runs set `regularization_width = 0`
+## Regularization: off by default, and never needed by kinetic runs
 
 `[PerturbedEquilibrium] regularization_width` smooths the displacements before they drive the NTV
 integrand, multiplying ``\xi^{\psi\prime}`` and ``\xi^\alpha`` by
@@ -172,7 +172,9 @@ for smoothing.) Logan, *Electromagnetic Torque in Tokamaks with Toroidal Asymmet
 Princeton, 2015), Ch. 7 Eq. 7.46 states the same resonance shifting/splitting result in GPEC's own
 notation.
 
-GPEC therefore **forces `regularization_width = 0` whenever `kinetic_factor > 0`**, logging the override.
+GPEC therefore **defaults `regularization_width` to 0** for every result, as the Fortran `reg_flag = f`
+template does. Set it, typically to 0.05, only to compute the NTV torque of an ideal solve; a kinetic run
+should leave it at 0.
 Leaving it on suppresses a finite physical response and does so inconsistently — ``\xi^\psi``
 is never regularized, so damping the other two breaks their near-resonance cancellation in
 ``\delta B/B`` and leaves a spurious residue driving the NTV integrand.

diff --git a/docs/src/perturbed_equilibrium.md b/docs/src/perturbed_equilibrium.md
@@ -1,18 +1,21 @@
 # Perturbed Equilibrium
 
-The `PerturbedEquilibrium` module computes the plasma response to external magnetic perturbations.
+The `PerturbedEquilibrium` module computes the plasma response to external magnetic perturbations,
+in two stages: `plasma_response(ffs)` builds the source-independent response and singular-coupling
+operators of a force-free-states result once, and `perturbed_equilibrium(resp, rmp)` drives them with
+one external field. A scan over fields reuses one response.
 
 ## Types
 
 ```@docs
-GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.PerturbedEquilibriumControl
-GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.PerturbedEquilibriumInternal
+GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.PlasmaResponse
 GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.PerturbedEquilibriumState
 ```
 
 ## Functions
 
 ```@docs
+GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.plasma_response
 GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.compute_perturbed_equilibrium
 GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.write_outputs_to_HDF5
 ```
@@ -53,7 +56,7 @@ the same way from a finished run in memory or from its `gpec.h5`:
 ```julia
 using GeneralizedPerturbedEquilibrium.PerturbedEquilibrium
 
-rc = ResonantCoupling("gpec.h5")                 # post hoc; or ResonantCoupling(pe_state, ffs) in memory
+rc = ResonantCoupling("gpec.h5")                 # post hoc; or ResonantCoupling(resp) in memory
 dom = dominant_coupling(rc; psi_low=0.0, psi_high=0.9)   # SVD over the surfaces in the window
 
 b̃ = rootarea_field(rc, coil_modes)               # unit-norm forcing modes → root-area-weighted field
@@ -78,16 +81,17 @@ GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.CORE_PSI_HIGH
 GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.rootarea_field
 GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.coupling_overlap
 GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.check_mode_basis
-GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.compute_dominant_coupling!
+GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.build_dominant_coupling!
 ```
 
 ## Torque response matrices
 
 A kinetic forward solve (`integrator = "forward"`, `kinetic_factor > 0`) carries the
 drift-kinetic response inside its Euler-Lagrange energy: the complex plasma energy inside a
 flux surface, `δW(ψ) = ξ†·U₂(ψ)·U₁(ψ)⁻¹·ξ/(2μ₀)`, has an anti-Hermitian part that is the
-neoclassical toroidal viscosity torque (Logan & Park 2013, eq. 19). With
-`compute_torque_response = true` the response stage refers that quadratic form back to the
+neoclassical toroidal viscosity torque (Logan & Park 2013, eq. 19). On request,
+`plasma_response(ffs; compute_torque_response=true)` (deck: `compute_torque_response = true`
+in `[PerturbedEquilibrium]`) refers that quadratic form back to the
 applied root-area-weighted control-surface field b̃ — through the permeability, the
 boundary flux-to-displacement factor `χ₁(m − n·q_lim)·2πi` and the b̃ → Φ conform `R` — and
 stores the result at every stored node as the **torque response matrix** `T_xe(ψ)`, cumulative
@@ -105,8 +109,9 @@ Two conventions are worth stating explicitly:
     its eigenvectors are the applied spectra of extremal torque and its eigenvalues their torque
     per unit `‖b̃‖²/2`. The anti-Hermitian part carries the energy.
 
-The forcing-source form `T_source(ψ) = M†·T·M` is written alongside, with one column of `M` per
-forcing source of the run's `[ForcingTerms]`, each as built: every coil set at its deck currents,
+Each `perturbed_equilibrium(resp, rmp)` call contracts it with that field into a `TorqueResponse`:
+the applied torque `T_applied` and the forcing-source form `T_source(ψ) = M†·T·M`, with one column
+of `M` per forcing source, each as built: every coil set at its deck currents,
 or the forcing file's spectrum as a single source. Driving the sources at scale factors `s` gives
 `s†·T_source·s/2`.
 
@@ -124,10 +129,11 @@ T_h = Hermitian((tr.T_xe[:, :, end] + tr.T_xe[:, :, end]') / 2)
 
 ```@docs
 GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.TorqueResponse
-GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.torque_response_matrices
+GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.torque_response_matrix
 GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.torque_profile
 GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.coil_flux_spectra
-GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.compute_torque_response!
+GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.build_torque_response!
+GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.apply_torque_response!
 ```
 
 ## Plotting per-surface results against ψ or q

diff --git a/docs/src/workflow.md b/docs/src/workflow.md
@@ -113,14 +113,16 @@ The single `gpec.toml` file supplies user-selected options to every module. The
 
 **Module**: `PerturbedEquilibrium`
 
-**Purpose**: Compute the self-consistent plasma response to the external perturbation specified in Stage 4. Uses the force-free eigenfunctions and vacuum matrices to project the response into mode space, then evaluates singular coupling diagnostics at each rational surface.
+**Purpose**: Compute the self-consistent plasma response to the external perturbation specified in Stage 4, in two steps:
+- `plasma_response(ffs)` builds the source-independent operators once from the force-free eigenfunctions and vacuum matrices: the permeability, inductance and reluctance matrices on the control surface, and the singular coupling matrices at each rational surface (`PlasmaResponse`).
+- `perturbed_equilibrium(resp, rmp)` applies them to one forcing field: it projects the response into mode space and evaluates the singular coupling diagnostics for that field. One response serves any number of fields.
 
 **Configuration**: `[PerturbedEquilibrium]` section of `gpec.toml`
 
 **Inputs**:
-- ForceFreeStates results from Stage 3 (eigenfunctions, singular surface data)
-- Vacuum response matrices from Stage 2
-- `ForcingMode` array from Stage 4
+- ForceFreeStates results from Stage 3 (eigenfunctions, singular surface data), for `plasma_response`
+- Vacuum response matrices from Stage 2, for `plasma_response`
+- The forcing field from Stage 4, for `perturbed_equilibrium`
 
 **Outputs** (`PerturbedEquilibriumState`):
 - Mode-space normal displacement ξ(m, ψ) across the plasma

diff --git a/examples/DIIID-like_error_field_example/gpec.toml b/examples/DIIID-like_error_field_example/gpec.toml
@@ -206,7 +206,6 @@ height = -1.5780                        # Hoop height: centroid Z of the pack cr
 currents = [55000.0]                    # 1 kA per turn × 55 turns; sensitivities scale linearly with current [A]
 
 [PerturbedEquilibrium]
-compute_response = true                 # Compute plasma response to forcing
 compute_singular_coupling = true        # Compute singular layer coupling metrics
 verbose = false                         # Enable verbose logging
 write_outputs_to_HDF5 = true            # Write perturbed equilibrium outputs to HDF5

diff --git a/examples/DIIID-like_gal_resistive_pe_example/gpec.toml b/examples/DIIID-like_gal_resistive_pe_example/gpec.toml
@@ -1,9 +1,9 @@
 # DIII-D-like H-mode DRIVEN resistive example: gal Δ′ + inner-layer matching, with DIII-D
-# C-coil forcing (n=1 cosine phasing) staged for PerturbedEquilibrium. The response and the
-# singular coupling are switched off: they require the free-boundary δW, which the Galerkin
-# formalism does not yet produce, and PE refuses them on a Galerkin or matched result. The
-# gal-matched resistive ξ and penetrated field are carried on the result for when gal-side
-# δW lands. Same equilibrium and Galerkin settings as the sibling
+# C-coil forcing (n=1 cosine phasing) staged for PerturbedEquilibrium. The run stops at the
+# PerturbedEquilibrium stage for now: the plasma response needs the free-boundary δW, which the
+# Galerkin formalism does not yet produce, and a follow-up makes it run. The gal-matched
+# resistive ξ and penetrated field are carried on the result for when gal-side δW lands.
+# Same equilibrium and Galerkin settings as the sibling
 # DIIID-like_gal_resistive_example; differences start at [ForcingTerms].
 [Equilibrium]
 eq_filename = "TkMkr_D3Dlike_Hmode.geqdsk" # Path to equilibrium file
@@ -103,8 +103,7 @@ name = "c"                              # DIII-D C-coil (6 coils); resolves to c
 currents = [1000.0, 500.0, -500.0, -1000.0, -500.0, 500.0]  # n=1 cosine phasing, 1 kA peak [A]
 
 [PerturbedEquilibrium]
-compute_response = false                # Compute plasma response to forcing
-compute_singular_coupling = false       # Compute singular layer coupling metrics
+compute_singular_coupling = true        # Compute singular layer coupling metrics
 verbose = true                          # Enable verbose logging
 write_outputs_to_HDF5 = true            # Write perturbed equilibrium outputs to HDF5
-regularization_width = 0.05             # Regularization width for singular surfaces (0 = disabled)
+regularization_width = 0.0              # Regularization width for singular surfaces (0 = disabled)
diff --git a/examples/DIIID-like_ideal_example/gpec.toml b/examples/DIIID-like_ideal_example/gpec.toml
@@ -75,7 +75,6 @@ name = "c"                              # Coil set; resolves to coil_geometries/
 currents = [1000.0, 500.0, -500.0, -1000.0, -500.0, 500.0]  # n=1 cosine phasing, 1 kA peak [A]
 
 [PerturbedEquilibrium]
-compute_response = true                 # Compute plasma response to forcing
 compute_singular_coupling = true        # Compute singular layer coupling metrics
 verbose = true                          # Enable verbose logging
 write_outputs_to_HDF5 = true            # Write perturbed equilibrium outputs to HDF5

diff --git a/examples/DIIID-like_ideal_multi_n_example/gpec.toml b/examples/DIIID-like_ideal_multi_n_example/gpec.toml
@@ -77,8 +77,7 @@ name = "c"                              # Coil set; resolves to coil_geometries/
 currents = [1800.0, -50.0, -450.0, -800.0, -450.0, -50.0]  # Sum of n=1, 2, 3 cosine phasings with 1000, 500, 300 A peaks; n=3 lies outside nn_high [A]
 
 [PerturbedEquilibrium]
-compute_response = true                 # Compute plasma response to forcing
 compute_singular_coupling = true        # Compute singular layer coupling metrics
 verbose = false                         # Enable verbose logging
 write_outputs_to_HDF5 = true            # Write perturbed equilibrium outputs to HDF5
-regularization_width = 0.05             # Regularization width for singular surfaces (0 = disabled)
+regularization_width = 0.0              # Regularization width for singular surfaces (0 = disabled)