diff --git a/benchmarks/benchmark_threads.jl b/benchmarks/benchmark_threads.jl
index 4831f6988..048c64c66 100644
--- a/benchmarks/benchmark_threads.jl
+++ b/benchmarks/benchmark_threads.jl
@@ -32,7 +32,7 @@ function run_ffs(ex; use_parallel, use_riccati=false)
     metric = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_metric(equil, intr.mpert)
     ffit = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_matrix(equil, intr, metric)
     odet, _, _, _ = GeneralizedPerturbedEquilibrium.ForceFreeStates.eulerlagrange_integration(ctrl, equil, ffit, intr)
-    vac = GeneralizedPerturbedEquilibrium.ForceFreeStates.free_run!(odet, ctrl, equil, ffit, intr)
+    vac = GeneralizedPerturbedEquilibrium.ForceFreeStates.free_run(odet, ctrl, equil, ffit, intr)
     return real(vac.et[1]), intr.numpert_total
 end
 
diff --git a/benchmarks/plot_xi_eigenmode.jl b/benchmarks/plot_xi_eigenmode.jl
index e7c675207..f39682caa 100644
--- a/benchmarks/plot_xi_eigenmode.jl
+++ b/benchmarks/plot_xi_eigenmode.jl
@@ -1,6 +1,6 @@
 # Plot the ξ(ψ) profile of the eigenmode with the highest total-energy eigenvalue.
 #
-# The total energy operator W = W_plasma + W_vacuum (free_run!, Free.jl); its eigenvectors are the
+# The total energy operator W = W_plasma + W_vacuum (free_run, Free.jl); its eigenvectors are the
 # free-boundary edge displacement patterns and the eigenvalues are δW. This picks the eigenmode with the
 # largest Re(eigenvalue) and reconstructs its radial profile by projecting the EL fundamental matrix
 # (integration/xi_psi) onto that edge eigenvector:  c = U_edge \ w,  ξ(ψ) = U(ψ)·c.
diff --git a/docs/src/stability.md b/docs/src/stability.md
index f74a41649..4e7e72b7b 100644
--- a/docs/src/stability.md
+++ b/docs/src/stability.md
@@ -311,7 +311,7 @@ ffit   = FFS.make_matrix(equil, intr, metric)
 # and always returns a 4-tuple (odet, propagators, chunks, S_at_surface_left).
 odet, _, _, _ = FFS.eulerlagrange_integration(ctrl, equil, ffit, intr)
 
-vac = FFS.free_run!(odet, ctrl, equil, ffit, intr)
+vac = FFS.free_run(odet, ctrl, equil, ffit, intr)
 println("Energy eigenvalue et[1] = ", real(vac.et[1]))
 ```
 
@@ -356,5 +356,5 @@ end
 
 - `docs/src/galerkin.md` — RDCON outer-region Galerkin Δ′ solver (part of this module)
 - `docs/src/equilibrium.md` — build the `PlasmaEquilibrium` object required by this module
-- `docs/src/vacuum.md` — vacuum response computed from the EL solution in `free_run!`
+- `docs/src/vacuum.md` — vacuum response computed from the EL solution in `free_run`
 - `docs/src/perturbed_equilibrium.md` — downstream singular coupling analysis using Δ'
diff --git a/docs/src/vacuum.md b/docs/src/vacuum.md
index 46ca5f44d..40b7c08fb 100644
--- a/docs/src/vacuum.md
+++ b/docs/src/vacuum.md
@@ -73,8 +73,8 @@ wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings(
     equal_arc_wall = true     # Use equal arc length spacing
 )
 
-# Compute vacuum response matrix
-wv, grri, xzpts = GeneralizedPerturbedEquilibrium.Vacuum.compute_vacuum_response(inputs, wall_settings)
+# Compute vacuum response; returns a VacuumResponse with wv, grri, grre, plasma_pts, wall_pts
+vac = GeneralizedPerturbedEquilibrium.Vacuum.compute_vacuum_response(inputs, wall_settings)
 ```
 
 ### Vacuum Field Calculation at Observation Points
diff --git a/examples/DIIID-like_ideal_example_IMAS/run_imas_example.jl b/examples/DIIID-like_ideal_example_IMAS/run_imas_example.jl
index df9641531..f7f6cfb9e 100644
--- a/examples/DIIID-like_ideal_example_IMAS/run_imas_example.jl
+++ b/examples/DIIID-like_ideal_example_IMAS/run_imas_example.jl
@@ -61,7 +61,7 @@ try
         TOML.print(io, config_imas);
     end
     result_imas = GPEC.main([tmpdir_imas]; dd=dd)
-    global et_imas = real(result_imas.vac_data.et[1])
+    global et_imas = real(result_imas.free_energies.et[1])
     global mpert_imas = result_imas.intr.mpert
     GPEC.write_imas(dd, result_imas)
     @assert dd.mhd_linear.time_slice[1].toroidal_mode[1].energy_perturbed ≈ et_imas
@@ -80,7 +80,7 @@ try
         TOML.print(io, config_gfile);
     end
     result_gfile = GPEC.main([tmpdir_gfile])
-    global et_gfile = real(result_gfile.vac_data.et[1])
+    global et_gfile = real(result_gfile.free_energies.et[1])
     global mpert_gfile = result_gfile.intr.mpert
 finally
     rm(tmpdir_gfile; recursive=true)
diff --git a/src/ForceFreeStates/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl
index abcaa2bb1..69ea26cdc 100644
--- a/src/ForceFreeStates/EulerLagrange.jl
+++ b/src/ForceFreeStates/EulerLagrange.jl
@@ -253,7 +253,7 @@ function serial_eulerlagrange_integration(ctrl::ForceFreeStatesControl, equil::E
     end
     odet.nzero = evaluate_stability_criterion!(odet, equil.profiles)
 
-    # Undo Gaussian reduction to get true solution vectors (for free_run! eigenvector use)
+    # Undo Gaussian reduction to get true solution vectors
     transform_u!(odet, intr)
 
     return (odet, nothing, nothing, nothing)
diff --git a/src/ForceFreeStates/ForceFreeStatesStructs.jl b/src/ForceFreeStates/ForceFreeStatesStructs.jl
index a3d783812..7c2a70d87 100644
--- a/src/ForceFreeStates/ForceFreeStatesStructs.jl
+++ b/src/ForceFreeStates/ForceFreeStatesStructs.jl
@@ -256,7 +256,7 @@ gpec.toml.
   - `use_riccati::Bool` - Use the dual Riccati reformulation S = U₁·U₂⁻¹ instead of the standard U₁/U₂ ODE. Reduces stiffness for faster integration. See Glasser (2018) Phys. Plasmas 25, 032507.
   - `use_parallel::Bool` - Parallel fundamental matrix (propagator) integration using `Threads.@threads`. Each chunk is integrated independently from identity IC and assembled serially. Requires `singfac_min != 0`. Uses the same chunk bounds as the standard path but sub-divides chunks for load balancing. Crossings use the Riccati-style algorithm (no Gaussian reduction).
   - `parallel_threads::Int` - Cap on the number of threads the parallel BVP uses. **Default `2`** parallelises the FM chunks across two threads (the BVP has ~10 chunks; 2 threads is enough to amortize them — speedup saturates here, raising to 4 adds scheduling overhead). Set `parallel_threads = 1` to run the FM chunks SERIALLY (no `Threads.@threads`), which is bit-deterministic and immune to the thread-schedule sensitivity that can cause intermittent BVP divergence on numerically delicate equilibria. The parallel path produces bit-identical Δ′ across thread counts; `parallel_threads = 2` is about 20% faster than serial and saturates the speedup. If a parallel run diverges, drop to `parallel_threads = 1` rather than switching `use_parallel = false` — the latter is silently wrong. Capped at `Threads.nthreads()`.
-  - `populate_dense_xi::Bool` - When `use_parallel = true`, append a serial Euler-Lagrange pass at the end of the propagator BVP and let it replace the `odet` returned to the main pipeline.  This populates `u_store` / `du_store` / `xi_s_store` densely in the axis (EL) basis — the only convention the PerturbedEquilibrium / FieldReconstruction downstream code consumes correctly.  Without it the parallel path stores only chunk-endpoint Riccati S matrices with diagnostic derivatives (see Riccati.jl docstring caveats), and HDF5 `integration/xi_psi`/`dxi_psi`/`xi_s` are unusable.  Δ' (`singular/delta_prime_matrix`) is computed from the parallel BVP and is bit-identical between `populate_dense_xi=true` and `false`.  Energies (`vacuum/ep`/`ev`/`et`) are computed by `free_run!` from `odet`, so with `populate_dense_xi=true` they match what a pure serial run (`use_parallel=false`) would produce; with `populate_dense_xi=false` they use the parallel-pass Riccati `odet.u` instead (differs by the ~0.12 % Riccati-vs-axis algorithmic gap on DIIID-class cases).  **Default `false`** to avoid paying the dense-pass cost on Δ'/vacuum/ideal-stability-only runs; **PerturbedEquilibrium-using configs must set `populate_dense_xi = true` explicitly** when `use_parallel = true` (otherwise PE silently reads Riccati-basis garbage).  Auto-disabled when `force_termination = true` regardless of the user setting, since the dense pass has no downstream consumer in that case.  Approximate cost when enabled: one extra serial EL integration (~1× the parallel BVP wall-clock for typical N).
+  - `populate_dense_xi::Bool` - When `use_parallel = true`, append a serial Euler-Lagrange pass at the end of the propagator BVP and let it replace the `odet` returned to the main pipeline.  This populates `u_store` / `du_store` / `xi_s_store` densely in the axis (EL) basis — the only convention the PerturbedEquilibrium / FieldReconstruction downstream code consumes correctly.  Without it the parallel path stores only chunk-endpoint Riccati S matrices with diagnostic derivatives (see Riccati.jl docstring caveats), and HDF5 `integration/xi_psi`/`dxi_psi`/`xi_s` are unusable.  Δ' (`singular/delta_prime_matrix`) is computed from the parallel BVP and is bit-identical between `populate_dense_xi=true` and `false`.  Energies (`vacuum/ep`/`ev`/`et`) are computed by `free_run` from `odet`, so with `populate_dense_xi=true` they match what a pure serial run (`use_parallel=false`) would produce; with `populate_dense_xi=false` they use the parallel-pass Riccati `odet.u` instead (differs by the ~0.12 % Riccati-vs-axis algorithmic gap on DIIID-class cases).  **Default `false`** to avoid paying the dense-pass cost on Δ'/vacuum/ideal-stability-only runs; **PerturbedEquilibrium-using configs must set `populate_dense_xi = true` explicitly** when `use_parallel = true` (otherwise PE silently reads Riccati-basis garbage).  Auto-disabled when `force_termination = true` regardless of the user setting, since the dense pass has no downstream consumer in that case.  Approximate cost when enabled: one extra serial EL integration (~1× the parallel BVP wall-clock for typical N).
   - `extended_precision_bvp::Bool` - When `true` (default), promote the Δ' BVP linear system to `Complex{Double64}` (~31 digits) for the LU solve and PEST3 combination. Guards against catastrophic cancellation in the PEST3 four-term combination (dp_raw entries can be 10⁴–10⁵× larger than the result; the imaginary part of off-diagonal Δ' is particularly sensitive). Disabling (`false`) saves ~1.5–2× the BVP solve time but on DIIID-class equilibria the imaginary Δ' components can drift by factors of 2–5×; only disable for performance experiments on cases where Float64 has been validated against Double64.
 """
 @kwdef struct ForceFreeStatesControl
@@ -411,52 +411,39 @@ end
 FourFitVars(mpert::Int, numpert_total::Int) = FourFitVars(; mpert, numpert_total)
 
 """
