diff --git a/regression-harness/cases/diiid_n1_riccati.toml b/regression-harness/cases/diiid_n1_riccati.toml
index f41d9ed0a..811ed13c0 100644
--- a/regression-harness/cases/diiid_n1_riccati.toml
+++ b/regression-harness/cases/diiid_n1_riccati.toml
@@ -10,13 +10,25 @@ name = "diiid_n1_riccati"
 description = "DIII-D-like equilibrium, n=1, Riccati integrator Δ' matrix"
 example_dir = "examples/DIIID-like_riccati_deltaprime_example"
 
+# Which formalism produced the Δ′ below. The Δ′ matrix is written to one canonical
+# SingularSurfaces/ path whichever integrator computed it, so the numbers alone cannot say
+# whether a gold pinned from the BVP is being compared against a Galerkin run. Pinning the
+# deck's declared choice makes a formalism swap a loud failure instead of a silent re-gold.
+[quantities.integrator]
+h5path = "Input/gpec_toml_raw"
+type = "token"
+extract = "toml_key:ForceFreeStates.integrator"
+label = "integrator (Δ′ provenance)"
+noise_threshold = 0
+order = 1
+
 # Tearing stability Δ' — the msing×msing BVP matrix diagonal, and the raw side-major
 # outer-region matrix it is projected from.
 [quantities.delta_prime]
 h5path = "SingularSurfaces/Delta_prime_matrix"
 type = "complex_matrix"
 extract = "diagonal_complex"
-label = "delta prime (BVP diagonal)"
+label = "Δ′ STRIDE-like BVP (PEST3) diagonal"
 noise_threshold = 1e-8
 order = 10
 
@@ -24,7 +36,7 @@ order = 10
 h5path = "SingularSurfaces/Delta_prime_raw"
 type = "complex_matrix"
 extract = "all_complex"
-label = "delta prime (raw side-major)"
+label = "Δ′ STRIDE-like BVP (PEST3) raw side-major"
 noise_threshold = 1e-8
 order = 11
 
diff --git a/regression-harness/cases/diiid_slayer_n1.toml b/regression-harness/cases/diiid_slayer_n1.toml
index 540ca7e50..69bd7b656 100644
--- a/regression-harness/cases/diiid_slayer_n1.toml
+++ b/regression-harness/cases/diiid_slayer_n1.toml
@@ -9,6 +9,17 @@ description = "DIII-D-like H-mode equilibrium, n=1, SLAYER tearing-mode analysis
 example_dir = "examples/DIIID-like_SLAYER_example"
 
 # Per-surface SLAYER layer parameters (geometry + dimensionless)
+# Which formalism produced the Δ′ this case's layer physics consumes. The Δ′ matrix lands on
+# one canonical SingularSurfaces/ path whichever integrator computed it, so pinning the deck's
+# declared choice is what keeps a BVP-derived gold from being compared against a Galerkin run.
+[quantities.integrator]
+h5path = "Input/gpec_toml_raw"
+type = "token"
+extract = "toml_key:ForceFreeStates.integrator"
+label = "integrator (Δ′ provenance)"
+noise_threshold = 0
+order = 1
+
 [quantities.slayer_ising]
 h5path = "Tearing/PerSurface/rational_index"
 type = "real_vector"
@@ -132,6 +143,52 @@ label = "SLAYER no_root flags [2/1,3/1,4/1]"
 noise_threshold = 0
 order = 34
 
+# Outer-region Δ' consumed by the SLAYER matching — the same STRIDE-like BVP (PEST3) matrix
+# the riccati case pins at source; tracked here too so a change in what SLAYER actually
+# ingests is caught even if the source case is not run.
+[quantities.slayer_delta_prime_diag]
+h5path = "Tearing/PerSurface/Delta_prime_matrix"
+type = "complex_matrix"
+extract = "diagonal_complex"
+label = "Δ′ STRIDE-like BVP (PEST3) diagonal, as ingested by SLAYER"
+noise_threshold = 1e-8
+order = 40
+
+# Inner-layer solve outputs per surface (the resistive-layer Riccati solution itself, not the
+# outer matching). The adaptive-scan Δ(Q) sample arrays are deliberately NOT pinned: AMR sample
+# locations move under any refinement change, so they have no stable identity to pin.
+[quantities.slayer_delta_s]
+h5path = "Tearing/LayerWidths/delta_s"
+type = "complex_vector"
+extract = "all_complex"
+label = "SLAYER inner-layer thickness δ_s (complex, Riccati)"
+noise_threshold = 1e-10
+order = 41
+
+[quantities.slayer_delta_s_over_d_beta]
+h5path = "Tearing/LayerWidths/delta_s_over_d_beta"
+type = "complex_vector"
+extract = "all_complex"
+label = "SLAYER inner-layer δ_s/d_β (complex, dimensionless)"
+noise_threshold = 1e-10
+order = 42
+
+[quantities.slayer_delta_s_abs]
+h5path = "Tearing/LayerWidths/delta_s_abs"
+type = "real_vector"
+extract = "all_real"
+label = "SLAYER inner-layer physical thickness |δ_s| [m]"
+noise_threshold = 1e-12
+order = 43
+
+[quantities.slayer_d_beta]
+h5path = "Tearing/LayerWidths/d_beta"
+type = "real_vector"
+extract = "all_real"
+label = "SLAYER β-weighted ion drift scale d_β [m]"
+noise_threshold = 1e-12
+order = 44
+
 # Settings (catches accidental config drift)
 [quantities.slayer_enabled]
 h5path = "Tearing/enabled"
diff --git a/regression-harness/cases/gal_resistive_diiid.toml b/regression-harness/cases/gal_resistive_diiid.toml
index b13388d8a..bc0079848 100644
--- a/regression-harness/cases/gal_resistive_diiid.toml
+++ b/regression-harness/cases/gal_resistive_diiid.toml
@@ -9,6 +9,18 @@ name = "gal_resistive_diiid"
 description = "DIII-D-like, n=1, RDCON outer-region Galerkin Δ′ with rpec coil columns (delta_coil)"
 example_dir = "examples/DIIID-like_gal_resistive_example"
 
+# Which formalism produced the Galerkin Δ′ below. The Δ′ matrix is written to one canonical
+# SingularSurfaces/ path whichever integrator computed it, so the numbers alone cannot say whether
+# a gold pinned from the Galerkin solve is being compared against a BVP run. Pinning the deck's
+# declared choice makes a formalism swap a loud failure instead of a silent re-gold.
+[quantities.integrator]
+h5path = "Input/gpec_toml_raw"
+type = "token"
+extract = "toml_key:ForceFreeStates.integrator"
+label = "integrator (Δ′ provenance)"
+noise_threshold = 0
+order = 1
+
 # Number of resonant surfaces in the gal domain
 [quantities.gal_msing]
 h5path = "ForceFreeStates/Solutions/GalerkinIntegration/rational_count"
@@ -31,7 +43,7 @@ order = 11
 h5path = "SingularSurfaces/Delta_prime_matrix"
 type = "complex_matrix"
 extract = "diagonal_complex"
-label = "gal PEST3 Δ diagonal"
+label = "Δ′ Galerkin RDCON (PEST3) diagonal"
 noise_threshold = 1e-6
 order = 20
 
@@ -40,7 +52,7 @@ order = 20
 h5path = "SingularSurfaces/Delta_prime_raw"
 type = "complex_matrix"
 extract = "norm"
-label = "||gal Δ′ matrix||"
+label = "||Δ′ Galerkin RDCON (PEST3) matrix||"
 noise_threshold = 1e-6
 order = 21
 
diff --git a/regression-harness/cases/gal_resistive_pe.toml b/regression-harness/cases/gal_resistive_pe.toml
index 9a04b30ac..f943e8d49 100644
--- a/regression-harness/cases/gal_resistive_pe.toml
+++ b/regression-harness/cases/gal_resistive_pe.toml
@@ -13,6 +13,18 @@ example_dir = "examples/DIIID-like_gal_resistive_pe_example"
 # FINITE, per-rational-surface singular-coupling scalars (the physics-meaningful driven quantities) — not
 # any whole-profile ξ/ξ′ norm (those diverge at the rationals and are ill-posed as benchmarks).
 
+# Which formalism produced the Galerkin Δ′ below. The Δ′ matrix is written to one canonical
+# SingularSurfaces/ path whichever integrator computed it, so the numbers alone cannot say whether
+# a gold pinned from the Galerkin solve is being compared against a BVP run. Pinning the deck's
+# declared choice makes a formalism swap a loud failure instead of a silent re-gold.
+[quantities.integrator]
+h5path = "Input/gpec_toml_raw"
+type = "token"
+extract = "toml_key:ForceFreeStates.integrator"
+label = "integrator (Δ′ provenance)"
+noise_threshold = 0
+order = 1
+
 # Surface locations / count (sanity)
 [quantities.pe_rational_q]
 h5path = "PerturbedEquilibrium/SingularCoupling/rational_q"
