diff --git a/.github/workflows/test.yaml b/.github/workflows/test.yaml
index a5e3561d5..2afb787a7 100644
--- a/.github/workflows/test.yaml
+++ b/.github/workflows/test.yaml
@@ -86,7 +86,9 @@ jobs:
           fi
 
   test:
-    name: runtests ${{ matrix.version }} - ${{ matrix.os }}
+    # Branch protection lists the single-threaded job names as required checks, so those must
+    # stay byte-identical; the multi-threaded leg gets a distinct name and is purely additive.
+    name: runtests ${{ matrix.version }} - ${{ matrix.os }}${{ matrix.threads != '1' && format(' ({0} threads)', matrix.threads) || '' }}
     needs: changes
     if: needs.changes.outputs.julia == 'true'
     runs-on: ${{ matrix.os }}
@@ -99,8 +101,17 @@ jobs:
           - '1.x'  # latest (currently 1.12)
         os:
           - ubuntu-latest
+        threads:
+          - '1'
+        include:
+          # Threads.@threads executes serially at one thread, so a single-threaded matrix never
+          # exercises threaded scheduling; this leg covers real concurrent execution.
+          - version: '1.11'
+            os: ubuntu-latest
+            threads: '4'
 
     env:
+      JULIA_NUM_THREADS: ${{ matrix.threads }}
       DEPOT_PATHS: |
         ~/.julia/artifacts
         ~/.julia/packages
diff --git a/CLAUDE.md b/CLAUDE.md
index 9417068b5..99886fe00 100644
--- a/CLAUDE.md
+++ b/CLAUDE.md
@@ -22,16 +22,22 @@ GPEC (Generalized Perturbed Equilibrium Code, Julia implementation) is a compreh
 # Run all tests
 julia --project=. -e 'using Pkg; Pkg.activate("."); Pkg.instantiate(); include("test/runtests.jl")'
 