-    VacuumData
+    FreeBoundaryResult
 
-A struct containing relevant data from the vacuum calculation.
-Populated in `Free.jl`.
+Result of the free-boundary calculation, returned by `free_run`. All matrices are in the ξ Fourier
+basis and are `numpert_total × numpert_total`; the energies are generalized (W, N) pencil values,
+power-normalized and invariant to the working (Jacobian) coordinate.
 
 ## Fields
 
-  - `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 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.
+  - `wt::Matrix{ComplexF64}` - Eigenvector matrix of W·v = λ·N·v. Columns are eigenmodes sorted most-unstable first, normalized to unit power norm v†·N·v = 1.
+  - `wt0::Matrix{ComplexF64}` - Total-energy matrix W = wp + wv before diagonalisation
+  - `wp::Matrix{ComplexF64}` - Plasma energy matrix
+  - `wv::Matrix{ComplexF64}` - Vacuum energy matrix, singfac-scaled at `qlim`
   - `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`
+  - `et::Vector{ComplexF64}` - Total energy eigenvalues of the pencil (W, N); et = ep + ev per mode
+  - `n_tor_idx::Vector{Int}` - 0-based toroidal mode number index of each sorted eigenvalue
   - `vacuum_eigenvalue::Float64` - Least stable (minimum) eigenvalue of the pencil (wv, N), clamped to zero
-  - `grri::Array{ComplexF64, 2}` - Interior Green's function matrices (2 * mthvac * nzvac × numpert_total)
-  - `grre::Array{ComplexF64, 2}` - Exterior Green's function matrices (2 * mthvac * nzvac × numpert_total)
-  - `plasma_pts::Array{Float64, 3}` - Cartesian coordinates of plasma points, shape (mthvac * nzvac) × 3 for (x, y, z)
-  - `wall_pts::Array{Float64, 3}` - Cartesian coordinates of wall points, shape (mthvac * nzvac) × 3 for (x, y, z)
+  - `plasma_pts`, `wall_pts::Matrix{Float64}` - Cartesian (x, y, z) surface coordinates, `numpoints × 3`, retained for HDF5 output
 """
-@kwdef mutable struct VacuumData
-    numpoints::Int
-    numpert_total::Int
-    mthvac::Int # this is only needed to not break GPEC functionality currently
-
-    wt::Array{ComplexF64,2} = Array{ComplexF64}(undef, numpert_total, numpert_total)
-    wt0::Array{ComplexF64,2} = Array{ComplexF64}(undef, numpert_total, numpert_total)
-    wp::Array{ComplexF64,2} = Array{ComplexF64}(undef, numpert_total, numpert_total)
-    wv::Array{ComplexF64,2} = Array{ComplexF64}(undef, numpert_total, numpert_total)
-    ep::Vector{ComplexF64} = Vector{ComplexF64}(undef, numpert_total)
-    ev::Vector{ComplexF64} = Vector{ComplexF64}(undef, numpert_total)
-    et::Vector{ComplexF64} = Vector{ComplexF64}(undef, numpert_total)
-    n_tor_idx::Vector{Int} = zeros(Int, numpert_total)
-    vacuum_eigenvalue::Float64 = NaN
-    grri::Array{ComplexF64,2} = Array{ComplexF64}(undef, 2 * numpoints, numpert_total)
-    grre::Array{ComplexF64,2} = Array{ComplexF64}(undef, 2 * numpoints, numpert_total)
-    plasma_pts::Array{Float64,2} = Array{Float64}(undef, numpoints, 3)
-    wall_pts::Array{Float64,2} = Array{Float64}(undef, numpoints, 3)
+struct FreeBoundaryResult
+    wt::Matrix{ComplexF64}
+    wt0::Matrix{ComplexF64}
+    wp::Matrix{ComplexF64}
+    wv::Matrix{ComplexF64}
+    ep::Vector{ComplexF64}
+    ev::Vector{ComplexF64}
+    et::Vector{ComplexF64}
+    n_tor_idx::Vector{Int}
+    vacuum_eigenvalue::Float64
+    plasma_pts::Matrix{Float64}
+    wall_pts::Matrix{Float64}
 end
 
-VacuumData(numpoints::Int, numpert_total::Int, mthvac::Int) = VacuumData(; numpoints, numpert_total, mthvac)
-
 """
 EdgeScanState
 
diff --git a/src/ForceFreeStates/Free.jl b/src/ForceFreeStates/Free.jl
index 42a4a36f4..51efa7610 100644
--- a/src/ForceFreeStates/Free.jl
+++ b/src/ForceFreeStates/Free.jl
@@ -9,7 +9,7 @@ Fourier band `jmat` (length 2·mpert−1, evaluated from the `ffit.jmats` spline
 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
+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.
 """
@@ -26,21 +26,44 @@ function power_norm_matrix!(Nmat::AbstractMatrix{ComplexF64}, jmat::AbstractVect
 end
 
 """
-    free_run!(odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal) -> VacuumData
+    normalize_eigenfunctions!(odet::OdeState, wt::AbstractMatrix{ComplexF64}, psio::Float64) -> OdeState
+
+Rescale the stored EL solution vectors in `odet.u_store`, `du_store`, and `xi_s_store` so the
+edge displacement matches the free-boundary eigenvectors `wt` (scaled by `2π·psio·1e-3`).
+Modifies `odet` in place. Call after `free_run` when downstream code consumes the stored ξ profiles.
+"""
+@with_pool pool function normalize_eigenfunctions!(odet::OdeState, wt::AbstractMatrix{ComplexF64}, psio::Float64)
+    N = size(wt, 1)
+    tmp_mat = zeros!(pool, ComplexF64, N, N)
+    coeffs = odet.u[:, :, 1, end] \ (wt .* (2π * psio * 1e-3))
+    @views for istep in 1:odet.step
+        mul!(tmp_mat, odet.u_store[:, :, 1, istep], coeffs)
+        odet.u_store[:, :, 1, istep] .= tmp_mat
+        mul!(tmp_mat, odet.u_store[:, :, 2, istep], coeffs)
+        odet.u_store[:, :, 2, istep] .= tmp_mat
+        mul!(tmp_mat, odet.du_store[:, :, 1, istep], coeffs)
+        odet.du_store[:, :, 1, istep] .= tmp_mat
+        mul!(tmp_mat, odet.du_store[:, :, 2, istep], coeffs)
+        odet.du_store[:, :, 2, istep] .= tmp_mat
+        mul!(tmp_mat, odet.xi_s_store[:, :, istep], coeffs)
+        odet.xi_s_store[:, :, istep] .= tmp_mat
+    end
+    return odet
+end
+
+"""
+    free_run(odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal) -> FreeBoundaryResult
 
 Compute the free boundary energies using the Julia port of the VACUUM code. Performs the same function as `free_run`
-in the Fortran code, except now all data is passed in memory instead of via files. This
-modifies `odet` in place to normalize the eigenfunctions stored in `u_store`, `du_store`, and `xi_s_store`,
-and returns a `VacuumData` struct containing the data needed for perturbed equilibrium calculations
-and data dumping.
+in the Fortran code.
+
+Returns a `FreeBoundaryResult` struct containing the data needed for perturbed equilibrium
+calculations and data dumping.
 """
-@with_pool pool function free_run!(odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal)
+@with_pool pool function free_run(odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal)
 
     # Initializations and allocations
     (; mpert, mlow, mhigh, numpert_total, psilim, qlim, npert, nlow, nhigh, wall_settings) = intr
-    vac_data = VacuumData(ctrl.mthvac * ctrl.nzvac, intr.numpert_total, ctrl.mthvac)
-    etemp = zeros!(pool, ComplexF64, numpert_total)
-    wp = zeros!(pool, ComplexF64, numpert_total, numpert_total)
     wpt = zeros!(pool, ComplexF64, numpert_total, numpert_total)
     wvt = zeros!(pool, ComplexF64, numpert_total, numpert_total)
     tmp_mat = zeros!(pool, ComplexF64, numpert_total, numpert_total)
@@ -49,15 +72,17 @@ and data dumping.
     dV_dpsi = equil.profiles.dVdpsi_spline(psilim)
 
     # Compute plasma response matrix W = U₂ * U₁⁻¹
+    wp = zeros(ComplexF64, numpert_total, numpert_total)
     @views wp .= (odet.u[:, :, 2] / odet.u[:, :, 1]) ./ equil.psio^2
 
-    # Compute vacuum response matrix in-place (handles 2D single-n, 2D multi-n block-diagonal, and 3D)
+    # Compute vacuum response (handles 2D single-n, 2D multi-n block-diagonal, and 3D)
     vac_inputs = Vacuum.VacuumInput(equil, psilim, ctrl.mthvac, ctrl.nzvac, mlow:mhigh, nlow:nhigh)
-    Vacuum.compute_vacuum_response!(vac_data, vac_inputs, wall_settings)
+    vac = Vacuum.compute_vacuum_response(vac_inputs, wall_settings)
+    wv = vac.wv
 
     # Scale by (m - n*q)(m' - n'*q) [Chance Phys. Plasmas 1997 2161 eq. 126]
     singfac = vec((mlow:mhigh) .- qlim .* (nlow:nhigh)')
-    vac_data.wv .*= singfac .* singfac'
+    wv .*= singfac .* singfac'
 
     # 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.
@@ -68,81 +93,64 @@ and data dumping.
 
     # 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
+    vacuum_eigenvalue = max(0.0, minimum(real.(eigvals(Hermitian(wv), Hermitian(Nmat)))))
 
     # 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, Nmat)
-    vac_data.et .= Ev.values
-    eindex = sortperm(real.(vac_data.et); rev=true)
+    wt = wp .+ wv
+    wt0 = copy(wt)
+    Ev = eigen(wt, Nmat)
+    et = Ev.values
+    eindex = sortperm(real.(et); rev=true)
 
-    etemp .= vac_data.et
     # Rearrange wt columns for ascending real eigenvalues (most unstable first)
+    etemp = copy(et)
+    n_tor_idx = zeros(Int, numpert_total)
     for ipert in 1:numpert_total
         orig = eindex[numpert_total+1-ipert]
-        vac_data.wt[:, ipert] .= Ev.vectors[:, orig]
-        vac_data.et[ipert] = etemp[orig]
+        wt[:, ipert] .= Ev.vectors[:, orig]
+        et[ipert] = etemp[orig]
         # Store which n this eigenvector corresponds to (needed to write IMAS data)
         # This relies on the block diagonal matrix structure due to n decoupling in tokamaks
         imax = argmax(abs.(Ev.vectors[:, orig]))
-        vac_data.n_tor_idx[ipert] = (imax - 1) ÷ mpert
+        n_tor_idx[ipert] = (imax - 1) ÷ mpert
     end
 
     # 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
-        v = @view vac_data.wt[:, isol]
+        v = @view wt[:, isol]
         v ./= sqrt(real(dot(v, Nmat, v)))
     end
 
     # Normalize phase
     imax = 0
     for isol in 1:numpert_total
-        imax = argmax(abs.(vac_data.wt[:, isol]))
-        phase = abs(vac_data.wt[imax, isol]) / vac_data.wt[imax, isol]
-        vac_data.wt[:, isol] .*= phase
+        imax = argmax(abs.(wt[:, isol]))
+        phase = abs(wt[imax, isol]) / wt[imax, isol]
+        wt[:, isol] .*= phase
     end
 
     # Project W_p and W_v into the eigenmode basis. Diagonal entries give
     # the plasma/vacuum energy split for each mode: et[i] = ep[i] + ev[i]
-    mul!(tmp_mat, wp, vac_data.wt)
-    mul!(wpt, vac_data.wt', tmp_mat)
-    mul!(tmp_mat, vac_data.wv, vac_data.wt)
-    mul!(wvt, vac_data.wt', tmp_mat)
-    vac_data.ep .= diag(wpt)
-    vac_data.ev .= diag(wvt)
-
-    # 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
-        mul!(tmp_mat, odet.u_store[:, :, 1, istep], coeffs)
-        odet.u_store[:, :, 1, istep] .= tmp_mat
-        mul!(tmp_mat, odet.u_store[:, :, 2, istep], coeffs)
-        odet.u_store[:, :, 2, istep] .= tmp_mat
-        mul!(tmp_mat, odet.du_store[:, :, 1, istep], coeffs)
-        odet.du_store[:, :, 1, istep] .= tmp_mat
-        mul!(tmp_mat, odet.du_store[:, :, 2, istep], coeffs)
-        odet.du_store[:, :, 2, istep] .= tmp_mat
-        mul!(tmp_mat, odet.xi_s_store[:, :, istep], coeffs)
-        odet.xi_s_store[:, :, istep] .= tmp_mat
-    end
+    mul!(tmp_mat, wp, wt)
+    mul!(wpt, wt', tmp_mat)
+    mul!(tmp_mat, wv, wt)
+    mul!(wvt, wt', tmp_mat)
+    ep = diag(wpt)
+    ev = diag(wvt)
 
     # Write energies to screen
     if ctrl.verbose
         @info "Least Stable Eigenmode Energies:\n" *
-              "  Plasma = $((@sprintf "%+.3e %+.3ei" real(vac_data.ep[1]) imag(vac_data.ep[1])))\n" *
-              "  Vacuum = $((@sprintf "%+.3e %+.3ei" real(vac_data.ev[1]) imag(vac_data.ev[1])))\n" *
-              "  Total  = $((@sprintf "%+.3e %+.3ei" real(vac_data.et[1]) imag(vac_data.et[1])))"
+              "  Plasma = $((@sprintf "%+.3e %+.3ei" real(ep[1]) imag(ep[1])))\n" *
+              "  Vacuum = $((@sprintf "%+.3e %+.3ei" real(ev[1]) imag(ev[1])))\n" *
+              "  Total  = $((@sprintf "%+.3e %+.3ei" real(et[1]) imag(et[1])))"
     end
 
-    return vac_data
+    return FreeBoundaryResult(wt, wt0, wp, wv, ep, ev, et, n_tor_idx, vacuum_eigenvalue, vac.plasma_pts, vac.wall_pts)
 end
 
 """
@@ -177,8 +185,8 @@ q-window minimum.
 
         # Compute raw vacuum matrix at the actual scan psi (singfac NOT applied; free_compute_total applies it analytically)
         vac_inputs = Vacuum.VacuumInput(equil, psi_array[i], ctrl.mthvac, ctrl.nzvac, intr.mlow:intr.mhigh, intr.nlow:intr.nhigh)
-        wv, _, _, _, _ = Vacuum.compute_vacuum_response(vac_inputs, intr.wall_settings)
-        @views wv_array[i, :, :] .= wv
+        vac = Vacuum.compute_vacuum_response(vac_inputs, intr.wall_settings)
+        @views wv_array[i, :, :] .= vac.wv
     end
 
     # Flatten 3D array to (npsi+1 × numpert_total^2) for series interpolant
@@ -234,7 +242,7 @@ wv matrix spline to `free_compute_wv_spline` and pass it in `odet.edge_scan.wvma
     power_norm_matrix!(Nmat, jmat_local, intr.mpert, intr.npert, dV_dpsi)
 
     # Total energy matrix and generalized eigen-decomposition of the pencil (W, N) — the