@@ -66,12 +78,14 @@ label = "penetrated b field"
 noise_threshold = 1e-8
 order = 32
 
-# Per-surface Δ' from the PE singular-coupling analysis
+# Per-surface driven-response Δ' from the PE singular-coupling jump (∂b^ψ/∂ψ across the
+# rational surface over 2π·χ₁, forced by the applied flux). This is a driven-response
+# diagnostic, NOT the tearing-stability Δ' of the STRIDE-like BVP or Galerkin RDCON solvers.
 [quantities.pe_delta_prime]
 h5path = "PerturbedEquilibrium/SingularCoupling/Delta_prime"
 type = "complex_vector"
 extract = "all_complex"
-label = "PE Δ' per surface"
+label = "PE driven-response Δ′ (SingularCoupling jump)"
 noise_threshold = 1e-8
 order = 33
 
diff --git a/regression-harness/cases/slayer_dels_fitzpatrick.toml b/regression-harness/cases/slayer_dels_fitzpatrick.toml
new file mode 100644
index 000000000..cc3c226bf
--- /dev/null
+++ b/regression-harness/cases/slayer_dels_fitzpatrick.toml
@@ -0,0 +1,68 @@
+# Regression case: GPEC's del_s Riccati layer-width solver against Fitzpatrick, "Tearing Mode
+# Dynamics in Tokamak Plasmas" (IOP 2023), figures 6.2 and 6.3 (p. 6-10).
+#
+# This is the suite's only case pinned to an EXTERNAL published reference rather than to GPEC's own
+# previous output: a failure means we no longer reproduce the textbook. Verified at eight points
+# against the published contour bands, e.g. (Q,P) = (4.0, 0.5) -> Re 2.008 against the 1.9-2.1
+# contours, (0.5, 4.0) -> Re 1.537 against the 1.5 contour. The magnitudes also match the book's
+# table 6.2 reactor estimates (|delta_s| 2.95-4.69 mm) at DIII-D-like parameters.
+#
+# Pure ODE: no equilibrium, no geqdsk, no kinetic file. Sub-second, so it is fast-tier eligible, and
+# it carries no equilibrium dependence -- its golden values survive equilibrium-side re-baselining.
+# Complements slayer_delta_probe, which pins the dispersion solve Delta(Q); this pins the width
+# solve delta_s/d_beta over (Q_hat, P_hat).
+[case]
+name = "slayer_dels_fitzpatrick"
+description = "del_s Riccati layer width vs Fitzpatrick 2023 figs 6.2/6.3 over (Q_hat, P_hat)"
+kind = "computed"
+
+# The scan axes, pinned so a silently reshaped grid is a failure rather than a re-gold.
+[quantities.q_hat]
+h5path = "fitzpatrick/Q_hat"
+type = "real_vector"
+extract = "all_real"
+label = "grid Q_hat (= Q_*/D^4)"
+noise_threshold = 0
+order = 10
+
+[quantities.p_hat]
+h5path = "fitzpatrick/P_hat"
+type = "real_vector"
+extract = "all_real"
+label = "grid P_hat (= P/D^6)"
+noise_threshold = 0
+order = 11
+
+# tau is not stated in the book's figure captions; tau = 1 (table 5.1, low-field reactor) is the
+# assumption this reproduction rests on, so it is pinned rather than left implicit.
+[quantities.tau]
+h5path = "fitzpatrick/tau"
+type = "real_scalar"
+extract = "value"
+label = "tau (assumed; not given in the figure captions)"
+noise_threshold = 0
+order = 12
+
+[quantities.dels_db_re]
+h5path = "fitzpatrick/dels_db_re"
+type = "real_vector"
+extract = "all_real"
+label = "Re(delta_s/d_beta) — Fitzpatrick fig 6.2"
+noise_threshold = 1e-10
+order = 20
+
+[quantities.dels_db_im]
+h5path = "fitzpatrick/dels_db_im"
+type = "real_vector"
+extract = "all_real"
+label = "Im(delta_s/d_beta) — Fitzpatrick fig 6.3"
+noise_threshold = 1e-10
+order = 21
+
+[quantities.runtime]
+h5path = ""
+type = "runtime"
+extract = "value"
+label = "Runtime (s)"
+noise_threshold = 0.0
+order = 999
diff --git a/regression-harness/cases/slayer_delta_probe.toml b/regression-harness/cases/slayer_delta_probe.toml
new file mode 100644
index 000000000..41a7d8798
--- /dev/null
+++ b/regression-harness/cases/slayer_delta_probe.toml
@@ -0,0 +1,50 @@
+# Regression case: fixed-Q probe of the SLAYER inner-layer dispersion Δ(Q).
+# Evaluates Δ(Q) on a fixed 4×4 grid over Re(Q), Im(Q) ∈ [-10, 10] using the DIII-D-like
+# SLAYER deck's own 2/1-surface layer parameters, quoted verbatim in the runner's computed
+# script. Fixed Q values are refinement-stable, unlike the AMR scan samples in gpec.h5, so
+# this pins the dispersion curve itself; the parameter chain producing those layer inputs
+# is pinned separately by diiid_slayer_n1. Seconds to run — fast-tier eligible.
+[case]
+name = "slayer_delta_probe"
+description = "SLAYER inner-layer Δ(Q) on a fixed 4×4 Q grid, DIII-D-like 2/1 layer parameters"
+kind = "computed"
+
+[quantities.q_re]
+h5path = "slayer_probe/Q_re"
+type = "real_vector"
+extract = "all_real"
+label = "probe grid Re(Q)"
+noise_threshold = 0
+order = 10
+
+[quantities.q_im]
+h5path = "slayer_probe/Q_im"
+type = "real_vector"
+extract = "all_real"
+label = "probe grid Im(Q)"
+noise_threshold = 0
+order = 11
+
+[quantities.delta_re]
+h5path = "slayer_probe/Delta_re"
+type = "real_vector"
+extract = "all_real"
+label = "SLAYER inner-layer Re Δ(Q) on the fixed grid"
+noise_threshold = 1e-10
+order = 20
+
+[quantities.delta_im]
+h5path = "slayer_probe/Delta_im"
+type = "real_vector"
+extract = "all_real"
+label = "SLAYER inner-layer Im Δ(Q) on the fixed grid"
+noise_threshold = 1e-10
+order = 21
+
+[quantities.runtime]
+h5path = ""
+type = "runtime"
+extract = "value"
+label = "Runtime (s)"
+noise_threshold = 0.0
+order = 999
diff --git a/regression-harness/golden/slayer_dels_fitzpatrick.toml b/regression-harness/golden/slayer_dels_fitzpatrick.toml
new file mode 100644
index 000000000..ddb9ce848
--- /dev/null
+++ b/regression-harness/golden/slayer_dels_fitzpatrick.toml
@@ -0,0 +1,59 @@
+[meta]
+blas_threads = 2
+case = "slayer_dels_fitzpatrick"
+commit = "af3e4492"
+generated_at = "2026-08-18"
+golden_version = 1
+julia_version = "1.11.6"
+manifest_sha = "7e5c34ad1722fa6addf94f7387b972fa8ac17ff86166793819b158496d495132"
+nthreads = 2
+os_arch = "arm64-apple-darwin24.0.0"
+reason = "Initial pin: validated against Fitzpatrick 2023 figs 6.2/6.3 at eight contour points"
+
+[values.dels_db_im]
+atol = 0.0
+class = "physics_converged"
+rtol = 1.0e-6
+tolerance_basis = "class-default (provisional)"
+value_text = "[-0.03903669609227284,-0.7097932061530505,-1.0865009301031268,-1.5690999603536648,-2.2242340056733148,-0.011118192607611124,-0.35424388489008746,-0.776484631542084,-1.416939644030615,-2.2250603072320665,-0.008402000357094586,-0.23928032159546678,-0.5666784684880879,-1.241924594739638,-2.1750391797697133,-0.006297613829213364,-0.1657149360631757,-0.3716832232708168,-0.9092532037784165,-1.9732264855758332,-0.004682490543220957,-0.11916296976842101,-0.2502562086693129,-0.578522722631012,-1.4658508790642786]"
+value_type = "json_array"
+
+[values.dels_db_re]
+atol = 0.0
+class = "physics_converged"
+rtol = 1.0e-6
+tolerance_basis = "class-default (provisional)"
+value_text = "[0.46252746222290736,0.7691732891864712,1.0821020639186198,1.5537655127977916,2.2130417549691974,0.9838377906344368,0.9426353501706125,1.017567466663986,1.35072260354625,2.0080394342246684,1.1484429925403958,1.1069495858506195,1.0754748077370269,1.2328881423647,1.8273926389542299,1.336215433387759,1.3118923623802783,1.259193827055807,1.2120466712964955,1.540405108811276,1.5494641775774003,1.5373751336175125,1.5040680534472435,1.4084024975856595,1.350103015944749]"
+value_type = "json_array"
+
+[values.p_hat]
+atol = 0.0
+class = "physics_converged"
+rtol = 1.0e-6
+tolerance_basis = "class-default (provisional)"
+value_text = "[0.02,0.02,0.02,0.02,0.02,0.5,0.5,0.5,0.5,0.5,1.0,1.0,1.0,1.0,1.0,2.0,2.0,2.0,2.0,2.0,4.0,4.0,4.0,4.0,4.0]"
+value_type = "json_array"
+
+[values.q_hat]
+atol = 0.0
+class = "physics_converged"
+rtol = 1.0e-6
+tolerance_basis = "class-default (provisional)"
+value_text = "[0.02,0.5,1.0,2.0,4.0,0.02,0.5,1.0,2.0,4.0,0.02,0.5,1.0,2.0,4.0,0.02,0.5,1.0,2.0,4.0,0.02,0.5,1.0,2.0,4.0]"
+value_type = "json_array"
+
+[values.runtime]
+atol = +inf
+class = "diagnostic"
+rtol = +inf
+tolerance_basis = "class-default (provisional)"
+value = 1.2057960033416748
+value_type = "real"
+
+[values.tau]
+atol = 0.0
+class = "physics_converged"
+rtol = 1.0e-6
+tolerance_basis = "class-default (provisional)"
+value = 1.0
+value_type = "real"
diff --git a/regression-harness/regress.jl b/regression-harness/regress.jl
index 651a3df8d..371a69c7f 100644
--- a/regression-harness/regress.jl
+++ b/regression-harness/regress.jl
@@ -6,6 +6,7 @@ const HARNESS_DIR = @__DIR__
 const REPO_ROOT = abspath(joinpath(HARNESS_DIR, ".."))
 const DEFAULT_DB_PATH = joinpath(HARNESS_DIR, ".regress_cache.sqlite")
 const CASES_DIR = joinpath(HARNESS_DIR, "cases")