-# Run specific test file
-julia --project=. test/runtests.jl test/runtests_solovev.jl
-
-# Available test files:
-
-# - test/runtests_vacuum_julia.jl   # Julia vacuum module
-# - test/runtests_solovev.jl        # Analytical equilibrium
-# - test/runtests_ode.jl            # ODE integration
-# - test/runtests_sing.jl           # Singular surface handling
-# - test/runtests_fullruns.jl       # End-to-end tests
+# Run specific test file — the argument is included relative to test/, so pass the bare
+# filename, not a path prefixed with test/
+julia --project=. test/runtests.jl runtests_sing.jl
+
+# Run the suite multi-threaded (the parallel FM/BVP paths reduce to their serial form at one
+# thread, so a single-threaded run never exercises threaded execution)
+julia -t 4 --project=. test/runtests.jl
+
+# A few of the available test files (see test/runtests.jl for the full list):
+
+# - test/runtests_vacuum.jl             # Vacuum module
+# - test/runtests_equil.jl              # Equilibrium reconstruction
+# - test/runtests_sing.jl               # Singular surface handling
+# - test/runtests_parallel_integration.jl  # Parallel FM integration and BVP Delta'
+# - test/runtests_decomposition_invariance.jl  # Riccati Δ' chunk-decomposition invariance
+# - test/runtests_fullruns.jl           # End-to-end tests
 ```
 
 ### Building Documentation
diff --git a/test/runtests.jl b/test/runtests.jl
index 48a53e2d6..9a0eeafa2 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -32,6 +32,7 @@ else
     include("./runtests_parallel_integration.jl")
     include("./runtests_result_struct.jl")
     include("./runtests_solve_api.jl")
+    include("./runtests_decomposition_invariance.jl")
     include("./runtests_sing.jl")
     include("./runtests_innerlayer.jl")
     include("./runtests_tj_analytic.jl")
diff --git a/test/runtests_decomposition_invariance.jl b/test/runtests_decomposition_invariance.jl
new file mode 100644
index 000000000..097af281d
--- /dev/null
+++ b/test/runtests_decomposition_invariance.jl
@@ -0,0 +1,77 @@
+using Test
+using TOML
+
+# Decomposition invariance of the Riccati Δ′ path.
+#
+# The chunked propagator driver reassociates the fundamental-matrix products whenever the chunk
+# decomposition changes, and floating-point products do not reassociate exactly in general — so
+# Δ′ reproducibility is a real end-to-end claim, beyond the structural chunk-count invariance
+# pinned by the unit tests in runtests_parallel_integration.jl. This file asserts it: the Δ′
+# matrix is bit-identical when the same integration is cut into a genuinely different set of
+# chunks. Driven through the public solve API so it exercises the production path itself.
+
+const GP_TI = GeneralizedPerturbedEquilibrium
+
+"""
+Solve the DIII-D-like deck with the Riccati integrator at a given chunk count (`nchunks = 0` is
+the msing-derived auto target) and return the published Δ′ matrix alongside the leading energy
+eigenvalue, which is used only as a witness that the two decompositions really differed.
+"""
+function _solve_at_nchunks(dir::String, nchunks::Int)
+    inputs = TOML.parsefile(joinpath(dir, "gpec.toml"))
+    ffs_in = inputs["ForceFreeStates"]
+    eq_config = GP_TI.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], dir)
+    equil = GP_TI.Equilibrium.setup_equilibrium(eq_config, nothing)
+    if GP_TI.Equilibrium.wants_two_pass(eq_config)
+        mand = GP_TI.ForceFreeStates.rational_psi_nodes(equil; nlow=ffs_in["nn_low"], nhigh=ffs_in["nn_high"])
+        psi_nodes = GP_TI.Equilibrium.refined_psi_grid(equil; tau=eq_config.psi_accuracy, mandatory=mand)
+        rerun_input = GP_TI.Equilibrium.build_direct_from_ingest(eq_config, equil.ingest)
+        equil = GP_TI.Equilibrium.setup_equilibrium(eq_config, rerun_input; override_psi_nodes=psi_nodes)
+    end
+
+    # Every ForceFreeStatesControl key from the deck except the ones the problem or the
+    # integrator owns: nn_low/nn_high come from `nn`, nchunks and the formalism from the alg.
+    ctrl_kwargs = Dict(Symbol(k) => v for (k, v) in ffs_in
+                       if !(k in ("nn_low", "nn_high", "nchunks", "integrator")))
+    ctrl_kwargs[:verbose] = false
+    ctrl_kwargs[:write_outputs_to_HDF5] = false
+
+    wall = GP_TI.Vacuum.WallShapeSettings(; (Symbol(k) => v for (k, v) in inputs["Wall"])...)
+    prob = GP_TI.EulerLagrangeProblem(equil; nn=ffs_in["nn_low"], wall=wall, dir_path=dir, ctrl_kwargs...)
+    return GP_TI.solve(prob, GP_TI.ForceFreeStates.Riccati(; nchunks=nchunks))
+end
+
+@testset "Decomposition invariance of the Riccati Δ′ path" begin
+    # The deck whose Δ′ the regression harness pins, and where the BVP Δ′ is well-conditioned
+    # (Solovev sits near marginal stability and its BVP Δ′ is pathological there).
+    dir = joinpath(@__DIR__, "..", "examples", "DIIID-like_ideal_example")
+    auto = _solve_at_nchunks(dir, 0)
+    @test auto.delta_prime !== nothing
+    msing = size(auto.delta_prime.matrix, 1)
+
+    # The nchunks=0 target, mirroring balance_integration_chunks' internal formula (as
+    # runtests_parallel_integration.jl does). Invariance is asserted ABOVE this floor only:
+    # fewer chunks than the msing-derived minimum is not a different decomposition but a
+    # structurally deficient one (the floor gives the rational-surface crossings room).
+    auto_target = max(2 * msing + 3, 8 * (msing + 1) + msing)
+    finer = _solve_at_nchunks(dir, auto_target + 11)
+    @test finer.delta_prime !== nothing
+
+    # Premise check: the two runs must be genuinely different computations, or the equality
+    # below is vacuous. et[1] is decomposition-SENSITIVE on the unified driver, so its differing
+    # witnesses that the decompositions differed; a failure here means the chunk steering
+    # stopped taking effect (or et[1] invariance was restored) and wants a look. Its value is
+    # pinned by the regression harness, not here.
+    @test auto.free_boundary !== nothing && finer.free_boundary !== nothing
+    @test finer.free_boundary.et[1] != auto.free_boundary.et[1]
+
+    @test size(finer.delta_prime.matrix) == size(auto.delta_prime.matrix)
+
+    # Bit-identical, not approximate: a tolerance would hide exactly the reassociation drift
+    # this test exists to catch. The element-wise `===` pairs with the whole-matrix `==`
+    # because NaN === NaN — the `==` is what fails if the computation degrades to NaN.
+    for j in 1:size(auto.delta_prime.matrix, 1)
+        @test finer.delta_prime.matrix[j, j] === auto.delta_prime.matrix[j, j]
+    end
+    @test finer.delta_prime.matrix == auto.delta_prime.matrix
+end