-    # eigenvalues are power-normalized, Jacobian-invariant mode energies (see free_run!)
+    # eigenvalues are power-normalized, Jacobian-invariant mode energies
     wt .= wp .+ wv
     Ev = eigen(wt, Nmat)
 
diff --git a/src/ForceFreeStates/Galerkin/GalerkinSolve.jl b/src/ForceFreeStates/Galerkin/GalerkinSolve.jl
index ace955237..8389ad9bc 100644
--- a/src/ForceFreeStates/Galerkin/GalerkinSolve.jl
+++ b/src/ForceFreeStates/Galerkin/GalerkinSolve.jl
@@ -48,17 +48,17 @@ end
 
 """
     galerkin_solve(ctrl::ForceFreeStatesControl, equil, ffit::FourFitVars,
-                   intr::ForceFreeStatesInternal; vac_data=nothing) -> GalerkinResult
+                   intr::ForceFreeStatesInternal; wv=nothing) -> GalerkinResult
 
 Compute the outer-region Δ′ matching matrix by the singular Galerkin method. Port of `gal_solve`
 (gal.f). Single toroidal mode only (`intr.npert == 1`). Returns a `GalerkinResult`; if there
 are no resonant surfaces in the domain it returns an empty result.
 
-`vac_data` (a `VacuumData` from `free_run!`) supplies the free-boundary edge term
-`wv_edge = vac_data.wv · psio²`; pass `nothing` for a fixed-boundary edge.
+`wv` is the vacuum energy matrix from `free_run`, which supplies the free-boundary edge term
+`wv_edge = wv · psio²`; pass `nothing` for a fixed-boundary edge.
 """
 function galerkin_solve(ctrl::ForceFreeStatesControl, equil, ffit::FourFitVars,
-    intr::ForceFreeStatesInternal; vac_data=nothing)
+    intr::ForceFreeStatesInternal; wv=nothing)
 
     intr.npert == 1 || error("galerkin_solve: only single-n (npert == 1) is supported")
     ctrl.gal_solver in ("LU", "cholesky") || error("galerkin_solve: gal_solver must be \"LU\" or \"cholesky\"")