+const GOLDEN_DIR = joinpath(HARNESS_DIR, "golden")
 
 include("src/types.jl")
 include("src/env.jl")
@@ -15,6 +16,7 @@ include("src/utils.jl")
 include("src/extractor.jl")
 include("src/runner.jl")
 include("src/reporter.jl")
+include("src/golden.jl")
 
 function parse_args(args)
     cases = String[]
@@ -30,6 +32,9 @@ function parse_args(args)
     no_pin_manifest = false
     allow_env_mismatch = false
     fail_on_change = false
+    check = false
+    update_golden = false
+    reason = nothing
     help = false
 
     i = 1
@@ -59,6 +64,15 @@ function parse_args(args)
         elseif arg == "--fail-on-change"
             fail_on_change = true
             i += 1
+        elseif arg == "--check"
+            check = true
+            i += 1
+        elseif arg == "--update-golden"
+            update_golden = true
+            i += 1
+        elseif arg == "--reason" && i < length(args)
+            reason = args[i+1]
+            i += 2
         elseif arg == "--cases" && i < length(args)
             cases = split(args[i+1], ",") |> collect .|> strip
             i += 2
@@ -83,7 +97,7 @@ function parse_args(args)
     end
 
     return CLIOptions(cases, refs, ref_range, force, list_cases, show_qty, show_case, db_path, verbose,
-        no_instantiate, no_pin_manifest, allow_env_mismatch, fail_on_change, help)
+        no_instantiate, no_pin_manifest, allow_env_mismatch, fail_on_change, check, update_golden, reason, help)
 end
 
 const HELP_TEXT = """
@@ -117,15 +131,26 @@ Options:
                            re-running them
     --fail-on-change       Exit non-zero if any tracked quantity changed (for CI use; a failed
                            run always exits non-zero regardless)
+    --check                Run the working tree and compare against the committed golden values
+                           (regression-harness/golden/<case>.toml). This is the mode CI gates on.
+    --update-golden        Regenerate the golden values from a fresh run. Requires --reason.
+    --reason "..."         Why the goldens changed. Mandatory with --update-golden: a golden
+                           change is a claim about physics that a reviewer has to evaluate.
     --help                 Print this help message
 
 Environment:
     GPEC_REGRESS_THREADS   Threads for GPEC subprocesses (default "auto" = all cores);
                            the actual count is recorded in each run's env fingerprint
 
+Golden values:
+    A golden file pins each quantity together with the tolerance it must be reproduced within
+    and the evidence behind that tolerance. Tolerances are derived from measured convergence
+    plateau drift and cross-platform spread — never widened to make a check pass. A failing
+    --check is fixed by explaining the physics or fixing the regression, not by loosening rtol.
+
 Exit status:
     0  all runs completed (and, with --fail-on-change, nothing changed)
-    1  a run failed, or a quantity changed under --fail-on-change
+    1  a run failed, a quantity changed under --fail-on-change, or a golden check failed
 
 Examples:
     # Compare two refs
@@ -150,6 +175,13 @@ Examples:
     # Compare current working tree against develop
     julia --project=regression-harness regression-harness/regress.jl \\
         --cases solovev_n1 --refs develop,local
+
+    # Gate the working tree against the committed golden values
+    julia --project=regression-harness regression-harness/regress.jl --check --cases solovev_n1
+
+    # Regenerate goldens after an intended physics change
+    julia --project=regression-harness regression-harness/regress.jl \\
+        --update-golden --cases solovev_n1 --reason "two-point b1 fix changes resonant field"
 """
 
 function main(args=ARGS)
@@ -189,14 +221,29 @@ function main(args=ARGS)
             [all_cases[name] for name in opts.cases if haskey(all_cases, name) || error("Unknown case: $name")]
         end
 
-        # Resolve refs
+        golden_mode = opts.check || opts.update_golden
+        if opts.check && opts.update_golden
+            error("--check and --update-golden are mutually exclusive: one asserts the goldens, the other rewrites them")
+        end
+        if opts.update_golden && (opts.reason === nothing || isempty(strip(opts.reason)))
+            error("--update-golden requires --reason \"...\": a golden change is a claim about physics, " *
+                  "and the reviewer needs to know what changed and why")
+        end
+
+        # Resolve refs. Golden modes judge the working tree, so they default to it rather than
+        # requiring the caller to spell out --refs local.
         resolved_refs = if opts.ref_range !== nothing
             expand_ref_range(opts.ref_range, REPO_ROOT)
         elseif !isempty(opts.refs)
             [resolve_ref(ref, REPO_ROOT) for ref in opts.refs]
+        elseif golden_mode
+            [resolve_ref(LOCAL_REF, REPO_ROOT)]
         else
             error("Must specify --refs or --ref-range")
         end
+        if golden_mode && length(resolved_refs) != 1
+            error("Golden modes compare one ref against the committed values; got $(length(resolved_refs)) refs")
+        end
 
         if isempty(resolved_refs)
             error("No commits resolved from the given refs")
@@ -220,6 +267,10 @@ function main(args=ARGS)
 
         n_failed = 0
         n_changed = 0
+        n_golden_fail = 0
+        n_untracked = 0
+        n_checked = 0
+        n_no_golden = 0
 
         # Run every case against each ref, grouped by ref so cases sharing a commit share
         # one worktree (and one Pkg.instantiate/precompile) instead of paying for it per case.
@@ -236,20 +287,46 @@ function main(args=ARGS)
             println("Case: $(case_spec.name) — $(case_spec.description)")
             println("="^64)
 
-            summary = if length(resolved_refs) == 2
-                report_two_ref_comparison(db, case_spec,
-                    resolved_refs[1], resolved_refs[2])
+            if opts.update_golden
+                update_golden_from_run(db, case_spec, resolved_refs[1].commit_hash,
+                    String(opts.reason), REPO_ROOT)
+            elseif opts.check
+                summary = report_golden_check(db, case_spec, resolved_refs[1].commit_hash)
+                n_golden_fail += summary.n_fail
+                n_untracked += summary.n_untracked
+                has_golden(case_spec.name) ? (n_checked += 1) : (n_no_golden += 1)
             else
-                report_multi_ref(db, case_spec, resolved_refs)
+                summary = if length(resolved_refs) == 2
+                    report_two_ref_comparison(db, case_spec,
+                        resolved_refs[1], resolved_refs[2])
+                else
+                    report_multi_ref(db, case_spec, resolved_refs)
+                end
+                n_failed += summary.n_failed
+                n_changed += summary.n_changed
             end
-            n_failed += summary.n_failed
-            n_changed += summary.n_changed
         end
 
         if n_failed > 0
             @error "$n_failed run(s) failed — see the reports above"
             exit(1)
         end
+        if n_golden_fail > 0
+            @error "$n_golden_fail quantity/quantities are outside their golden tolerance"
+            exit(1)
+        end
+        if opts.check && n_checked == 0
+            # A green gate that checked nothing is worse than a red one: a deleted or
+            # typo-named golden file must not be indistinguishable from a passing check.
+            @error "--check ran $n_no_golden case(s) but found no golden file for any of them — nothing was actually gated"
+            exit(1)
+        end
+        if opts.check && n_no_golden > 0
+            @warn "$n_no_golden of $(n_no_golden + n_checked) requested case(s) have no golden file and were not gated"
+        end
+        if n_untracked > 0
+            @info "$n_untracked tracked quantity/quantities have no golden value yet (not gating)"
+        end
         if opts.fail_on_change && n_changed > 0
             @error "$n_changed quantity/quantities changed (--fail-on-change)"
             exit(1)
diff --git a/regression-harness/src/extractor.jl b/regression-harness/src/extractor.jl
index 33a328d7c..2d1d632f2 100644
--- a/regression-harness/src/extractor.jl
+++ b/regression-harness/src/extractor.jl
@@ -112,6 +112,24 @@ function apply_extraction(spec::QuantitySpec, raw)::ExtractedQuantity
         json_str = JSON.json(pairs; allownan=true)
         return ExtractedQuantity(name, label, nothing, nothing, json_str, "json_array", threshold)
 
+    elseif startswith(spec.extract, "toml_key:")
+        # "toml_key:<dotted.path>": pin a deck's declared control value, read out of the
+        # `Input/gpec_toml_raw` rerun snapshot. Used to record which formalism a case asked
+        # for, because the Δ′ matrix is written to one canonical `SingularSurfaces/` path
+        # whichever integrator produced it, so the numbers alone no longer say. Reading it
+        # from `Input/` respects the schema rule that control values live nowhere else.
+        # A key the deck leaves unset pins as "<unset>", so a deck that starts or stops
+        # declaring one is itself a change rather than being silently absorbed.
+        keypath = spec.extract[(length("toml_key:")+1):end]
+        table = TOML.parse(raw isa AbstractString ? raw : String(raw))
+        node = table
+        for key in split(keypath, ".")
+            node = (node isa AbstractDict && haskey(node, key)) ? node[key] : nothing
+            node === nothing && break
+        end
+        token = node === nothing ? "<unset>" : string(node)
+        return ExtractedQuantity(name, label, nothing, nothing, token, "token", threshold)
+
     elseif spec.extract == "checksum"
         bytes = reinterpret(UInt8, vec(collect(raw)))
         hash = bytes2hex(sha256(bytes))
@@ -218,7 +236,7 @@ function compare_values(q1::NamedTuple, q2::NamedTuple)
         status = max_diff <= threshold ? "OK" : "CHANGED"
         return (max_diff, rel_diff, status)
 
-    elseif vtype == "checksum"
+    elseif vtype == "checksum" || vtype == "token"
         t1 = q1.value_text
         t2 = q2.value_text
         if t1 === nothing || t2 === nothing