@@ -121,12 +121,12 @@ function galerkin_solve(ctrl::ForceFreeStatesControl, equil, ffit::FourFitVars,
     ws.rhs = zeros(ComplexF64, ws.ndim, nsol)
     ws.sol = zeros(ComplexF64, ws.ndim, nsol)
 
-    # Free-boundary edge term wvac·psio² (vac_data.wv is already singfac-scaled at qlim in free_run!).
+    # Free-boundary edge term wvac·psio² (wv is already singfac-scaled at qlim in free_run).
     # rpec passes wv_edge=nothing; gal_set_boundary! then applies the identity edge AND injects the coil
     # unit sources (its three-way branch). The vacuum block is only built for the non-rpec free case.
     wv_edge = nothing
-    if ctrl.vac_flag && vac_data !== nothing && ncoil == 0
-        wv_edge = Matrix{ComplexF64}(vac_data.wv .* equil.psio^2)
+    if ctrl.vac_flag && wv !== nothing && ncoil == 0
+        wv_edge = Matrix{ComplexF64}(wv .* equil.psio^2)
     end
 
     gal_make_arrays!(ws, ctrl, equil, ffit, intr, asymps, sings, nn, wv_edge)
diff --git a/src/ForceFreeStates/Riccati.jl b/src/ForceFreeStates/Riccati.jl
index f67214c64..66d1361e7 100644
--- a/src/ForceFreeStates/Riccati.jl
+++ b/src/ForceFreeStates/Riccati.jl
@@ -1672,7 +1672,7 @@ function parallel_eulerlagrange_integration(
 
     chunks, propagators = _handle_edge_dW_scan!(odet, chunks, propagators, ctrl, equil, ffit, intr)
 
-    # compute_delta_prime_matrix! is called from the main pipeline (after free_run!) so
+    # compute_delta_prime_matrix! is called from the main pipeline (after free_run) so
     # that vacuum response wv is available for the edge BC. With self-consistent truncation,
     # the propagators/chunks returned here match intr.psilim exactly, so Δ' is well-defined
     # for both truncate_at_dW_peak=false (full domain) and =true (peak).
@@ -1907,8 +1907,9 @@ end
 Replace the propagator-BVP's `odet` with a fresh serial-EL `odet` that has
 dense `u_store` / `du_store` populated in axis basis (the PerturbedEquilibrium
 convention).  The caller's `odet` is fully replaced by the fresh one because
-`free_run!` downstream uses `odet.u[:,:,1,end]` to normalize `odet.u_store`,
-so both must be in the same basis.  The parallel BVP results that survive
+`free_run` / `normalize_eigenfunctions!` downstream use `odet.u[:,:,1,end]` to
+normalize `odet.u_store`, so both must be in the same basis.  The parallel BVP
+results that survive
 downstream are stored in `intr` (psilim/qlim, sing[*].delta_prime, …) and in
 the externally-returned `propagators` / `chunks` / `S_at_surface_left` —
 none of those live on `odet`, so replacing `odet` is safe.
diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl
index fde4148d7..9e6cd0aa2 100755
--- a/src/GeneralizedPerturbedEquilibrium.jl
+++ b/src/GeneralizedPerturbedEquilibrium.jl
@@ -62,12 +62,12 @@ export Analysis
 include("Rerun.jl")
 
 # Import ForceFreeStates types and functions needed for main
-using .ForceFreeStates: ForceFreeStatesInternal, ForceFreeStatesControl, DebugSettings, VacuumData, OdeState, FourFitVars
+using .ForceFreeStates: ForceFreeStatesInternal, ForceFreeStatesControl, DebugSettings, FreeBoundaryResult, OdeState, FourFitVars
 using .ForceFreeStates: sing_lim!, sing_min!, sing_find!, resist_eval_all!, resist_geometry, ResistGeometry
 using .ForceFreeStates: compute_local_stability, compute_ballooning_stability!, ballooning_alpha_boundary, ballooning_alpha_boundaries
 using .ForceFreeStates: make_metric, make_matrix, make_kinetic_matrix
 using .ForceFreeStates: find_kinetic_singular_surfaces!
-using .ForceFreeStates: eulerlagrange_integration, free_run!
+using .ForceFreeStates: eulerlagrange_integration, free_run, normalize_eigenfunctions!
 using .ForceFreeStates: galerkin_solve, write_galerkin!, GalerkinResult, gal_matched_odestate
 
 const _DEPRECATED_FFS_KEYS = ("mer_flag", "force_wv_symmetry", "ode_flag", "cyl_flag", "mat_flag")
@@ -446,14 +446,16 @@ function main_from_inputs(
         @warn "Fixed-boundary mode unstable for n = $nstring"
     end
 
-    # Compute free boundary energies
+    # Compute free boundary energies.
+    free_energies = nothing
     if ctrl.vac_flag && !(ctrl.ksing > 0 && ctrl.ksing <= intr.msing + 1)
         if ctrl.verbose
             wall_desc = intr.wall_settings.shape == "nowall" ? "no wall" : intr.wall_settings.shape
             @info "Computing free boundary energies ($wall_desc)"
         end
-        vac_data = free_run!(odet, ctrl, equil, ffit, intr)
-        if real(vac_data.et[1]) < 0
+        free_energies = free_run(odet, ctrl, equil, ffit, intr)
+        normalize_eigenfunctions!(odet, free_energies.wt, equil.psio)
+        if real(free_energies.et[1]) < 0
             if ctrl.verbose
                 @warn "Free-boundary mode unstable for n = $nstring"
             end
@@ -464,13 +466,13 @@ function main_from_inputs(
         end
 
         # Compute inter-surface Δ' matrix (STRIDE BVP) using vacuum edge BC.
-        # Requires propagators from parallel FM path and wv from free_run!.
+        # Requires propagators from parallel FM path and wv from free_run.
         if ctrl.kinetic_factor == 0 && intr.msing > 0 && fm_propagators !== nothing
             if ctrl.verbose
                 @info "Computing Δ' matrix (STRIDE BVP with vacuum coupling)"
             end
             ForceFreeStates.compute_delta_prime_matrix!(intr, fm_propagators, fm_chunks;
-                wv=vac_data.wv, psio=equil.psio, debug=ctrl.verbose,
+                wv=free_energies.wv, psio=equil.psio, debug=ctrl.verbose,
                 S_at_surface_left=fm_S_left,
                 ctrl=ctrl, equil=equil, ffit=ffit)
         end
@@ -480,7 +482,7 @@ function main_from_inputs(
     gal_data = nothing
     if ctrl.gal_flag
         gal_start = time()
-        gal_data = galerkin_solve(ctrl, equil, ffit, intr; vac_data=ctrl.vac_flag ? vac_data : nothing)
+        gal_data = galerkin_solve(ctrl, equil, ffit, intr; wv=free_energies !== nothing ? free_energies.wv : nothing)
         @info "Galerkin solve completed in $(@sprintf("%.3f", time() - gal_start)) s"
     end
 
@@ -490,7 +492,7 @@ function main_from_inputs(
             equil,
             intr,
             odet,
-            ctrl.vac_flag ? vac_data : nothing,
+            free_energies,
             ffit,
             git_version,
             inputs,
@@ -543,7 +545,7 @@ function main_from_inputs(
         slayer_result = _run_slayer_stage(nothing)
         @info "\n$_BANNER\n  GPEC completed successfully in $(@sprintf("%.3f", time() - total_start)) s\n$_BANNER"
         return (ctrl=ctrl, equil=equil, intr=intr, ffit=ffit, odet=odet,
-            vac_data=ctrl.vac_flag ? vac_data : nothing,
+            free_energies=free_energies,
             slayer=slayer_result)
     end
 
@@ -604,9 +606,9 @@ function main_from_inputs(
         end
 
         # Run perturbed equilibrium calculations
-        # Pass vac_data and intr for response matrix calculations
+        # Free-boundary wt0 drives the plasma inductance; mthvac sizes the Green's-function solves
         pe_state = PerturbedEquilibrium.compute_perturbed_equilibrium(
-            equil, pe_odet, ctrl.vac_flag ? vac_data : nothing, intr, ft_ctrl, pe_ctrl, pe_intr,
+            equil, pe_odet, free_energies !== nothing ? free_energies.wt0 : nothing, ctrl.mthvac, intr, ft_ctrl, pe_ctrl, pe_intr,
             metric, ffit
         )
 
@@ -677,7 +679,7 @@ function main_from_inputs(
     # TODO: Do not allow perturbed equilibrium calculations if zero crossings are found
 
     return (ctrl=ctrl, equil=equil, intr=intr, ffit=ffit, odet=odet,
-        vac_data=ctrl.vac_flag ? vac_data : nothing,
+        free_energies=free_energies,
         slayer=slayer_result)
 
 end
@@ -700,7 +702,7 @@ function write_outputs_to_HDF5(
     equil::Equilibrium.PlasmaEquilibrium,
     intr::ForceFreeStatesInternal,
     odet::OdeState,
-    vac_data::Union{VacuumData,Nothing},
+    free_energies::Union{FreeBoundaryResult,Nothing},
     ffit::Union{FourFitVars,Nothing}=nothing,
     git_version::String="unknown",
     inputs::Union{Nothing,Dict{String,Any}}=nothing,
@@ -906,23 +908,23 @@ function write_outputs_to_HDF5(
         # 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
+        out_h5["FreeBoundaryStability/W_freeboundary"] = free_energies !== nothing ? free_energies.wt0 : ComplexF64[]
+        out_h5["FreeBoundaryStability/W_plasma"] = free_energies !== nothing ? free_energies.wp : ComplexF64[]
+        out_h5["FreeBoundaryStability/W_vacuum"] = free_energies !== nothing ? free_energies.wv : ComplexF64[]
+        out_h5["FreeBoundaryStability/W_freeboundary_eigenmodes"] = free_energies !== nothing ? free_energies.wt : ComplexF64[]
+        out_h5["FreeBoundaryStability/eigenmode_energies"] = free_energies !== nothing ? free_energies.et : ComplexF64[]
+        out_h5["FreeBoundaryStability/eigenmode_plasma_energies"] = free_energies !== nothing ? free_energies.ep : ComplexF64[]
+        out_h5["FreeBoundaryStability/eigenmode_vacuum_energies"] = free_energies !== nothing ? free_energies.ev : ComplexF64[]
+        out_h5["FreeBoundaryStability/vacuum_eigenvalue"] = free_energies !== nothing ? free_energies.vacuum_eigenvalue : NaN
 
         # Cartesian surface point clouds used downstream for visualisation and
         # perturbed-equilibrium plotting.
-        out_h5["SurfaceGeometries/Plasma/x"] = ctrl.vac_flag ? vac_data.plasma_pts[:, 1] : Float64[]
-        out_h5["SurfaceGeometries/Plasma/y"] = ctrl.vac_flag ? vac_data.plasma_pts[:, 2] : Float64[]
-        out_h5["SurfaceGeometries/Plasma/z"] = ctrl.vac_flag ? vac_data.plasma_pts[:, 3] : Float64[]
-        out_h5["SurfaceGeometries/Wall/x"] = ctrl.vac_flag ? vac_data.wall_pts[:, 1] : Float64[]
-        out_h5["SurfaceGeometries/Wall/y"] = ctrl.vac_flag ? vac_data.wall_pts[:, 2] : Float64[]
-        out_h5["SurfaceGeometries/Wall/z"] = ctrl.vac_flag ? vac_data.wall_pts[:, 3] : Float64[]
+        out_h5["SurfaceGeometries/Plasma/x"] = free_energies !== nothing ? free_energies.plasma_pts[:, 1] : Float64[]
+        out_h5["SurfaceGeometries/Plasma/y"] = free_energies !== nothing ? free_energies.plasma_pts[:, 2] : Float64[]
+        out_h5["SurfaceGeometries/Plasma/z"] = free_energies !== nothing ? free_energies.plasma_pts[:, 3] : Float64[]
+        out_h5["SurfaceGeometries/Wall/x"] = free_energies !== nothing ? free_energies.wall_pts[:, 1] : Float64[]
+        out_h5["SurfaceGeometries/Wall/y"] = free_energies !== nothing ? free_energies.wall_pts[:, 2] : Float64[]
+        out_h5["SurfaceGeometries/Wall/z"] = free_energies !== nothing ? free_energies.wall_pts[:, 3] : Float64[]
 
         # Write kinetic parameters when kinetic mode is enabled
         if ctrl.kinetic_factor > 0
@@ -1012,9 +1014,9 @@ receives the correct least-stable δW regardless of how modes are interleaved in
 The `result` argument is the named tuple returned by `main`.
 """
 function write_imas(dd, result)
-    result.vac_data === nothing && return
+    result.free_energies === nothing && return
 
-    vac_data = result.vac_data
+    free_energies = result.free_energies
     intr = result.intr
 
     # Top-level metadata
@@ -1029,10 +1031,10 @@ function write_imas(dd, result)
     # n_tor_idx[i] (0-based) identifies which n-block eigenvalue i belongs to.
     resize!(ts.toroidal_mode, intr.npert)
     for j in 0:(intr.npert-1)
-        n_indices = findall(==(j), vac_data.n_tor_idx) # indices of eigenvalues in the j-th n-block
+        n_indices = findall(==(j), free_energies.n_tor_idx) # indices of eigenvalues in the j-th n-block
         mode = ts.toroidal_mode[j+1]
         mode.n_tor = intr.nlow + j
-        mode.energy_perturbed = minimum(real.(vac_data.et[n_indices])) # least-stable energy for this n-toroidal mode
+        mode.energy_perturbed = minimum(real.(free_energies.et[n_indices])) # least-stable energy for this n-toroidal mode
     end
 
     return dd
diff --git a/src/PerturbedEquilibrium/PerturbedEquilibrium.jl b/src/PerturbedEquilibrium/PerturbedEquilibrium.jl
index 659f979ba..5de559c58 100644
--- a/src/PerturbedEquilibrium/PerturbedEquilibrium.jl
+++ b/src/PerturbedEquilibrium/PerturbedEquilibrium.jl
@@ -11,7 +11,7 @@ using FastInterpolations
 # Import parent modules
 import ..Equilibrium
 import ..ForceFreeStates
-import ..ForceFreeStates: OdeState, VacuumData, ForceFreeStatesInternal, FourFitVars, MetricData
+import ..ForceFreeStates: OdeState, ForceFreeStatesInternal, FourFitVars, MetricData
 import ..Vacuum
 import ..ForcingTerms
 import ..ForcingTerms: ForcingMode, CoilSet, load_forcing_data!, convert_forcing_normalization!
@@ -39,7 +39,7 @@ export write_outputs_to_HDF5
 
 """
     compute_perturbed_equilibrium(
-        equil, ForceFreeStates_results, vac_data, ffs_intr,
+        equil, ForceFreeStates_results, wt0, mthvac, ffs_intr,
         ft_ctrl, ctrl, intr, metric, ffit
     )::PerturbedEquilibriumState
 
@@ -52,7 +52,8 @@ coupling metrics.
 
   - `equil`: Equilibrium solution from Equilibrium module
   - `ForceFreeStates_results`: Stability calculation results from ForceFreeStates module
-  - `vac_data`: Vacuum response data from ForceFreeStates free boundary calculation
+  - `wt0`: Free-boundary total-energy matrix W = wp + wv from `free_run`, or `nothing` when the free-boundary calculation was not run (response and singular coupling are then skipped)
+  - `mthvac`: Vacuum poloidal grid resolution to reuse for the Green's-function solves
   - `ffs_intr`: ForceFreeStates internal state with mode information
   - `ft_ctrl`: Forcing terms control parameters from [ForcingTerms] section
   - `ctrl`: Control parameters from [PerturbedEquilibrium] section
@@ -67,7 +68,8 @@ coupling metrics.
 function compute_perturbed_equilibrium(
     equil::Equilibrium.PlasmaEquilibrium,
     ForceFreeStates_results::OdeState,
-    vac_data::Union{VacuumData,Nothing},
+    wt0::Union{Matrix{ComplexF64},Nothing},
+    mthvac::Int,
     ffs_intr::ForceFreeStates.ForceFreeStatesInternal,
     ft_ctrl::ForcingTerms.ForcingTermsControl,
     ctrl::PerturbedEquilibriumControl,
@@ -119,19 +121,19 @@ function compute_perturbed_equilibrium(
 
     # Step 2: Compute plasma response
     if ctrl.compute_response
-        if vac_data === nothing
+        if wt0 === nothing
             @warn "Vacuum data not available. Skipping plasma response calculation. Set vac_flag=true in [ForceFreeStates] section."
         else
-            compute_plasma_response!(state, equil, ForceFreeStates_results, vac_data, ffs_intr, intr, ctrl, metric, ffit)
+            compute_plasma_response!(state, equil, ForceFreeStates_results, wt0, mthvac, ffs_intr, intr, ctrl, metric, ffit)
         end
     end
 
     # Step 3: Compute singular coupling metrics
     if ctrl.compute_singular_coupling
-        if vac_data === nothing
+        if wt0 === nothing
             @warn "Vacuum data not available. Skipping singular coupling calculation. Set vac_flag=true in [ForceFreeStates] section."
         else
-            compute_singular_coupling_metrics!(state, equil, ForceFreeStates_results, vac_data, ffs_intr, intr, ctrl, ffit)
+            compute_singular_coupling_metrics!(state, equil, ForceFreeStates_results, mthvac, ffs_intr, intr, ctrl, ffit)
         end
     end
 
diff --git a/src/PerturbedEquilibrium/Response.jl b/src/PerturbedEquilibrium/Response.jl
index ba2827246..e501bad29 100644
--- a/src/PerturbedEquilibrium/Response.jl
+++ b/src/PerturbedEquilibrium/Response.jl
@@ -1,6 +1,6 @@
 """
     compute_plasma_response!(
-        state, equil, ForceFreeStates_results, vac_data, ffs_intr,
+        state, equil, ForceFreeStates_results, wt0, mthvac, ffs_intr,
         intr, ctrl, metric, ffit
     )
 
@@ -18,7 +18,8 @@ function compute_plasma_response!(
     state::PerturbedEquilibriumState,
     equil::Equilibrium.PlasmaEquilibrium,
     ForceFreeStates_results::OdeState,
-    vac_data::VacuumData,
+    wt0::Matrix{ComplexF64},
+    mthvac::Int,
     ffs_intr::ForceFreeStatesInternal,
     intr::PerturbedEquilibriumInternal,
     ctrl::PerturbedEquilibriumControl,
@@ -30,19 +31,19 @@ function compute_plasma_response!(
     end
 
     # Build flux matrix from ForceFreeStates eigenmodes [mode × eigenmode]
-    flux_matrix = build_flux_matrix(equil, ForceFreeStates_results, vac_data, ffs_intr)
+    flux_matrix = build_flux_matrix(equil, ForceFreeStates_results, ffs_intr)
 
     # Plasma inductance Lambda (wt0 formula, Fortran resp_induct_flag=TRUE default)
-    plasma_inductance = calc_plasma_inductance(vac_data, ffs_intr, equil.psio)
+    plasma_inductance = calc_plasma_inductance(wt0, ffs_intr, equil.psio)
 
     # Surface inductance L from Green's functions at psilim.
     # Requires a 2D (nzvac=1) vacuum response so rows are theta points only,
     nn = ffs_intr.nlow
-    vac_input_2d = Vacuum.VacuumInput(equil, ffs_intr.psilim, vac_data.mthvac, 1, ffs_intr.mlow:ffs_intr.mhigh, [nn])
+    vac_input_2d = Vacuum.VacuumInput(equil, ffs_intr.psilim, mthvac, 1, ffs_intr.mlow:ffs_intr.mhigh, [nn])
     wall_nowall = Vacuum.WallShapeSettings(; shape="nowall")
-    _, grri_2d_raw, grre_2d_raw, _, _ = Vacuum.compute_vacuum_response(vac_input_2d, wall_nowall)
-    grri_2d = Matrix{ComplexF64}(grri_2d_raw)
-    grre_2d = Matrix{ComplexF64}(grre_2d_raw)
+    vac_2d = Vacuum.compute_vacuum_response(vac_input_2d, wall_nowall)
+    grri_2d = Matrix{ComplexF64}(vac_2d.grri)
+    grre_2d = Matrix{ComplexF64}(vac_2d.grre)
     ν_vac = Vacuum.PlasmaGeometry(vac_input_2d).ν
     surface_inductance = compute_surface_inductance_from_greens(grri_2d, grre_2d, ffs_intr, nn, ν_vac)
     permeability = calc_permeability(plasma_inductance, surface_inductance)
diff --git a/src/PerturbedEquilibrium/ResponseMatrices.jl b/src/PerturbedEquilibrium/ResponseMatrices.jl
index a4b798a6d..791ed9c57 100644
--- a/src/PerturbedEquilibrium/ResponseMatrices.jl
+++ b/src/PerturbedEquilibrium/ResponseMatrices.jl
@@ -161,7 +161,6 @@ end
     build_flux_matrix(
         equil::Equilibrium.PlasmaEquilibrium,
         ForceFreeStates_results::OdeState,
-        vac_data::VacuumData,
         intr::ForceFreeStatesInternal
     )::Matrix{ComplexF64}
 
@@ -181,7 +180,6 @@ The flux matrix relates eigenmode displacements to vacuum poloidal flux:
 
   - `equil`: Equilibrium solution containing flux surfaces and q-profile
   - `ForceFreeStates_results`: ForceFreeStates ODE integration results containing eigenmodes
-  - `vac_data`: Vacuum response data from free boundary calculation
   - `intr`: ForceFreeStates internal state with mode information
 
 ## Returns
@@ -192,7 +190,6 @@ The flux matrix relates eigenmode displacements to vacuum poloidal flux:
 function build_flux_matrix(
     equil::Equilibrium.PlasmaEquilibrium,
     ForceFreeStates_results::OdeState,
-    vac_data::VacuumData,
     intr::ForceFreeStatesInternal
 )::Matrix{ComplexF64}
 
@@ -209,7 +206,7 @@ end
 
 """
     calc_plasma_inductance(
-        vac_data::VacuumData,
+        wt0::Matrix{ComplexF64},
         ffs_intr::ForceFreeStatesInternal,
         psio::Float64
     )::Matrix{ComplexF64}
@@ -221,13 +218,13 @@ Calculate plasma inductance matrix Λ using the wt0-based energy formula
 
 where t₁ = im/(χ₁·s_i·2π), t₂ = -im/(χ₁·s_j·2π), s_i = m_i - n·q_lim.
 
-Note: `vac_data.wt0` already contains singfac² factors (s_i·s_j) baked into the
-vacuum term via the scaling in `free_run!`. The t₁/t₂ factors divide by s_i·s_j,
+Note: `wt0` already contains singfac² factors (s_i·s_j) baked into the
+vacuum term via the scaling in `free_run`. The t₁/t₂ factors divide by s_i·s_j,
 correctly recovering the properly-normalized inductance.
 
 ## Arguments
 
-  - `vac_data`: Vacuum data containing wt0 (total energy matrix before eigenvector sorting)
+  - `wt0`: Total energy matrix W = wp + wv, before eigenvector sorting
   - `ffs_intr`: ForceFreeStates internal state with mode info (mlow, mpert, nlow, qlim)
   - `psio`: Total toroidal flux [Wb/rad] from equilibrium (equil.psio)
 
@@ -236,7 +233,7 @@ correctly recovering the properly-normalized inductance.
   - Plasma inductance matrix Lambda [numpert_total × numpert_total]
 """
 function calc_plasma_inductance(
-    vac_data::VacuumData,
+    wt0::Matrix{ComplexF64},
     ffs_intr::ForceFreeStatesInternal,
     psio::Float64
 )::Matrix{ComplexF64}
@@ -250,9 +247,9 @@ function calc_plasma_inductance(
     s = [((i-1) % ffs_intr.mpert + ffs_intr.mlow) - n * qlim for i in 1:mpert]
 
     # Fortran idcon_norm: wt0 = wt0/(mu0*2)*psio^2
-    # Julia's vac_data.wt0 is raw wp+wv; Fortran additionally scales by psio^2/(mu0*2)
+    # Julia's wt0 is raw wp+wv; Fortran additionally scales by psio^2/(mu0*2)
     mu0 = 4π * 1e-7
-    wt0_norm = vac_data.wt0 .* (psio^2 / (mu0 * 2))
+    wt0_norm = wt0 .* (psio^2 / (mu0 * 2))
 
     # Build temp2[i,j] = 2·t1_i·wt0[i,j]·t2_j (matches Fortran gpresp_induct)
     temp2 = Matrix{ComplexF64}(undef, mpert, mpert)
diff --git a/src/PerturbedEquilibrium/SingularCoupling.jl b/src/PerturbedEquilibrium/SingularCoupling.jl
index 4106be093..94f6aed96 100644
--- a/src/PerturbedEquilibrium/SingularCoupling.jl
+++ b/src/PerturbedEquilibrium/SingularCoupling.jl
@@ -209,7 +209,7 @@ end
         state::PerturbedEquilibriumState,
         equil::Equilibrium.PlasmaEquilibrium,
         ForceFreeStates_results::OdeState,
-        vac_data::VacuumData,
+        mthvac::Int,
         ffs_intr::ForceFreeStatesInternal,
         intr::PerturbedEquilibriumInternal,
         ctrl::PerturbedEquilibriumControl,
@@ -241,7 +241,7 @@ function compute_singular_coupling_metrics!(
     state::PerturbedEquilibriumState,
     equil::Equilibrium.PlasmaEquilibrium,
     ForceFreeStates_results::OdeState,
-    vac_data::VacuumData,
+    mthvac::Int,
     ffs_intr::ForceFreeStatesInternal,
     intr::PerturbedEquilibriumInternal,
     ctrl::PerturbedEquilibriumControl,
@@ -263,7 +263,7 @@ function compute_singular_coupling_metrics!(
 
     chi1 = 2π * equil.psio
     twopi = 2π
-    mtheta = vac_data.mthvac
+    mtheta = mthvac
     wall_settings = Vacuum.WallShapeSettings(; shape="nowall")
 
     # Phase 1: Collect all resonant (surface, n) pairs in psi order
@@ -348,9 +348,9 @@ function compute_singular_coupling_metrics!(
 
             # Compute Green's functions at this surface for this n (once per pair)
             vac_input = Vacuum.VacuumInput(equil, sing_surf.psifac, mtheta, 1, mlow:mhigh, [nn])
-            _, grri_raw, grre_raw, _, _ = Vacuum.compute_vacuum_response(vac_input, wall_settings)
-            grri = Matrix{ComplexF64}(grri_raw)
-            grre = Matrix{ComplexF64}(grre_raw)
+            vac = Vacuum.compute_vacuum_response(vac_input, wall_settings)
+            grri = Matrix{ComplexF64}(vac.grri)
+            grre = Matrix{ComplexF64}(vac.grre)
 
             # Get ν on the vacuum theta grid (same ν used in the vacuum Fourier basis computation)
             ν_vac = Vacuum.PlasmaGeometry(vac_input).ν
diff --git a/src/Vacuum/DataTypes.jl b/src/Vacuum/DataTypes.jl
index 06ee4a6a0..63dec0b9e 100644
--- a/src/Vacuum/DataTypes.jl
+++ b/src/Vacuum/DataTypes.jl
@@ -143,6 +143,43 @@ function expand_field_periods(inputs::VacuumInput)
     )
 end
 
+"""
+    VacuumResponse
+
+Output of [`compute_vacuum_response`](@ref): the vacuum energy matrix and the surface data the
+boundary-integral solve produces along the way.
+
+## Fields
+
+  - `wv::Matrix{ComplexF64}`: Vacuum energy matrix Wᵛ (`num_modes × num_modes`), block-diagonal in n for 2D
+  - `grri`, `grre::Matrix{ComplexF64}`: Interior/exterior Green's functions (`2·num_points × num_modes`), zeroed on every 3D path
+  - `plasma_pts`, `wall_pts::Matrix{Float64}`: Cartesian surface coordinates (`num_points × 3`)
+"""
+struct VacuumResponse
+    wv::Matrix{ComplexF64}
+    grri::Matrix{ComplexF64}
+    grre::Matrix{ComplexF64}
+    plasma_pts::Matrix{Float64}
+    wall_pts::Matrix{Float64}
+end
+
+"""
+    VacuumResponse(inputs::VacuumInput) -> VacuumResponse
+
+Allocate zeroed output arrays sized for `inputs` (full torus, so `nfp` field periods).
+"""
+function VacuumResponse(inputs::VacuumInput)
+    num_points = inputs.mtheta * inputs.nzeta * inputs.nfp
+    num_modes = length(inputs.m_modes) * length(inputs.n_modes)
+    return VacuumResponse(
+        zeros(ComplexF64, num_modes, num_modes),
+        zeros(ComplexF64, 2 * num_points, num_modes),
+        zeros(ComplexF64, 2 * num_points, num_modes),
+        zeros(num_points, 3),
+        zeros(num_points, 3)
+    )
+end
+
 """
     WallShapeSettings
 
diff --git a/src/Vacuum/Vacuum.jl b/src/Vacuum/Vacuum.jl
index b060f59fc..30af03b9a 100644
--- a/src/Vacuum/Vacuum.jl
+++ b/src/Vacuum/Vacuum.jl
@@ -18,7 +18,7 @@ include("Kernel2D.jl")
 include("Kernel3D.jl")
 include("Field.jl")
 
-export VacuumInput, WallShapeSettings
+export VacuumInput, VacuumResponse, WallShapeSettings
 export compute_vacuum_response, compute_vacuum_response!, compute_vacuum_field
 export extract_plasma_surface_at_psi
 export PlasmaGeometry
@@ -53,7 +53,7 @@ function _symmetrize_vacuum_energy!(wv::AbstractMatrix)
 end
 
 """
-    _compute_vacuum_response_2d!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings)
+    _compute_vacuum_response_2d!(vac_data::VacuumResponse, inputs::VacuumInput, wall_settings::WallShapeSettings)
 
 2D (axisymmetric) vacuum response calculation.
 
@@ -62,12 +62,14 @@ building the double-/single-layer operators, solving the exterior and interior s
 filling the corresponding diagonal block of the response matrix and the matching column block
 of the Green's functions.
 """
-@with_pool pool function _compute_vacuum_response_2d!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings)
+@with_pool pool function _compute_vacuum_response_2d!(vac_data::VacuumResponse, inputs::VacuumInput, wall_settings::WallShapeSettings)
 
     mpert = length(inputs.m_modes)
     num_points_surf = inputs.mtheta
 
-    vac_data.wv .= 0
+    fill!(vac_data.wv, 0)
+    fill!(vac_data.grri, 0)
+    fill!(vac_data.grre, 0)
 
     # Form the plasma and wall geometries
     plasma_surf = PlasmaGeometry(inputs)
@@ -147,7 +149,7 @@ of the Green's functions.
 end
 
 """
-    _compute_vacuum_response_3d!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings)
+    _compute_vacuum_response_3d!(vac_data::VacuumResponse, inputs::VacuumInput, wall_settings::WallShapeSettings)
 
 3D (`inputs.nzeta > 1`) vacuum response via block-circulant field-period reduction. For `nfp == 1`
 the block-circulant assembly and residue-class loop are skipped in favour of a more efficient
@@ -169,7 +171,7 @@ Only `wv` is produced currently; `grri`/`grre` are returned zeroed and are not y
 Extension point: per residue class, apply the per-period basis to `D̂ₖ⁻¹Ŝₖ` for the exterior columns and the
 interior variant `-D + 2I` for the interior columns, then scatter back into the `[2N × 2·num_modes]` arrays.
 """
-@with_pool pool function _compute_vacuum_response_3d!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings)
+@with_pool pool function _compute_vacuum_response_3d!(vac_data::VacuumResponse, inputs::VacuumInput, wall_settings::WallShapeSettings)
 
     (; mtheta, nzeta, nfp, m_modes, n_modes) = inputs
     fill!(vac_data.wv, 0)
@@ -257,57 +259,26 @@ interior variant `-D + 2I` for the interior columns, then scatter back into the
 end
 
 """
-    compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSettings)
-
-Allocate and return the vacuum response matrix and Green's functions for the given vacuum
-inputs. Thin allocating wrapper around the in-place [`compute_vacuum_response!`]: it sizes the
-output arrays for the full torus and forwards to the same 2D/3D workers. For performance-critical
-paths that already own preallocated storage (e.g. `ForceFreeStates.VacuumData`), prefer the
-in-place method to avoid extra heap allocations.
+    compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSettings) -> VacuumResponse
 
-# Returns
-
-  - `wv`: complex vacuum response matrix (`num_modes × num_modes`).
-  - `grri`, `grre`: interior/exterior Green's functions (zeroed on the 3D nowall path).
-  - `plasma_pts`, `wall_pts`: surface coordinate arrays.
+Compute the vacuum response for the given inputs. Allocating wrapper around
+[`compute_vacuum_response!`](@ref); pass a preallocated [`VacuumResponse`](@ref) to that method
+instead when reusing storage across calls.
 """
 function compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSettings)
-
-    num_points = inputs.mtheta * inputs.nzeta * inputs.nfp # mtheta for 2D
-    num_modes = length(inputs.m_modes) * length(inputs.n_modes)
-
-    vac = (
-        wv=zeros(ComplexF64, num_modes, num_modes),
-        grri=zeros(ComplexF64, 2 * num_points, num_modes),
-        grre=zeros(ComplexF64, 2 * num_points, num_modes),
-        plasma_pts=zeros(num_points, 3),
-        wall_pts=zeros(num_points, 3)
-    )
+    vac = VacuumResponse(inputs)
     compute_vacuum_response!(vac, inputs, wall_settings)
-
-    return vac.wv, vac.grri, vac.grre, vac.plasma_pts, vac.wall_pts
+    return vac
 end
 
 """
-    compute_vacuum_response!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings)
-
-In-place variant that computes the vacuum response and directly populates the arrays stored in
-`vac_data`. Dispatches on dimensionality only: 2D (`inputs.nzeta == 1`) routes to
-[`_compute_vacuum_response_2d!`], 3D to [`_compute_vacuum_response_3d!`].
-
-The `vac_data` argument is expected to provide the following writable fields with compatible
-sizes:
-
-  - `wv::AbstractMatrix{ComplexF64}`             – vacuum response matrix
-  - `grri::AbstractMatrix{ComplexF64}`           – interior Green's functions
-  - `grre::AbstractMatrix{ComplexF64}`           – exterior Green's functions
-  - `plasma_pts::AbstractMatrix{Float64}`        – plasma surface coordinates
-  - `wall_pts::AbstractMatrix{Float64}`          – wall surface coordinates
+    compute_vacuum_response!(vac_data::VacuumResponse, inputs::VacuumInput, wall_settings::WallShapeSettings)
 
-This is designed to work with `ForceFreeStates.VacuumData` but does not depend on its concrete
-type (duck-typed on field names only).
+In-place variant that populates the arrays of an existing [`VacuumResponse`](@ref). Dispatches on
+dimensionality only: 2D (`inputs.nzeta == 1`) routes to [`_compute_vacuum_response_2d!`], 3D to
+[`_compute_vacuum_response_3d!`].
 """
-function compute_vacuum_response!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings)
+function compute_vacuum_response!(vac_data::VacuumResponse, inputs::VacuumInput, wall_settings::WallShapeSettings)
     if inputs.nzeta == 1
         _compute_vacuum_response_2d!(vac_data, inputs, wall_settings)
     else
diff --git a/test/runtests_imas.jl b/test/runtests_imas.jl
index d2afd3091..774fdbfe8 100644
--- a/test/runtests_imas.jl
+++ b/test/runtests_imas.jl
@@ -115,7 +115,7 @@ using GeneralizedPerturbedEquilibrium.Equilibrium
         mock_et      = [0.4+0im, 0.7+0im, 1.1+0im]
         mock_n_idx   = [0, 0, 0]   # all belong to n=1 (j=0)
         mock_result  = (
-            vac_data = (et=mock_et, n_tor_idx=mock_n_idx),
+            free_energies = (et=mock_et, n_tor_idx=mock_n_idx),
             intr     = (numpert_total=3, npert=1, nlow=1),
         )
 
@@ -145,7 +145,7 @@ using GeneralizedPerturbedEquilibrium.Equilibrium
         mock_et    = [0.3+0im, 0.5+0im, 0.6+0im, 0.7+0im]
         mock_n_idx = [0, 1, 0, 1]
         mock_result = (
-            vac_data = (et=mock_et, n_tor_idx=mock_n_idx),
+            free_energies = (et=mock_et, n_tor_idx=mock_n_idx),
             intr     = (numpert_total=4, npert=2, nlow=1),
         )
 
@@ -169,7 +169,7 @@ using GeneralizedPerturbedEquilibrium.Equilibrium
         # Single n=1 run
         dd_single = IMASdd.dd()
         result_single = (
-            vac_data = (et=[0.3+0im, 0.6+0im], n_tor_idx=[0, 0]),
+            free_energies = (et=[0.3+0im, 0.6+0im], n_tor_idx=[0, 0]),
             intr     = (numpert_total=2, npert=1, nlow=1),
         )
         GeneralizedPerturbedEquilibrium.write_imas(dd_single, result_single)
@@ -177,7 +177,7 @@ using GeneralizedPerturbedEquilibrium.Equilibrium
         # Single n=2 run
         dd_single2 = IMASdd.dd()
         result_single2 = (
-            vac_data = (et=[0.5+0im, 0.7+0im], n_tor_idx=[0, 0]),
+            free_energies = (et=[0.5+0im, 0.7+0im], n_tor_idx=[0, 0]),
             intr     = (numpert_total=2, npert=1, nlow=2),
         )
         GeneralizedPerturbedEquilibrium.write_imas(dd_single2, result_single2)
@@ -185,7 +185,7 @@ using GeneralizedPerturbedEquilibrium.Equilibrium
         # Combined n=1,2 run (same eigenvalues, now interleaved)
         dd_multi = IMASdd.dd()
         result_multi = (
-            vac_data = (et=[0.3+0im, 0.5+0im, 0.6+0im, 0.7+0im], n_tor_idx=[0, 1, 0, 1]),
+            free_energies = (et=[0.3+0im, 0.5+0im, 0.6+0im, 0.7+0im], n_tor_idx=[0, 1, 0, 1]),
             intr     = (numpert_total=4, npert=2, nlow=1),
         )
         GeneralizedPerturbedEquilibrium.write_imas(dd_multi, result_multi)
diff --git a/test/runtests_parallel_integration.jl b/test/runtests_parallel_integration.jl
index b045cf028..7f6c57394 100644
--- a/test/runtests_parallel_integration.jl
+++ b/test/runtests_parallel_integration.jl
@@ -254,7 +254,7 @@ using TOML
             metric = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_metric(equil, intr.mpert)
             ffit = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_matrix(equil, intr, metric)
             odet, _, _, _ = GeneralizedPerturbedEquilibrium.ForceFreeStates.eulerlagrange_integration(ctrl, equil, ffit, intr)
-            vac = GeneralizedPerturbedEquilibrium.ForceFreeStates.free_run!(odet, ctrl, equil, ffit, intr)
+            vac = GeneralizedPerturbedEquilibrium.ForceFreeStates.free_run(odet, ctrl, equil, ffit, intr)
             return real(vac.et[1]), intr
         end
 
@@ -313,7 +313,7 @@ using TOML
             metric = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_metric(equil, intr.mpert)
             ffit = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_matrix(equil, intr, metric)
             odet, _, _, _ = GeneralizedPerturbedEquilibrium.ForceFreeStates.eulerlagrange_integration(ctrl, equil, ffit, intr)
-            vac = GeneralizedPerturbedEquilibrium.ForceFreeStates.free_run!(odet, ctrl, equil, ffit, intr)
+            vac = GeneralizedPerturbedEquilibrium.ForceFreeStates.free_run(odet, ctrl, equil, ffit, intr)
             return real(vac.et[1]), intr
         end
 
@@ -522,7 +522,7 @@ using TOML
         ffit = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_matrix(equil, intr, metric)
         odet, fm_propagators, fm_chunks, fm_S_left =
             GeneralizedPerturbedEquilibrium.ForceFreeStates.eulerlagrange_integration(ctrl, equil, ffit, intr)
-        vac = GeneralizedPerturbedEquilibrium.ForceFreeStates.free_run!(odet, ctrl, equil, ffit, intr)
+        vac = GeneralizedPerturbedEquilibrium.ForceFreeStates.free_run(odet, ctrl, equil, ffit, intr)
         GeneralizedPerturbedEquilibrium.ForceFreeStates.compute_delta_prime_matrix!(
             intr, fm_propagators, fm_chunks;
             wv=vac.wv, psio=equil.psio,
diff --git a/test/runtests_riccati.jl b/test/runtests_riccati.jl
index 84194d948..57b345702 100644
--- a/test/runtests_riccati.jl
+++ b/test/runtests_riccati.jl
@@ -118,14 +118,14 @@ end
     # (which overwrites intr_ric.sing[s].delta_prime)
     delta_prime_inline = [copy(intr_ric.sing[s].delta_prime) for s in 1:intr_ric.msing]
 
-    vac_ric = FFS.free_run!(odet_ric, ctrl, equil, ffit, intr_ric)
+    vac_ric = FFS.free_run(odet_ric, ctrl, equil, ffit, intr_ric)
     et_ric = real(vac_ric.et[1])
 
     # Standard integration (needed only for energy comparison).  eulerlagrange_integration
     # returns (odet, propagators, chunks, S_at_surface_left); only odet is used here.
     intr_std = make_solovev_intr(inputs, ctrl, equil, ex)
     odet_std, _, _, _ = FFS.eulerlagrange_integration(ctrl, equil, ffit, intr_std)
-    vac_std = FFS.free_run!(odet_std, ctrl, equil, ffit, intr_std)
+    vac_std = FFS.free_run(odet_std, ctrl, equil, ffit, intr_std)
     et_std = real(vac_std.et[1])
 
     # ─────────────────────────────────────────────────────────────────────────
@@ -215,7 +215,7 @@ end
         # result must be bit-for-bit identical (not just approximately equal).
         #
         # Note: this call overwrites intr_ric.sing[s].delta_prime; delta_prime_inline was
-        # saved before free_run! above so it holds the original inline values.
+        # saved before free_run above so it holds the original inline values.
         #
         # See benchmarks/benchmark_delta_prime_methods.jl for the extended version.
         FFS.compute_delta_prime_from_ca!(odet_ric, intr_ric, equil)
diff --git a/test/runtests_vacuum.jl b/test/runtests_vacuum.jl
index 7e8b51b46..630c60436 100644
--- a/test/runtests_vacuum.jl
+++ b/test/runtests_vacuum.jl
@@ -368,7 +368,7 @@
             @testset "nowall" begin
                 inputs = _make_inputs()
                 wall_settings = WallShapeSettings(shape="nowall")
-                wv, grri, grre, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings)
+                (; wv, grri, grre, plasma_pts, wall_pts) = compute_vacuum_response(inputs, wall_settings)
 
                 numpoints = inputs.mtheta * inputs.nzeta
                 num_modes = length(inputs.m_modes) * length(inputs.n_modes)
@@ -391,7 +391,7 @@
             @testset "conformal wall" begin
                 inputs = _make_inputs()
                 wall_settings = WallShapeSettings(shape="conformal", a=0.5)
-                wv, grri, grre, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings)
+                (; wv, grri, grre, plasma_pts, wall_pts) = compute_vacuum_response(inputs, wall_settings)
 
                 numpoints = inputs.mtheta * inputs.nzeta
                 num_modes = length(inputs.m_modes) * length(inputs.n_modes)
@@ -408,7 +408,7 @@
             @testset "edge: single poloidal mode mpert=1" begin
                 inputs = _make_inputs(m_modes=[1], n_modes=[1])
                 wall_settings = WallShapeSettings(shape="nowall")
-                wv, grri, grre, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings)
+                (; wv, grri, grre, plasma_pts, wall_pts) = compute_vacuum_response(inputs, wall_settings)
                 @test size(wv) == (1, 1)
                 @test all(isfinite, wv)
                 @test size(grri, 2) == 1
@@ -418,7 +418,7 @@
                 # Keep mtheta_eq=17 so boundary has enough points for periodic spline
                 inputs = _make_inputs(mtheta=16, mtheta_eq=17)
                 wall_settings = WallShapeSettings(shape="nowall")
-                wv, grri, grre, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings)
+                (; wv, grri, grre, plasma_pts, wall_pts) = compute_vacuum_response(inputs, wall_settings)
                 @test size(wv) == (2, 2)
                 @test size(grri) == (32, 2)  # 2*mtheta, num_modes=2
                 @test size(plasma_pts) == (16, 3)
@@ -426,24 +426,41 @@
 
             @testset "in-place compute_vacuum_response! matches wrapper" begin
                 # The allocating wrapper is a thin caller of the in-place routine; verify the
-                # in-place entry populates caller-owned duck-typed (NamedTuple) storage identically.
+                # in-place entry populates caller-owned storage identically.
                 for wall_settings in (WallShapeSettings(shape="nowall"), WallShapeSettings(shape="conformal", a=0.5))
                     inputs = _make_inputs()
-                    wv, grri, grre, pp, wp = compute_vacuum_response(inputs, wall_settings)
+                    ref = compute_vacuum_response(inputs, wall_settings)
 
-                    numpoints = inputs.mtheta * inputs.nzeta
-                    num_modes = length(inputs.m_modes) * length(inputs.n_modes)
-                    vac = (wv=zeros(ComplexF64, num_modes, num_modes), grri=zeros(ComplexF64, 2 * numpoints, num_modes), grre=zeros(ComplexF64, 2 * numpoints, num_modes),
-                        plasma_pts=zeros(numpoints, 3), wall_pts=zeros(numpoints, 3))
+                    vac = VacuumResponse(inputs)
                     compute_vacuum_response!(vac, inputs, wall_settings)
 
-                    @test vac.wv ≈ wv
-                    @test vac.grri ≈ grri
-                    @test vac.grre ≈ grre
-                    @test vac.plasma_pts ≈ pp
-                    @test vac.wall_pts ≈ wp
+                    @test vac.wv ≈ ref.wv
+                    @test vac.grri ≈ ref.grri
+                    @test vac.grre ≈ ref.grre
+                    @test vac.plasma_pts ≈ ref.plasma_pts
+                    @test vac.wall_pts ≈ ref.wall_pts
                 end
             end
+
+            @testset "in-place compute_vacuum_response! clears a reused buffer" begin
+                # The nowall path writes only the plasma rows of grri/grre, so a buffer left
+                # over from a wall run must not leak its wall rows into the next result.
+                inputs = _make_inputs()
+                wall_rows = (inputs.mtheta+1):(2*inputs.mtheta)
+
+                vac = VacuumResponse(inputs)
+                compute_vacuum_response!(vac, inputs, WallShapeSettings(; shape="conformal", a=0.5))
+                @test any(!iszero, view(vac.grre, wall_rows, :))
+
+                compute_vacuum_response!(vac, inputs, WallShapeSettings(; shape="nowall"))
+                fresh = compute_vacuum_response(inputs, WallShapeSettings(; shape="nowall"))
+
+                @test all(iszero, view(vac.grri, wall_rows, :))
+                @test all(iszero, view(vac.grre, wall_rows, :))
+                @test vac.wv ≈ fresh.wv
+                @test vac.grri ≈ fresh.grri
+                @test vac.grre ≈ fresh.grre
+            end
         end
 
         @testset "extract_plasma_surface_at_psi" begin
@@ -581,7 +598,7 @@
         @testset "compute_vacuum_response 3D nowall" begin
             inputs = _make_3d_inputs(mtheta=32, nzeta=32, mtheta_eq=17)
             wall_settings = WallShapeSettings(shape="nowall")
-            wv, grri, grre, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings)
+            (; wv, grri, grre, plasma_pts, wall_pts) = compute_vacuum_response(inputs, wall_settings)
 
             numpoints = inputs.mtheta * inputs.nzeta
             num_modes = length(inputs.m_modes) * length(inputs.n_modes)
@@ -604,7 +621,7 @@
         @testset "compute_vacuum_response 3D nonaxisymmetric boundary" begin
             inputs = _make_3d_nonaxis_inputs(mtheta=24, nzeta=24, mtheta_in=12, nzeta_in=12, mpert=2, nlow=0, npert=2)
             wall_settings = WallShapeSettings(shape="nowall")
-            wv, grri, grre, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings)
+            (; wv, grri, grre, plasma_pts, wall_pts) = compute_vacuum_response(inputs, wall_settings)
 
             numpoints = inputs.mtheta * inputs.nzeta
             num_modes = length(inputs.m_modes) * length(inputs.n_modes)
@@ -625,7 +642,7 @@
         @testset "compute_vacuum_response 3D conformal wall" begin
             inputs = _make_3d_inputs(mtheta=32, nzeta=32, mtheta_eq=17)
             wall_settings = WallShapeSettings(shape="conformal", a=0.3)
-            wv, grri, grre, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings)
+            (; wv, grri, grre, plasma_pts, wall_pts) = compute_vacuum_response(inputs, wall_settings)
 
             numpoints = inputs.mtheta * inputs.nzeta
             num_modes = length(inputs.m_modes) * length(inputs.n_modes)
@@ -676,13 +693,14 @@
             wall_settings = WallShapeSettings(shape="nowall")
 
             # Reduced (block-circulant) path
-            wv_red, _, _, plasma_pts_red, _ = compute_vacuum_response(inputs_red, wall_settings)
+            vac_red = compute_vacuum_response(inputs_red, wall_settings)
+            wv_red, plasma_pts_red = vac_red.wv, vac_red.plasma_pts
 
             # Full-torus reference: pre-expand so nfp=1 forces the dense path
             inputs_full = GeneralizedPerturbedEquilibrium.Vacuum.expand_field_periods(inputs_red)
             @test inputs_full.nfp == 1
             @test inputs_full.nzeta == nzeta_p * nfp
-            wv_full, _, _, _, _ = compute_vacuum_response(inputs_full, wall_settings)
+            wv_full = compute_vacuum_response(inputs_full, wall_settings).wv
 
             num_modes = length(m_modes) * length(n_modes)
             @test size(wv_red) == (num_modes, num_modes)