diff --git a/regression-harness/src/golden.jl b/regression-harness/src/golden.jl
new file mode 100644
index 000000000..fb1dd1c30
--- /dev/null
+++ b/regression-harness/src/golden.jl
@@ -0,0 +1,528 @@
+"""
+Golden values: committed, reviewable reference numbers for a case.
+
+The `--refs` comparison answers "did this change?"; it cannot answer "is this right?", because
+both sides come from the same working copy of the code and nothing is recorded in git. Golden
+values close that gap: one TOML file per case, tracked alongside the source, holding the value
+of every quantity plus the tolerance it must be reproduced within and the evidence behind that
+tolerance.
+
+The tolerance policy is the load-bearing part. A tolerance is a claim about how much a quantity
+may legitimately move, so it is *derived* — from the residual drift at a convergence plateau and
+from the spread measured across platforms — never chosen to make a check pass. `update_golden`
+refuses to write a tolerance below the recorded platform spread, and a failing check is fixed by
+explaining the physics or fixing the regression, not by widening the bound.
+"""
+
+"""
+How strictly a quantity must reproduce, and why.
+
+  - `topological` — integer counts and mode numbers; exact equality
+  - `equilibrium_scalar` — spline/quadrature outputs with no adaptive branching; very tight
+  - `physics_converged` — the quantities golden values exist for (δW, Δ′, torque, growth
+    rates); tolerance comes from the measured convergence plateau and platform spread
+  - `diagnostic` — ODE step counts, runtimes; recorded and reported, never gating, because they
+    describe the numerics rather than the physics
+  - `unconverged` — measured and found to have no plateau; tracked differentially, never pinned,
+    so that a known-unconverged quantity is visibly excluded instead of quietly given a wide bound
+"""
+const TOLERANCE_CLASSES = ("topological", "equilibrium_scalar", "physics_converged", "diagnostic", "unconverged")
+
+"""Classes whose failure fails the run. `diagnostic` and `unconverged` are reported only."""
+const GATING_CLASSES = ("topological", "equilibrium_scalar", "physics_converged")
+
+"""
+Provisional per-class tolerances, used only when no measurement has been supplied.
+
+These are starting points, not results: a golden file written with them is stamped
+`tolerance_basis = "class-default (provisional)"` so that an unmeasured pin is never mistaken
+for a converged one.
+"""
+const CLASS_DEFAULT_TOLERANCE = Dict(
+    "topological" => (rtol=0.0, atol=0.0),
+    "equilibrium_scalar" => (rtol=1e-9, atol=0.0),
+    "physics_converged" => (rtol=1e-6, atol=0.0),
+    "diagnostic" => (rtol=Inf, atol=Inf),
+    "unconverged" => (rtol=Inf, atol=Inf)
+)
+
+"""
+One quantity's pinned value and the terms it is judged by.
+
+## Fields
+
+  - `name` — quantity key, matching the case's `[quantities.*]` table
+  - `value_type` — "real", "integer", "json_array" (mirrors `ExtractedQuantity`)
+  - `value_real` / `value_int` / `value_text` — the pinned value in its stored form
+  - `rtol` / `atol` — pass when `|x - gold| <= atol + rtol*|gold|`
+  - `class` — one of `TOLERANCE_CLASSES`
+  - `tolerance_basis` — how the tolerance was arrived at: "measured" or
+    "class-default (provisional)"
+  - `plateau_drift` — relative movement over the final refinement of the convergence scan; NaN
+    when not measured
+  - `platform_spread` — relative disagreement between reference platforms; NaN when not measured
+  - `converged_at` — the discretization the value was taken at, e.g.
+    "grid_type=ldp, mpsi=1024, mtheta=512"
+"""
+struct GoldenValue
+    name::String
+    value_type::String
+    value_real::Union{Float64,Nothing}
+    value_int::Union{Int,Nothing}
+    value_text::Union{String,Nothing}
+    rtol::Float64
+    atol::Float64
+    class::String
+    tolerance_basis::String
+    plateau_drift::Float64
+    platform_spread::Float64
+    converged_at::String
+end
+
+is_gating(g::GoldenValue) = g.class in GATING_CLASSES
+
+"""
+Provenance for a whole golden file: what produced these numbers and why they last changed.
+"""
+struct GoldenMeta
+    case::String
+    golden_version::Int
+    generated_at::String
+    commit::String
+    reason::String
+    julia_version::String
+    os_arch::String
+    manifest_sha::String
+    nthreads::Int
+    blas_threads::Int
+end
+
+golden_path(case_name::AbstractString) = joinpath(GOLDEN_DIR, "$(case_name).toml")
+
+has_golden(case_name::AbstractString) = isfile(golden_path(case_name))
+
+"""
+Read a golden file. Returns `(meta, Dict{name => GoldenValue})`, or `nothing` when absent.
+"""
+function load_golden(case_name::AbstractString)
+    path = golden_path(case_name)
+    isfile(path) || return nothing
+    data = TOML.parsefile(path)
+    m = get(data, "meta", Dict{String,Any}())
+    meta = GoldenMeta(
+        get(m, "case", String(case_name)),
+        Int(get(m, "golden_version", 0)),
+        get(m, "generated_at", ""),
+        get(m, "commit", ""),
+        get(m, "reason", ""),
+        get(m, "julia_version", ""),
+        get(m, "os_arch", ""),
+        get(m, "manifest_sha", ""),
+        Int(get(m, "nthreads", -1)),
+        Int(get(m, "blas_threads", -1))
+    )
+    values = Dict{String,GoldenValue}()
+    for (name, v) in get(data, "values", Dict{String,Any}())
+        class = get(v, "class", "physics_converged")
+        class in TOLERANCE_CLASSES || error("Golden file $path: quantity '$name' has unknown class '$class'")
+        # Enforced on every load, not only at write time: a hand edit or a merge taking the
+        # wrong side must not produce a gate quieter than the one save_golden refused to write.
+        if class in GATING_CLASSES
+            rt = Float64(get(v, "rtol", NaN))
+            (isfinite(rt) && rt >= 0) || error("Golden file $path: gating quantity '$name' has no finite rtol — malformed or hand-edited entry")
+            sp = Float64(get(v, "platform_spread", NaN))
+            isfinite(sp) && rt < sp && error("Golden file $path: quantity '$name' has rtol $rt below its recorded platform_spread $sp")
+        end
+        values[name] = GoldenValue(
+            name,
+            get(v, "value_type", "real"),
+            haskey(v, "value") && v["value"] isa Real ? Float64(v["value"]) : nothing,
+            haskey(v, "value_int") ? Int(v["value_int"]) : nothing,
+            haskey(v, "value_text") ? String(v["value_text"]) : nothing,
+            Float64(get(v, "rtol", NaN)),
+            Float64(get(v, "atol", 0.0)),
+            class,
+            get(v, "tolerance_basis", "unspecified"),
+            Float64(get(v, "plateau_drift", NaN)),
+            Float64(get(v, "platform_spread", NaN)),
+            get(v, "converged_at", "")
+        )
+    end
+    return (meta=meta, values=values)
+end
+
+"""
+Write a golden file.
+
+Refuses to emit a tolerance tighter than a recorded platform spread: a pin that no second
+platform can reproduce is a broken gate, and silently loosening it later is exactly the failure
+this whole mechanism exists to prevent. Raising the tolerance to cover a measured spread is a
+deliberate act the caller performs, not something this function does behind the caller's back.
+"""
+function save_golden(meta::GoldenMeta, values::Dict{String,GoldenValue})
+    for g in Base.values(values)
+        if g.value_real === nothing && g.value_int === nothing && g.value_text === nothing
+            error("Quantity '$(g.name)': refusing to write a golden entry with no value (the run " *
+                  "produced NaN or nothing) — a valueless gating entry fails forever and regenerating " *
+                  "reproduces it. Exclude the quantity or fix the extraction.")
+        end
+        if isfinite(g.platform_spread) && isfinite(g.rtol) && g.rtol < g.platform_spread
+            error("Quantity '$(g.name)': rtol $(g.rtol) is tighter than the measured platform " *
+                  "spread $(g.platform_spread). Widen it deliberately, with the measurement recorded, or " *
+                  "reduce the spread before pinning.")
+        end
+    end
+    mkpath(GOLDEN_DIR)
+    out = Dict{String,Any}(
+        "meta" => Dict{String,Any}(
+            "case" => meta.case,
+            "golden_version" => meta.golden_version,
+            "generated_at" => meta.generated_at,
+            "commit" => meta.commit,
+            "reason" => meta.reason,
+            "julia_version" => meta.julia_version,
+            "os_arch" => meta.os_arch,
+            "manifest_sha" => meta.manifest_sha,
+            "nthreads" => meta.nthreads,
+            "blas_threads" => meta.blas_threads
+        )
+    )
+    vals = Dict{String,Any}()
+    for (name, g) in values
+        entry = Dict{String,Any}(
+            "value_type" => g.value_type,
+            "class" => g.class,
+            "tolerance_basis" => g.tolerance_basis,
+            "rtol" => g.rtol,
+            "atol" => g.atol
+        )
+        g.value_real !== nothing && (entry["value"] = g.value_real)
+        g.value_int !== nothing && (entry["value_int"] = g.value_int)
+        g.value_text !== nothing && (entry["value_text"] = g.value_text)
+        isfinite(g.plateau_drift) && (entry["plateau_drift"] = g.plateau_drift)
+        isfinite(g.platform_spread) && (entry["platform_spread"] = g.platform_spread)
+        isempty(g.converged_at) || (entry["converged_at"] = g.converged_at)
+        vals[name] = entry
+    end
+    out["values"] = vals
+    open(golden_path(meta.case), "w") do io
+        TOML.print(io, out; sorted=true)
+    end
+    return golden_path(meta.case)
+end
+
+"""
+Classify a quantity from its case spec, so a new case gets sensible defaults without every
+tolerance having to be written by hand.
+
+Integer counts are topological; runtimes and step counts describe the numerics rather than the
+physics and so are diagnostic; everything else is assumed to be a physics quantity that must be
+converged, which is the conservative assumption — it gates.
+"""
+function infer_class(spec::QuantitySpec)::String
+    spec.type == "runtime" && return "diagnostic"
+    spec.name in ("nstep", "nstep_total") && return "diagnostic"
+    spec.type == "int_scalar" && return "topological"
+    # A declared control token is an exact-match gate, like a count: there is no tolerance
+    # that means anything between "riccati" and "galerkin".
+    startswith(spec.extract, "toml_key:") && return "topological"
+    name = spec.name
+    # sing_psi / sing_q are deliberately absent: singular-surface locations come from a root
+    # search, not pure spline/quadrature, so they take the measured physics_converged path.
+    equilibrium_names = ("q0", "q95", "betat", "betan", "betap1", "betap2", "betap3", "betaj",
+        "li1", "li2", "li3", "volume", "crnt", "bt0", "bwall", "aratio", "kappa")
+    name in equilibrium_names && return "equilibrium_scalar"
+    return "physics_converged"
+end
+
+"""
+Compare one extracted quantity against its golden value.
+
+Returns `(passed, deviation, detail)` where `deviation` is the relative deviation (absolute when
+the golden value is zero) and `detail` is a human-readable reason on failure. Array quantities
+pass only when every element is within tolerance; the reported deviation is the worst element.
+"""
+function compare_to_golden(q::NamedTuple, g::GoldenValue)
+    # Non-finite tolerances mean "recorded, never judged" (diagnostic/unconverged); without this
+    # guard, gold == 0 turns atol + rtol*|gold| into Inf + NaN and the comparison is false.
+    within = (x, gold) -> !isfinite(g.atol) || !isfinite(g.rtol) ||
+        abs(x - gold) <= g.atol + g.rtol * abs(gold)
+
+    if q.value_type != g.value_type
+        return (false, NaN, "type changed: golden $(g.value_type), got $(q.value_type)")
+    end
+
+    if g.value_type == "integer"
+        (q.value_int === nothing || g.value_int === nothing) && return (false, NaN, "missing integer value")
+        d = Float64(abs(q.value_int - g.value_int))
+        passed = g.class == "topological" ? q.value_int == g.value_int : within(Float64(q.value_int), Float64(g.value_int))
+        return (passed, g.value_int == 0 ? d : d / abs(g.value_int), passed ? "" : "$(g.value_int) → $(q.value_int)")
+
+    elseif g.value_type == "real"
+        (q.value_real === nothing || g.value_real === nothing) && return (false, NaN, "missing value")
+        gold = g.value_real
+        got = q.value_real
+        d = abs(got - gold)
+        rel = gold == 0.0 ? d : d / abs(gold)
+        return (within(got, gold), rel, within(got, gold) ? "" : @sprintf("%.9g → %.9g", gold, got))
+
+    elseif g.value_type == "json_array"
+        (q.value_text === nothing || g.value_text === nothing) && return (false, NaN, "missing array")
+        got = JSON.parse(q.value_text; allownan=true)
+        gold = JSON.parse(g.value_text; allownan=true)
+        length(got) == length(gold) && return _compare_arrays(got, gold, g, within)
+        return (false, NaN, "length $(length(gold)) → $(length(got))")
+
+    elseif g.value_type == "token"
+        # A token names a discrete choice (which integrator produced the Δ′), so there is no
+        # tolerance to apply: it matches or the run is answering a different question than the
+        # gold does. Deviation is reported as 1.0 on mismatch rather than NaN so it sorts as a
+        # real failure in reports.
+        (q.value_text === nothing || g.value_text === nothing) && return (false, NaN, "missing token")
+        matched = q.value_text == g.value_text
+        return (matched, matched ? 0.0 : 1.0, matched ? "" : "$(g.value_text) → $(q.value_text)")
+    end
+
+    return (false, NaN, "unsupported value type $(g.value_type)")
+end
+
+"""Worst-element comparison for array quantities, shared by the real and complex encodings."""
+function _compare_arrays(got, gold, g::GoldenValue, within)
+    worst_rel = 0.0
+    worst_idx = 0
+    all_ok = true
+    for i in eachindex(gold)
+        a = _json_element_abs(gold[i])
+        d = _json_element_diff(gold[i], got[i])
+        rel = a == 0.0 ? d : d / a
+        ok = d <= g.atol + g.rtol * a
+        ok || (all_ok = false)
+        if rel > worst_rel
+            worst_rel = rel
+            worst_idx = i
+        end
+    end
+    detail = all_ok ? "" : "worst element $(worst_idx): rel $(@sprintf("%.3e", worst_rel))"
+    return (all_ok, worst_rel, detail)
+end
+
+"""
+Compare a case's fresh run against its golden file and print the report.
+
+Returns `(n_pass, n_fail, n_untracked, n_informational)`. Only gating classes can fail;
+`diagnostic` and `unconverged` quantities are shown with their deviation and marked as
+informational, so a reader sees them move without the run failing on them.
+"""
+function report_golden_check(db::SQLite.DB, case_spec::CaseSpec, commit_hash::String)
+    golden = load_golden(case_spec.name)
+    println()
+    println("Golden Check: $(case_spec.name)")
+    if golden === nothing
+        println("  No golden file at $(golden_path(case_spec.name)) — nothing to check against.")
+        println("  Generate one with:  regress --update-golden --cases $(case_spec.name) --reason \"...\"")
+        return (n_pass=0, n_fail=0, n_untracked=0, n_informational=0)
+    end
+
+    quantities = get_quantities(db, commit_hash, case_spec.name)
+    info = get_run_info(db, commit_hash, case_spec.name)
+    if info !== nothing && !info.success
+        println("  RUN FAILED — nothing to compare (this is a crash, not a tolerance failure):")
+        println("  $(_short_err(info.error_msg))")
+        return (n_pass=0, n_fail=1, n_untracked=0, n_informational=0)
+    end
+
+    rows = Vector{Vector{String}}()
+    n_pass = n_fail = n_untracked = n_informational = 0
+
+    for spec in case_spec.quantities
+        g = get(golden.values, spec.name, nothing)
+        if g === nothing
+            # Runtime and checksums are structurally un-goldenable (no tolerance semantics), so
+            # they must not inflate the untracked count that flags genuinely unpinned physics.
+            (spec.type == "runtime" || spec.extract == "checksum") && continue
+            n_untracked += 1
+            continue
+        end
+        q_raw = get(quantities, spec.name, nothing)
+        if q_raw === nothing
+            push!(rows, [spec.label, g.class, "MISSING", "—", "FAIL"])
+            n_fail += 1
+            continue
+        end
+        # SQLite NULLs surface as , which the === nothing guards in compare_to_golden
+        # never match; normalize here as the update path already does.
+        q = (label=q_raw.label, value_real=_column(q_raw.value_real, nothing),
+             value_int=_column(q_raw.value_int, nothing), value_text=_column(q_raw.value_text, nothing),
+             value_type=q_raw.value_type, noise_threshold=q_raw.noise_threshold)
+        passed, deviation, detail = compare_to_golden(q, g)
+        gating = is_gating(g)
+        status = if !gating
+            "info"
+        elseif passed
+            "ok"
+        else
+            "** FAIL **"
+        end
+        if !gating
+            n_informational += 1
+        elseif passed
+            n_pass += 1
+        else
+            n_fail += 1
+        end
+        dev = isnan(deviation) ? "—" : @sprintf("%.2e", deviation)
+        tol = isfinite(g.rtol) ? @sprintf("%.1e", g.rtol) : "—"
+        push!(rows, [spec.label, g.class, dev, tol, status * (isempty(detail) ? "" : "  $detail")])
+    end
+
+    # Golden entries with no matching case quantity fail regardless of class: a renamed or
+    # removed quantity would otherwise silently delete its gate (the new name enters as merely
+    # "untracked"), and a check that quietly stops checking something is worse than one that
+    # fails. The fix is to regenerate the goldens for the reshaped case, with a reason.
+    spec_names = Set(spec.name for spec in case_spec.quantities)
+    for name in sort(collect(keys(golden.values)))
+        name in spec_names && continue
+        push!(rows, [name, golden.values[name].class, "ORPHANED", "—",
+            "** FAIL **  no matching quantity in the case — renamed or removed? Regenerate goldens."])
+        n_fail += 1
+    end
+
+    header = ["Quantity", "Class", "Deviation", "rtol", "Status"]
+    widths = [length(h) for h in header]
+    for row in rows, i in eachindex(row)
+        widths[i] = max(widths[i], length(row[i]))
+    end
+    total_w = sum(widths) + 2 * (length(header) - 1)
+
+    println("="^total_w)
+    println("Golden v$(golden.meta.golden_version), generated $(golden.meta.generated_at) @ $(golden.meta.commit)")
+    isempty(golden.meta.reason) || println("Reason: $(golden.meta.reason)")
+    println("Golden env: julia $(golden.meta.julia_version), $(golden.meta.os_arch), $(golden.meta.nthreads) thread/$(golden.meta.blas_threads) BLAS")
+    if info !== nothing && !isempty(info.fingerprint.julia_version)
+        println("This run:   $(describe_env(info.fingerprint))")
+        info.fingerprint.os_arch == golden.meta.os_arch ||
+            println("  NOTE: different platform from the one the goldens were measured on; deviations at or below the recorded platform_spread are expected.")
+    end
+    println("-"^total_w)
+    _print_row(header, widths)
+    println("-"^total_w)
+    for row in rows
+        _print_row(row, widths)
+    end
+    println("="^total_w)
+    parts = String[]
+    n_pass > 0 && push!(parts, "$n_pass pass")
+    n_fail > 0 && push!(parts, "$n_fail FAIL")
+    n_informational > 0 && push!(parts, "$n_informational informational")
+    n_untracked > 0 && push!(parts, "$n_untracked untracked")
+    println("Summary: ", join(parts, ", "))
+    println()
+    return (n_pass=n_pass, n_fail=n_fail, n_untracked=n_untracked, n_informational=n_informational)
+end
+
+"""
+Regenerate a case's golden file from a completed run, reporting what moved.
+
+Prints an old→new delta for every quantity whose value changed, so the diff a reviewer sees in
+git is accompanied by the size of each move.
+"""
+function update_golden_from_run(db::SQLite.DB, case_spec::CaseSpec, commit_hash::String,
+                                reason::String, repo_root::String)
+    info = get_run_info(db, commit_hash, case_spec.name)
+    (info === nothing || !info.success) && error("Cannot update goldens for '$(case_spec.name)': the run did not succeed")
+
+    # SQLite NULLs come back as `missing`, which the Union{...,Nothing} fields reject; `_column`
+    # normalizes both absent forms to nothing.
+    extracted = ExtractedQuantity[]
+    for (name, q) in get_quantities(db, commit_hash, case_spec.name)
+        push!(extracted, ExtractedQuantity(name, String(_column(q.label, name)),
+            _column(q.value_real, nothing), _column(q.value_int, nothing),
+            _column(q.value_text, nothing), q.value_type, q.noise_threshold))
+    end
+
+    previous = load_golden(case_spec.name)
+    existing = previous === nothing ? nothing : previous.values
+    values = build_golden_values(extracted, case_spec.quantities, existing)
+
+    fp = info.fingerprint
+    meta = GoldenMeta(case_spec.name,
+        previous === nothing ? 1 : previous.meta.golden_version + 1,
+        Dates.format(Dates.now(), "yyyy-mm-dd"),
+        strip(read(`git -C $repo_root rev-parse --short HEAD`, String)),
+        reason, fp.julia_version, fp.os_arch, fp.manifest_sha, fp.nthreads, fp.blas_threads)
+    if commit_hash == LOCAL_REF && !isempty(strip(read(`git -C $repo_root status --porcelain`, String)))
+        # The numbers came from uncommitted source; a clean HEAD checkout will not reproduce
+        # them, and the commit field is the one a reviewer uses to reproduce a disputed number.
+        meta = GoldenMeta(meta.case, meta.golden_version, meta.generated_at,
+            meta.commit * "-dirty", meta.reason, meta.julia_version, meta.os_arch,
+            meta.manifest_sha, meta.nthreads, meta.blas_threads)
+        @warn "Working tree is dirty: golden provenance recorded as $(meta.commit). Commit first if these values are meant to be reproducible."
+    end
+
+    println()
+    println("Golden update: $(case_spec.name)  (v$(previous === nothing ? 0 : previous.meta.golden_version) → v$(meta.golden_version))")
+    if existing !== nothing
+        for name in sort(collect(keys(existing)))
+            if !haskey(values, name)
+                println(@sprintf("  %-34s REMOVED — extraction returned missing (renamed h5 path?) or the quantity left the case. This deletes its gate; confirm it is intentional.", name))
+            end
+        end
+        for (name, g) in sort(collect(values); by=first)
+            old = get(existing, name, nothing)
+            old === nothing && (println(@sprintf("  %-34s NEW", name)); continue)
+            if g.value_real !== nothing && old.value_real !== nothing && g.value_real != old.value_real
+                rel = old.value_real == 0 ? Inf : abs(g.value_real - old.value_real) / abs(old.value_real)
+                println(@sprintf("  %-34s %.9g → %.9g   (rel %.2e)", name, old.value_real, g.value_real, rel))
+            elseif g.value_text !== nothing && old.value_text !== nothing && g.value_text != old.value_text
+                println(@sprintf("  %-34s array changed", name))
+            end
+        end
+    end
+    path = save_golden(meta, values)
+    println("Wrote $path  ($(length(values)) quantities)")
+    provisional = count(g -> startswith(g.tolerance_basis, "class-default"), Base.values(values))
+    provisional > 0 && println("  $provisional quantity/quantities still carry provisional class-default tolerances — " *
+                               "replace them with measured plateau drift and platform spread before relying on this as a gate.")
+    println()
+    return path
+end
+
+"""
+Build golden entries from a set of freshly extracted quantities.
+
+`existing` carries forward the tolerance, class, and evidence already recorded for a quantity, so
+regenerating values after a physics change does not silently reset hard-won measurements to
+provisional defaults.
+"""
+function build_golden_values(extracted::Vector{ExtractedQuantity}, specs::Vector{QuantitySpec},
+                             existing::Union{Dict{String,GoldenValue},Nothing})
+    spec_by_name = Dict(s.name => s for s in specs)
+    values = Dict{String,GoldenValue}()
+    for eq in extracted
+        eq.value_type == "missing" && continue
+        spec = get(spec_by_name, eq.name, nothing)
+        spec === nothing && continue
+        prior = existing === nothing ? nothing : get(existing, eq.name, nothing)
+        if prior !== nothing && prior.value_type != eq.value_type
+            # A type change invalidates the class and every measurement made under the old type;
+            # carrying a topological rtol=0 onto a float (or a float rtol onto a count) mis-gates.
+            @warn "Golden '$(eq.name)': value_type changed $(prior.value_type) → $(eq.value_type); resetting class and tolerances to provisional"
+            prior = nothing
+        end
+        class = prior === nothing ? infer_class(spec) : prior.class
+        # Checksums have no notion of "close", so they cannot carry a tolerance; they stay a
+        # same-machine differential tool rather than a golden gate.
+        eq.value_type == "checksum" && continue
+        if prior === nothing
+            defaults = CLASS_DEFAULT_TOLERANCE[class]
+            rtol, atol, basis = defaults.rtol, defaults.atol, "class-default (provisional)"
+            drift, spread, at = NaN, NaN, ""
+        else
+            rtol, atol, basis = prior.rtol, prior.atol, prior.tolerance_basis
+            drift, spread, at = prior.plateau_drift, prior.platform_spread, prior.converged_at
+        end
+        values[eq.name] = GoldenValue(eq.name, eq.value_type, eq.value_real, eq.value_int,
+            eq.value_text, rtol, atol, class, basis, drift, spread, at)
+    end
+    return values
+end
diff --git a/regression-harness/src/runner.jl b/regression-harness/src/runner.jl
index ed7e5f4d5..af8294ad1 100644
--- a/regression-harness/src/runner.jl
+++ b/regression-harness/src/runner.jl
@@ -121,6 +121,80 @@ end
 %RUNINFO%
 """
 
+
+# Fixed-Q probe of the SLAYER inner-layer dispersion Δ(Q) on the DIII-D-like 2/1 surface.
+# The adaptive AMR scan samples in gpec.h5 cannot be pinned (sample locations move under any
+# refinement change), so this evaluates Δ(Q) on a fixed 4×4 grid over Re(Q), Im(Q) ∈ [-10, 10]
+# instead — a refinement-stable pin of the dispersion curve itself. The layer parameters are the
+# DIII-D-like SLAYER deck's own 2/1 surface values (Tearing/PerSurface, develop @ 1f193a6c),
+# quoted verbatim so the case is self-contained and probes the SOLVER alone: the parameter
+# chain that produces these numbers is pinned separately by the diiid_slayer_n1 case.
+const COMPUTED_SLAYER_DELTA_PROBE_SCRIPT_TEMPLATE = """
+using Pkg
+%INSTANTIATE%
+using GeneralizedPerturbedEquilibrium
+using GeneralizedPerturbedEquilibrium.InnerLayer
+using HDF5
+p = SLAYERParameters(;
+    tau=1.1975430647804235, lu=6.086905739791344e6, c_beta=0.22094021004591707,
+    D_norm=4.191469284125091, P_perp=50.2972642329308, P_tor=34.74437891503841,
+    Q_e=1.0915286815773122, Q_i=-1.6558720999124832, iota_e=0.39729503206497013,
+    tauk=0.00010358784131467763, tau_r=3.453343933553279, delta_n=504.745127277822,
+    rs=0.3617373814196757, R0=1.7433359412007365, bt=1.0, sval_r=1.260093929519795,
+    eta=4.761642777337999e-8, d_beta=0.011249087118484661)
+axis = range(-10.0, 10.0; length=4)
+Q = ComplexF64[re + im_ * 1im for im_ in axis for re in axis]
+t_start = time()
+Δ = ComplexF64[solve_inner(SLAYERModel(), p, q).tearing for q in Q]
+elapsed = time() - t_start
+h5open(ARGS[1], "w") do fid
+    fid["slayer_probe/Q_re"]     = real.(Q)
+    fid["slayer_probe/Q_im"]     = imag.(Q)
+    fid["slayer_probe/Delta_re"] = real.(Δ)
+    fid["slayer_probe/Delta_im"] = imag.(Δ)
+end
+%RUNINFO%
+"""
+# External-reference validation: GPEC's del_s Riccati solver against Fitzpatrick, "Tearing Mode
+# Dynamics in Tokamak Plasmas" (IOP 2023), figures 6.2 and 6.3. Unlike every other case in this
+# suite -- which pins GPEC against its own previous output -- this one pins it against published
+# values, so a failure means "we no longer reproduce the textbook", not "we changed".
+#
+# riccati_del_s integrates the book's layer equations directly: E and F are (6.29)/(6.30), the
+# Riccati form is (6.38), and the returned dels_db is the plotted delta_s/d_beta. Prescribing the
+# normalized parameters (D_norm = 1, P_perp = P_tor = Phat, Q_e = Qhat/(1+1/tau)) makes the solver's
+# internal Q_hat equal the book's Qhat_*, so the grid below is exactly the figures' axes.
+#
+# tau is NOT stated in the figure captions; tau = 1 (table 5.1, low-field reactor) is pinned here as
+# an explicit assumption of the case, not inherited from a default that could drift.
+#
+# Phat = 0 is a singular edge of the model, not a hard case: alpha = sqrt(Phat_perp/(1+1/tau))
+# vanishes so the large-q boundary form degenerates, and F -> -i*Qhat at the origin. The grid starts
+# just inside so the pins record results rather than a modelling boundary.
+const COMPUTED_SLAYER_DELS_FITZPATRICK_SCRIPT_TEMPLATE = """
+using Pkg
+%INSTANTIATE%
+using GeneralizedPerturbedEquilibrium
+using GeneralizedPerturbedEquilibrium.InnerLayer
+using HDF5
+const TAU = 1.0
+probe(Q, P) = SLAYERParameters(; tau=TAU, lu=1.0, c_beta=0.0, D_norm=1.0,
+    P_perp=P, P_tor=P, Q_e=Q / (1 + 1/TAU), Q_i=0.0, iota_e=0.0, tauk=1.0,
+    tau_r=1.0, delta_n=1.0, rs=1.0, R0=1.0, bt=1.0, sval_r=1.0, eta=1.0, d_beta=1.0)
+axis = [0.02, 0.5, 1.0, 2.0, 4.0]
+QP = [(q, p) for p in axis for q in axis]
+t_start = time()
+dels = ComplexF64[riccati_del_s(probe(q, p)) for (q, p) in QP]
+elapsed = time() - t_start
+h5open(ARGS[1], "w") do fid
+    fid["fitzpatrick/Q_hat"]  = Float64[q for (q, _) in QP]
+    fid["fitzpatrick/P_hat"]  = Float64[p for (_, p) in QP]
+    fid["fitzpatrick/dels_db_re"] = real.(dels)
+    fid["fitzpatrick/dels_db_im"] = imag.(dels)
+    fid["fitzpatrick/tau"] = TAU
+end
+%RUNINFO%
+"""
 # Self-contained separatrix-finder regression (PR #296). Loads a fixed-boundary EFIT whose
 # computational box hugs the LCFS (eps=0.05 TokaMaker aspect-scan g-file): outside the prescribed
 # LCFS the coil-vacuum flux turns back above the boundary value before the grid edge, so the old
@@ -253,8 +327,12 @@ function _computed_script_template(case_spec::CaseSpec)
         return COMPUTED_GGJ_SCRIPT_TEMPLATE
     elseif case_spec.name == "ggj_ray_q500i"
         return COMPUTED_GGJ_RAY_SCRIPT_TEMPLATE
+    elseif case_spec.name == "slayer_delta_probe"
+        return COMPUTED_SLAYER_DELTA_PROBE_SCRIPT_TEMPLATE
     elseif case_spec.name == "efit_fixedbdy_separatrix"
         return COMPUTED_SEPARATRIX_SCRIPT_TEMPLATE
+    elseif case_spec.name == "slayer_dels_fitzpatrick"
+        return COMPUTED_SLAYER_DELS_FITZPATRICK_SCRIPT_TEMPLATE
     end
     error("No computed-script template registered for case '$(case_spec.name)'")
 end
diff --git a/regression-harness/src/types.jl b/regression-harness/src/types.jl
index 05c7159ce..4c7393b06 100644
--- a/regression-harness/src/types.jl
+++ b/regression-harness/src/types.jl
@@ -8,8 +8,9 @@ Specification for a single quantity to extract from gpec.h5.
 struct QuantitySpec
     name::String
     h5path::String          # HDF5 dataset path (e.g. "ForceFreeStates/FreeBoundaryStability/eigenmode_energies"), empty for runtime
-    type::String            # "complex_vector", "real_vector", "real_scalar", "int_scalar", "real_matrix", "runtime"
-    extract::String         # "value", "real_first", "imag_first", "abs_first", "norm", "all_real", "all_complex", "checksum"
+    type::String            # "complex_vector", "real_vector", "real_scalar", "int_scalar", "real_matrix", "token", "runtime"
+    extract::String         # "value", "real_first", "imag_first", "abs_first", "norm", "all_real", "all_complex",
+    # "diagonal_complex", "first_<N>", "first_<N>_complex", "checksum", "toml_key:<dotted.path>"
     label::String           # Human-readable label for reports
     noise_threshold::Float64 # Absolute changes below this are noise
     order::Int              # Display order in reports (lower = earlier)
@@ -53,6 +54,11 @@ Parsed CLI options.
 worktree (pinning is on by default, so that two refs differ only by source code).
 `allow_env_mismatch` lets a cached result from a different environment be reused instead of
 re-run. `fail_on_change` turns any changed quantity into a non-zero exit status, for CI use.
+
+`check` runs the working tree and compares it against the committed golden values, which is the
+mode CI gates on. `update_golden` regenerates those values from a fresh run; `reason` records
+why, and is mandatory because a golden change is a claim about physics that a reviewer has to be
+able to evaluate.
 """
 struct CLIOptions
     cases::Vector{String}
@@ -68,5 +74,8 @@ struct CLIOptions
     no_pin_manifest::Bool
     allow_env_mismatch::Bool
     fail_on_change::Bool
+    check::Bool
+    update_golden::Bool
+    reason::Union{String,Nothing}
     help::Bool
 end
diff --git a/src/InnerLayer/SLAYER/LayerInputs.jl b/src/InnerLayer/SLAYER/LayerInputs.jl
index 96177903d..766a689ef 100644
--- a/src/InnerLayer/SLAYER/LayerInputs.jl
+++ b/src/InnerLayer/SLAYER/LayerInputs.jl
@@ -90,8 +90,10 @@ profiles, without an intermediate file round-trip.
   - `zeff`      -- effective charge (default `1.0`).
   - `chi_perp`  -- perpendicular heat diffusivity [m²/s]. Scalar or a
     callable of `psi` (default `1.0`).
-  - `chi_tor`   -- toroidal heat diffusivity [m²/s]. Scalar or a callable
-    of `psi` (default `1.0`).
+  - `chi_tor`   -- Fitzpatrick's χ_φ: anomalous perpendicular ion MOMENTUM
+    diffusivity `[m²/s]` (a viscosity entering the vorticity and parallel-flow
+    equations — not a heat diffusivity). Scalar or a callable of `psi`
+    (default `1.0`).
   - `dr_val`    -- resistive interchange index `D_R = E + F + H²`
     (Glasser-Greene-Johnson 1975) feeding the critical-Δ formulas
     (`:lar`, `:rfitzp`, `:toroidal`). When `nothing` (default), Julia
diff --git a/src/InnerLayer/SLAYER/LayerParameters.jl b/src/InnerLayer/SLAYER/LayerParameters.jl
index 5dbf9163f..f8a36f87b 100644
--- a/src/InnerLayer/SLAYER/LayerParameters.jl
+++ b/src/InnerLayer/SLAYER/LayerParameters.jl
@@ -29,7 +29,7 @@ de-normalization. The parametrization uses `P_perp`, `P_tor`, and
 | `c_beta`   | Compressibility √(β_local / (1 + β_local))                        |
 | `D_norm`   | (d_β/r_s) · S^(1/3) · √(τ/(1+τ))  (Fitzpatrick normalized scale)  |
 | `P_perp`   | Perpendicular Prandtl number τ_R / τ_⊥                            |
-| `P_tor`    | Toroidal-direction Prandtl number τ_R / τ_‖tor                    |
+| `P_tor`    | Momentum (magnetic) Prandtl number τ_R / τ_φ, τ_φ = r_s²/χ_φ      |
 | `Q_e`      | Normalized electron diamagnetic: −tauk · ω_*e                     |
 | `Q_i`      | Normalized ion diamagnetic:      −tauk · ω_*i                     |
 | `iota_e`   | Q_e / (Q_e − Q_i)                                                 |
@@ -197,7 +197,10 @@ parametrization (P_perp/P_tor/D_norm; the older magnetic/electron Prandtl
   - `R0`      -- major radius [m]
   - `mu_i`    -- ion mass in proton-mass units (e.g. 2.0 for D)
   - `zeff`    -- effective charge
-  - `chi_perp`, `chi_tor` -- perpendicular / toroidal heat diffusivity [m²/s]
+  - `chi_perp` -- perpendicular energy diffusivity χ⊥ [m²/s]
+  - `chi_tor`  -- Fitzpatrick's χ_φ, the anomalous perpendicular ion momentum
+    diffusivity `[m²/s]` (a viscosity, not heat; TJ Layer.tex eqs. — enters the
+    vorticity and parallel-flow equations, τ_φ = r_s²/χ_φ)
   - `m`, `n`  -- poloidal / toroidal mode numbers at the surface
   - `dr_val`, `dgeo_val` -- inputs for the critical-Δ formula
   - `dc_type` -- one of `:none`, `:lar`, `:rfitzp`, `:toroidal`
diff --git a/src/Tearing/Runner/Control.jl b/src/Tearing/Runner/Control.jl
index 1a12ff7e6..f41777468 100644
--- a/src/Tearing/Runner/Control.jl
+++ b/src/Tearing/Runner/Control.jl
@@ -34,10 +34,13 @@ constructor.
     callable of `psi` overrides it
   - `mu_i`     -- ion mass in proton-mass units (default 2.0 for D)
   - `zeff`     -- effective charge
-  - `chi_perp`, `chi_tor` -- fallback perpendicular / toroidal heat
-    diffusivity [m²/s], used only when the kinetic file carries no usable
-    `chi_e`/`chi_phi` profile (dataset absent or all-zero); otherwise the
-    file's χ⊥(ψ)/χ_φ(ψ) take precedence
+  - `chi_perp`, `chi_tor` -- fallback scalars [m²/s], used only when the
+    kinetic file carries no usable `chi_e`/`chi_phi` profile (dataset absent
+    or all-zero); otherwise the file's χ⊥(ψ)/χ_φ(ψ) take precedence.
+    `chi_perp` is the perpendicular ENERGY diffusivity χ⊥; `chi_tor` is
+    Fitzpatrick's χ_φ, the anomalous perpendicular ion MOMENTUM diffusivity
+    (a viscosity: it enters the vorticity and parallel-flow equations, TJ
+    Layer.tex — not a heat diffusivity)
   - `dr_val`, `dgeo_val`  -- critical-Δ formula inputs. `nothing` (default)
     auto-derives them from the equilibrium: `dr_val` from the resistive
     interchange index `D_R = E + F + H²` at each surface, `dgeo_val` from the