Tearing - BUGFIX! - Convert Δ' to the r_s reference length before slab-layer matching

[d-burg]

Release note

• Audience: users
• Numerical impact: SLAYER tearing eigenvalues move by the reference-length conversion of Δ' — growth rates by up to 42% (max element, the outermost analyzed surface of the DIII-D-like case), rotation frequencies by ≤ 0.1%. All layer inputs (rs, shear, S, P_perp, tauk) are unchanged to roundoff, the raw Δ' matrix and the whole outer-region solve are bit-identical, and the GGJ path is unchanged. (harness @ 63dcc81), re-run against develop e097951 after the develop sync merge and a comment-only src/ cleanup moved the branch past the previous stamp
• Migration: no input changes. Re-baseline any stored SLAYER eigenvalues (omega_Hz, gamma_Hz): the dispersion now matches Δ' in the r_s reference. The HDF5 PerSurface group gains k_ref, mu_mercier, and delta_prime_conversion datasets documenting the transform, and the Delta_prime_matrix long_name notes the r_s reference.

SLAYER's dispersion relation compared the outer-region Δ' (referenced to unit ψ_N) against layer quantities referenced to the minor radius r_s, overstating or understating growth rates surface-by-surface. Δ' is now converted to the r_s reference length via the Mercier-exponent K^(2μ) transform before slab-layer matching, and the radial label defining r_s (five options, including Fitzpatrick's toroidal-flux label) is selectable programmatically, with the historical midplane default unchanged.

Regression report

regress --cases diiid_slayer_n1,diiid_n1_riccati --refs e09795119,63dcc8113 --force (develop tip vs branch head, raw shas so no stale local-branch resolution, both sides force-re-run so neither comes from the results cache; same pinned manifest, julia 1.11.6):

Regression Report: diiid_slayer_n1
Ref 1: e09795119  @ e09795119 (2026-08-27)
Ref 2: 63dcc8113  @ 63dcc8113 (2026-08-28)
Quantity                            e09795119  63dcc8113  Diff                Status
SLAYER surface indices              [6 elem]   [6 elem]   0.0e+00             OK
SLAYER poloidal m                   [6 elem]   [6 elem]   0.0e+00             OK
SLAYER toroidal n                   [6 elem]   [6 elem]   0.0e+00             OK
SLAYER minor radius rs              [6 elem]   [6 elem]   0.0e+00             OK
SLAYER r-based shear                [6 elem]   [6 elem]   6.009e-08 (0.00%)   ** CHANGED **
SLAYER Lundquist S                  [6 elem]   [6 elem]   6.967e-02 (0.00%)   ** CHANGED **
SLAYER D_norm                       [6 elem]   [6 elem]   3.127e-10 (0.00%)   ** CHANGED **
SLAYER P_perp                       [6 elem]   [6 elem]   0.0e+00             OK
SLAYER tauk                         [6 elem]   [6 elem]   5.5e-15             OK
SLAYER iota_e                       [6 elem]   [6 elem]   1.1e-16             OK
SLAYER Q_root [2/1,3/1,4/1]         [3 elem]   [3 elem]   8.889e-03 (1.27%)   ** CHANGED **
SLAYER ω_Hz [2/1,3/1,4/1]           [3 elem]   [3 elem]   1.143e+01 (0.10%)   ** CHANGED **
SLAYER γ_Hz [2/1,3/1,4/1]           [3 elem]   [3 elem]   3.317e+02 (42.14%)  ** CHANGED **
SLAYER no_root flags [2/1,3/1,4/1]  [3 elem]   [3 elem]   0.0e+00             OK
SLAYER enabled flag                 1          1          0.0e+00             OK
Runtime (s)                         130.4s     131.8s                         --
Summary: 6 changed, 9 unchanged

Regression Report: diiid_n1_riccati
Ref 1: e09795119  @ e09795119 (2026-08-27)
Ref 2: 63dcc8113  @ 63dcc8113 (2026-08-28)
Quantity                       e09795119      63dcc8113      Diff     Status
delta prime (BVP diagonal)     [5 elem]       [5 elem]       0.0e+00  OK
delta prime (raw side-major)   [10 elem]      [10 elem]      0.0e+00  OK
edge coil response delta_coil  [10 elem]      [10 elem]      0.0e+00  OK
total energy Re(et[1])         8.038607e-01   8.038607e-01   0.0e+00  OK
plasma energy Re(ep[1])        -1.344628e+00  -1.344628e+00  0.0e+00  OK
vacuum energy Re(ev[1])        2.148488e+00   2.148488e+00   0.0e+00  OK
total energy (all)             [35 elem]      [35 elem]      0.0e+00  OK
# singular surfaces            5              5              0.0e+00  OK
singular psi locations         [5 elem]       [5 elem]       0.0e+00  OK
singular q values              [5 elem]       [5 elem]       0.0e+00  OK
ODE steps (saved)              51             51             0.0e+00  OK
ODE steps (total)              1667           1667           0.0e+00  OK
mpert                          35             35             0.0e+00  OK
npert                          1              1              0.0e+00  OK
q0                             1.204212e+00   1.204212e+00   0.0e+00  OK
q95                            4.781723e+00   4.781723e+00   0.0e+00  OK
beta_n                         1.372511e+00   1.372511e+00   0.0e+00  OK
Runtime (s)                    113.5s         111.2s                  --
Summary: 17 unchanged

The moved quantities are exactly the intent: the dispersion roots move because Δ' now enters the layer matching in the r_s reference. The 1e-8–1e-9 relative movement in shear/S/D_norm is the truncation error of the finite-difference stencil this branch replaces with analytic interpolant derivatives, not a physics change. diiid_n1_riccati is bit-identical throughout, including the Δ' BVP diagonal and raw matrix — this PR converts Δ' downstream of the outer-region solve and does not perturb it. No tolerance or acceptance criterion was altered.

Summary

The SLAYER tearing dispersion relation Δ' − Δ_crit = S^{1/3}Δ̂(Q) mixed two reference-length conventions: the outer BVP Δ' is referenced to unit Δψ_N (STRIDE convention), while every term on the layer side — the Fitzpatrick slab Δ(Q), the Lundquist number S = (μ₀r_s²/η)/τ_H, and the rfitzp/lar/toroidal critical-Δ — is referenced to unit x̂ = (r−r_s)/r_s. This PR converts Δ' to the r_s reference at the matching point:

Δ̂_ij = K_i^(1/2+μ_i) · Δ'_ij · K_j^(μ_j−1/2),   K = r_s·(dψ_N/dr)|_s,  μ = √(−D_I)

whose diagonal is K^(2μ)·Δ'_kk, reducing to the textbook Δ̂ = r_s·Δ'_phys at the slab point D_I = −¼. K and μ are computed per surface from the equilibrium (surface_minor_radius/surface_da_dpsi, GGJ E+F+H−¼ from ResistGeometry).

Why this is correct

• Frobenius scaling: near a rational surface the tearing-parity outer solution is A_L|x|^(1/2−μ) + A_S|x|^(1/2+μ) (Glasser, Wang & Park 2016, Eq. 26); rescaling x_ψ = K·x̂ maps Â_L = A_L K^(1/2−μ), Â_S = A_S K^(1/2+μ), giving the transform above.
• The layer side is r_s-referenced end to end — verified against the TJ code's own definitions: its critical-Δ uses Ŵ_d = W_d/r_s (proven by its source multiplying the reported width by r_s to convert to a-units), its Δ carries r_s^(1+2ν) by definition (its docs' pressure-flattening derivation, 1+2ν = 2μ), and its dispersion assembly (Δ − Δ_c)/S^{1/3} is identical in form to ours. The dc_tmp formulas in LayerParameters.jl are a line-for-line port of that r_s-referenced formula set.
• Parameter-free empirical check: on the frozen TJ circular β/ε benchmark scans, K^(2μ) computed from TJ's own equilibrium files reproduces the measured Δ̂_TJ/Δ'_GPEC ratio to +1 % median over 23 pole-free points; in absolute values the 2/1 median |difference from TJ| drops 17.5 % → 2.5 % (β) and 17.6 % → 1.1 % (ε), inside the combined grid band. The factor is 0.84 at the 2/1 but 1.20 at the 3/1, and raw GPEC was correspondingly +18 % high / −28 % low — one formula corrects both directions.

What changes numerically

• SLAYER path only. γ, ω, and the effective Δ' move; on the shipped DIII-D-like SLAYER example the q=2 diagonal factor is K^(2μ) = 0.906 (~9 % stabilizing shift in Δ_eff; no sign change — the flip threshold there is 1.585).
• The BVP Δ' outputs do NOT change: SingularSurfaces/Delta_prime_matrix stays ψ_N-referenced; STRIDE benchmarking is unaffected. Tearing/PerSurface/Delta_prime_matrix now stores the r_s-referenced matrix actually used in the matching (long_name updated), plus new per-surface datasets k_ref, mu_mercier, delta_prime_conversion.
• GGJ path unchanged (genuinely toroidal/ψ-based; its rescale_delta handles inner→outer units natively).
• Hand-built SLAYERParameters default to k_ref = 1 → identity conversion, so existing unit tests and run_slayer_from_inputs workflows are unaffected unless the caller opts in.

Changes

• SLAYERParameters: new fields k_ref, mu_mercier (documented; default 1.0/0.5 = identity).
• build_slayer_inputs: fills both from the equilibrium; Mercier-unstable surfaces (D_I ≥ 0) get μ → 0 continuously (Δ' left raw) — documented in-line.
• run_slayer_from_inputs: applies delta_prime_to_rs_reference (new, documented) for the slab-layer path.
• surface_coupling docstring corrected: S^{1/3} maps the inner Δ to the r_s-referenced outer convention, not "ψ-units" — the previous wording was the error that concealed this mismatch.
• HDF5: three new Tearing/PerSurface datasets + annotations; SLAYERResult.dp_matrix docstring notes the convention.
• Tests: transform verified (diagonal K^(2μ), off-diagonal split factors, slab-point reduction to K, identity default).

Validation

• SLAYER/inner-layer test files pass (runtests_slayer_params / _runner / _inputs / _riccati, runtests_innerlayer — 188 tests, includes new transform testset)
• Regression harness develop vs this branch: only the SLAYER dispersion roots moved (γ up to 42.14 % max-element; ω 0.10 %); all layer inputs and all 17 Riccati/BVP quantities (incl. SingularSurfaces/Delta_prime_matrix) bit-identical — current tables above; earlier full report

Benchmark figures and the full theory audit (TJ source citations, ratio test, corrected Δ' comparison plots) are recorded in the maintainer's private benchmark repo: issue · figures + scripts (maintainer access).

Important

This changes every SLAYER γ (≈9 % shift in Δ_eff at q=2 on the shipped example). Do not merge without human review of the derivation and the regression report.

🤖 Generated with Claude Code

[d-burg]

Regression harness report — develop (9491f89) vs this branch (af67255)

Cases: diiid_slayer_n1 (SLAYER γ — expected to move) and diiid_n1_riccati (BVP Δ' — must not move). Same pinned manifest (7e5c34ad) on both refs.

Signature is exactly as intended: the only changed quantities are the SLAYER dispersion roots (Q_root 1.2%, ω 0.12%, γ up to 42% on the max element — the ~9% stabilizing Δ_eff shift amplified through the layer response); every layer input (rs, shear, S, D_norm, P_perp, tauk, iota_e) is bit-identical, and all 17 Riccati/BVP quantities — including the Δ' diagonal and raw side-major matrix — are unchanged at 0.0e+00.

Full report

================================================================
Case: diiid_slayer_n1 — DIII-D-like H-mode equilibrium, n=1, SLAYER tearing-mode analysis (uncoupled per-surface, AMR, validity-gated)
================================================================

Regression Report: diiid_slayer_n1
======================================================================================================================
Ref 1: develop  @ 9491f893 (2026-08-16)
       env: julia 1.11.6, arm64-apple-darwin24.0.0, manifest 7e5c34ad (pinned), 8 threads/8 BLAS
Ref 2: bugfix/slayer-dprime-reference-length  @ af672554 (2026-08-17)
       env: julia 1.11.6, arm64-apple-darwin24.0.0, manifest 7e5c34ad (pinned), 8 threads/8 BLAS
----------------------------------------------------------------------------------------------------------------------
Quantity                            develop   bugfix/slayer-dprime-reference-length  Diff                Status       
----------------------------------------------------------------------------------------------------------------------
SLAYER surface indices              [6 elem]  [6 elem]                               0.0e+00             OK           
SLAYER poloidal m                   [6 elem]  [6 elem]                               0.0e+00             OK           
SLAYER toroidal n                   [6 elem]  [6 elem]                               0.0e+00             OK           
SLAYER minor radius rs              [6 elem]  [6 elem]                               0.0e+00             OK           
SLAYER r-based shear                [6 elem]  [6 elem]                               0.0e+00             OK           
SLAYER Lundquist S                  [6 elem]  [6 elem]                               0.0e+00             OK           
SLAYER D_norm                       [6 elem]  [6 elem]                               0.0e+00             OK           
SLAYER P_perp                       [6 elem]  [6 elem]                               0.0e+00             OK           
SLAYER tauk                         [6 elem]  [6 elem]                               0.0e+00             OK           
SLAYER iota_e                       [6 elem]  [6 elem]                               0.0e+00             OK           
SLAYER Q_root [2/1,3/1,4/1]         [3 elem]  [3 elem]                               1.328e-02 (1.22%)   ** CHANGED **
SLAYER ω_Hz [2/1,3/1,4/1]           [3 elem]  [3 elem]                               1.359e+01 (0.12%)   ** CHANGED **
SLAYER γ_Hz [2/1,3/1,4/1]           [3 elem]  [3 elem]                               3.182e+02 (42.21%)  ** CHANGED **
SLAYER no_root flags [2/1,3/1,4/1]  [3 elem]  [3 elem]                               0.0e+00             OK           
SLAYER enabled flag                 1         1                                      0.0e+00             OK           
Runtime (s)                         119.7s    118.9s                                                     --           
======================================================================================================================
Summary: 3 changed, 12 unchanged


================================================================
Case: diiid_n1_riccati — DIII-D-like equilibrium, n=1, Riccati integrator Δ' matrix
================================================================

Regression Report: diiid_n1_riccati
====================================================================================================
Ref 1: develop  @ 9491f893 (2026-08-16)
       env: julia 1.11.6, arm64-apple-darwin24.0.0, manifest 7e5c34ad (pinned), 8 threads/8 BLAS
Ref 2: bugfix/slayer-dprime-reference-length  @ af672554 (2026-08-17)
       env: julia 1.11.6, arm64-apple-darwin24.0.0, manifest 7e5c34ad (pinned), 8 threads/8 BLAS
----------------------------------------------------------------------------------------------------
Quantity                       develop        bugfix/slayer-dprime-reference-length  Diff     Status
----------------------------------------------------------------------------------------------------
delta prime (BVP diagonal)     [5 elem]       [5 elem]                               0.0e+00  OK    
delta prime (raw side-major)   [10 elem]      [10 elem]                              0.0e+00  OK    
edge coil response delta_coil  [10 elem]      [10 elem]                              0.0e+00  OK    
total energy Re(et[1])         8.038607e-01   8.038607e-01                           0.0e+00  OK    
plasma energy Re(ep[1])        -1.344628e+00  -1.344628e+00                          0.0e+00  OK    
vacuum energy Re(ev[1])        2.148488e+00   2.148488e+00                           0.0e+00  OK    
total energy (all)             [35 elem]      [35 elem]                              0.0e+00  OK    
# singular surfaces            5              5                                      0.0e+00  OK    
singular psi locations         [5 elem]       [5 elem]                               0.0e+00  OK    
singular q values              [5 elem]       [5 elem]                               0.0e+00  OK    
ODE steps (saved)              51             51                                     0.0e+00  OK    
ODE steps (total)              1667           1667                                   0.0e+00  OK    
mpert                          35             35                                     0.0e+00  OK    
npert                          1              1                                      0.0e+00  OK    
q0                             1.204212e+00   1.204212e+00                           0.0e+00  OK    
q95                            4.781723e+00   4.781723e+00                           0.0e+00  OK    
beta_n                         1.372511e+00   1.372511e+00                           0.0e+00  OK    
Runtime (s)                    99.6s          98.3s                                           --    
====================================================================================================
Summary: 17 unchanged

🤖 Generated with Claude Code

[d-burg]

Shaped-case check: shipped DIII-D-like SLAYER example (auto grid, rfitzp, n=1)

End-to-end run on this branch; ratio is Tearing/PerSurface/Delta_prime_matrix diagonal over the raw SingularSurfaces/Delta_prime_matrix diagonal, compared against the recorded delta_prime_conversion dataset. All match K^(2μ) to 1e-10.

q	k_ref	μ	K^(2μ)	Δ' raw	Δ' converted
2	0.9125	0.540	0.9059	9.0969	8.2405
3	1.2245	0.529	1.2392	−6.1961	−7.6782
4	1.3597	0.558	1.4092	−16.375	−23.076
5	1.5051	0.831	1.9735	−2475.3	−4885.0
6	1.5502	0.521	1.5787	−194.65	−307.30
7	1.5541	0.504	1.5599	242.23	377.86

Notes for review:

• q=2 (the 2/1): −9.4 % in Δ'; with Δ_crit = 14.42 the effective drive moves −5.32 → −6.18. No sign changes on any surface (the q=2 flip threshold is 1.585).
• Shaping puts K > 1 outboard of q≈2.7, so the conversion amplifies |Δ'| at q ≥ 3 — this is why γ moved harder on outer surfaces in the regression report.
• q=5 has μ = 0.83 (strongly non-slab, D_I = −0.69); its ×1.97 factor is consistent bookkeeping, but the slab-layer match itself is a poor approximation at that surface regardless of convention.

🤖 Generated with Claude Code

[d-burg]

Benchmark figures: Δ' alignment with TJ, before and after the K^(2μ) conversion

Frozen TJ circular benchmark scans (q_edge = 3.6): TJ = hollow black rings (its own r_s-referenced Δ̂), julia_GPEC raw = dashed line (ψ_N-referenced BVP Δ', the pre-PR pairing), julia_GPEC ×K^(2μ) and Fortran STRIDE ×K^(2μ) = filled markers. K and μ are computed per point/surface from TJ's own equilibrium files — nothing fitted. Grey vertical band = pole region (each code approaches its own ideal-kink limit there; TJ's β-limit sits at p_f = 0.1675 vs 0.1851 for both STRIDE codes, so the comparison is not meaningful past the band edge).

Same data as percentage difference from TJ (grey horizontal band = ±3 % combined grid uncertainty of the two STRIDE codes):

Reading guide:

• 2/1 (left panels): median |difference from TJ| drops 17.5 % → 2.5 % (β scan) and 17.6 % → 1.1 % (ε scan) — into the grid band, i.e. the residual is no longer resolvable against either code's own grid noise.
• 3/1 (right panels): the conversion factor is > 1 here (1.20 vs 0.84 at the 2/1) and the raw error had the opposite sign (−28 % low vs +18 % high) — one parameter-free formula corrects both directions, which a fitted constant cannot do. It improves (≈22 % → ≈9 %) without fully closing: the β-3/1 window is contaminated by that surface's own pole, and the ε-3/1 retains an unexplained −15 % → 0 drift at low ε (possibly the 3/1 sitting at r̂_s ≈ 0.909 against the q_edge truncation — untested).
• Fortran STRIDE (orange) collapses identically to julia_GPEC, as it must — both share the ψ_N convention; this is a convention fix, not a solver change.

Figures live on the non-merging assets/dprime-convention-figures branch (not part of this PR's diff); generating scripts are in the maintainer's benchmark repo linked in the PR body.

🤖 Generated with Claude Code

[d-burg]

Same comparison on fixed grids: julia ldp mpsi=256 (vs Fortran ldp mpsi=128)

The figures in the previous comment used julia's auto@1e-3 grid, which re-tessellates knots per equilibrium (npsi varies 734–1003 across neighboring scan points) and therefore carries ~1 % point-to-point realization jitter — visible as scatter relative to the smoother Fortran fixed-grid series. Re-plotting with julia's fixed ldp mpsi=256 scan removes that jitter without changing any conclusion; roughness metric (median |second difference| of the ratio-to-TJ, 2/1 pole-free windows):

series	β roughness (med/max)	ε roughness (med/max)
julia auto@1e-3	1.38 % / 7.94 %	0.91 % / 1.89 %
julia ldp@256	0.71 % / 2.45 %	0.16 % / 0.60 %
fortran ldp@128	0.34 % / 1.75 %	0.20 % / 1.57 %
fortran ldp@256	3.64 % / 6.83 %	2.12 % / 10.24 %

(The noise is a property of the auto grid's per-equilibrium adaptivity, not of the solver — julia's own fixed-grid series is the smoothest of the four. Fortran mpsi=256 is excluded from the figures because ldp knot oversampling above mpsi≈256 injects solver noise.)

With the fixed grid, the corrected 2/1 medians are 1.4 % (β) and 1.4 % (ε) from TJ — spreads [−6.1 %, +0.2 %] and [−3.7 %, −1.0 %] — and the two STRIDE implementations are nearly indistinguishable point-for-point after the conversion. Median offsets are grid-robust (16.7–18.5 % raw across all four grid choices), so the K^(2μ) result is unchanged; only the scatter tightens.

Generated with JULIA_SRC=ldp via the same scripts (see PR body links).

🤖 Generated with Claude Code

[d-burg]

Cross-reference: the "unexplained low-ε 3/1 drift" flagged in this PR's benchmark comments is now root-caused and fixed — it was a Riccati shooting error-control bug (#412, draft PR #413), not toroidal physics. With #413's branch the ε-scan 3/1 median |difference from TJ| after the K^(2μ) conversion drops 5.7 % → 2.3 % (see the scan-scale validation on #413), while the 2/1 panels here are unaffected (Δ'_21 moves ≤ 0.11 % across all 196 scan points). The two PRs are independent; merging #413 first and rebasing this branch is the suggested order.

🤖 Generated with Claude Code

[d-burg]

Rebased onto develop (c42d558e, includes #399 physical-B_T and #339 multi-ion NTV) + radial-label options added

Branch is now two commits on the develop tip:

1. The K^(2μ) reference-length conversion (unchanged; one trivial docstring conflict resolved against the new SLAYERResult fields).
2. New: rs_method options :halfwidth (midplane half-chord (R_out−R_in)/2 — shift-free stand-in for the circular-theory flux label; reproduces it to ~1% on circular benchmarks) and :volume (cylinder-equivalent √(V/2π²R₀), Rutherford-literature convention), alongside :midplane (default, unchanged) and :fsa. Every label drives S, the r-based shear, W_d, and k_ref together, so each choice is self-consistent by construction. Programmatic API only — deliberately not exposed via TOML (the label moves γ by O(1) at near-edge surfaces of shaped equilibria, a specialist decision; defaults and all TOML-driven results are unchanged, which the regression below confirms).

Tests on the rebased branch: LayerInputs 50/50 (includes new label testsets), Runner 79/79, LayerParameters 49/49, slayer-riccati pass.

Regression harness, develop c42d558e vs this branch (8315b5a1):

• diiid_slayer_n1: only the dispersion roots moved — γ up to 42.2 % max-element (the K^(2μ) Δ_eff shift through the layer response, same magnitude as the pre-rebase A/B, now on top of the physical-B_T fix), ω 0.10 %; every layer input (rs, shear, S, D_norm, P_perp, tauk, iota_e) bit-identical — the label commit is provably inert at defaults.
• diiid_n1_riccati: 17/17 unchanged including the BVP Δ' diagonal and raw matrix.

Radial-label sensitivity evidence (four-arm A/B on the shipped shaped example + circular-TJ label validation) is recorded in the maintainer's benchmark repo issue linked in the PR body.

🤖 Generated with Claude Code

[d-burg]

Radial-label unification (commit 06be788)

A parallel investigation into a resistive-layer-overlap truncation criterion (psihigh cap, PR to follow; #298 context) independently ran into the same radial-label question this PR addresses, and the two efforts collided in LayerInputs.jl — both rewrote the rs_method closures. This commit unifies the shared infrastructure here, so the cap PR can stack on this branch. What changed:

• One exported radial_label(equil; rs_method, theta) returning (r, dr/dpsi) closures, replacing the inline closures in build_slayer_inputs. Five labels: :midplane (default, unchanged), :halfwidth, :fsa, :volume, and the new :flux — Fitzpatrick's toroidal-flux surface label (Nucl. Fusion 2025, Eq. 30), r = sqrt(2 psi_t / B0) with the g = F/(B0 R0) correction carried (g departs from 1 by ~3 percent on the DIII-D-like deck).
• Analytic derivatives everywhere. The finite-difference stencils are gone; da/dpsi comes from the interpolants' own psi-derivatives ((∂r²/∂ψ)/(2a) per theta; dV/dψ/(4π²R₀r) for :volume; dr/dψ ∝ q for :flux). The old stencil silently clamped within 1e-4 of the flux boundaries — the limiting factor on how close to the separatrix a width could be evaluated. Agreement with the stencil is ~1e-10 away from the boundaries, so default results move only at roundoff (regression refresh below).
• k_ref and the K^(2mu) Delta-prime conversion consume the same closures, so any label — including :flux — drives S, the r-based shear, W_d, and the Delta-prime reference length together with no extra plumbing.
• :fsa gains the max(r^2, 0) guard needed on extrapolated surfaces (it previously threw DomainError past the boundary).

Findings on record from the cross-validation of the two investigations

1. Two independent lines of evidence agree that :midplane is the odd label out for Fitzpatrick-formalism quantities: this PR's TJ benchmark (flux-type labels reproduce TJ's reference to <= 1.2 percent on the circular eps = 0.15 case, midplane 1.9 percent off) and the overlap work's non-covariance argument (Eq. 100 scales as J^(1/2) under relabeling, so it is anchored to the Eq. 30 flux label; the label choice moves the overlap cut by one rational surface).
2. Scope of the circular validation: on circular equilibria all flux-type labels coincide to O(eps^2), so that test certifies the family but cannot discriminate within it. On the shaped deck, per-surface r values show :fsa tracks :flux to ~1 percent at psi = 0.5 (degrading to ~12 percent by psi = 0.99) while :halfwidth sits ~16 percent from :flux in the core — :halfwidth is its own convention on shaped equilibria, not an Eq. (30) proxy. Docstrings updated to say so.
3. The "which label should be the SLAYER-wide default" question is deliberately out of scope for this PR and the cap PR — it is a results-moving !-tagged change of its own, tracked in InnerLayer.SLAYER - API - Choose one default radial label for the Fitzpatrick layer formalism #417, which now frames the edge-criterion label (settled: Eq. 30 flux) separately from the interior conversion default (bounded convention choice, ~ +/- 3 percent in the 2/1 growth rate at defaults).

Tests: LayerInputs 56 (incl. a five-label analytic-vs-FD derivative testset), Runner 79/79, params green. Regression harness refresh of 8315b5a -> 06be788 on diiid_slayer_n1 is running; table to follow in a separate comment.

---

NO MERGE WITHOUT HUMAN REVIEW — this PR remains a draft and requires a named human reviewer and assignee before it can be considered for merge. This is non-negotiable.

[d-burg]

Regression refresh: 8315b5a -> 06be788 (diiid_slayer_n1)

Both refs run as git worktrees under the same pinned manifest (f72f3b3b), julia 1.11.6, 8 threads.

Quantity	Diff	Status
surface indices, m, n, minor radius rs	0.0e+00	OK
r-based shear	6.0e-08 (0.00%)	changed
Lundquist S	7.0e-02 abs (~1e-9 rel)	changed
D_norm	3.1e-10 (0.00%)	changed
P_perp, tauk, iota_e	0 / 5.5e-15 / 1.1e-16	OK
Q_root [2/1,3/1,4/1]	5.8e-06	OK
omega_Hz	1.9e-01	OK
gamma_Hz	1.5e-01 (0.01% max element)	changed
no_root flags, enabled	0.0e+00	OK

Interpretation: rs is bit-identical (the default :midplane radius formula is untouched); the 1e-8-1e-9 relative movement in shear/S/D_norm is precisely the h^2 ~ 1e-10 truncation error of the removed finite-difference stencil, i.e. the analytic derivative correcting the FD approximation, and it amplifies to at most 0.01 percent in the extracted growth rates through root-finding. The refactor is inert at defaults to within the old stencil's own numerical noise. No tolerance or acceptance criterion was changed anywhere in this commit.

Post-hoc annotation: the omega_Hz (1.9e-01) and gamma_Hz (0.01%) lines in this table were later found to be contaminated by a harness cached-sample artifact (see the forced-rerun verification comment below) and are upper bounds, not measured effects; the true effect of the derivative swap is at or below the ~1e-5 Hz sample-noise ceiling. The interpretation paragraph's conclusion (inert at defaults) stands, conservatively.

[d-burg]

Forced-rerun verification of the regression stamp. A harness caching subtlety was found today: the results DB keeps one cached sample per commit (UNIQUE(commit_hash, case_name), database.jl), and the threaded SLAYER root search is sample-dependent at the ~1e-5 Hz level, so a comparison against a cached row can inherit that row's jitter as an apparent change. (Diagnosis from this PR's session by cross-comparing two unrelated diffs sharing one cached row; confirmed by the layer-overlap session's forced experiment, which collapsed a nil-mechanism comparison's residuals 5,000–17,000x. tracked in #418.)

To rule it out for this PR, the stamped c42d558e → 06be78884 comparison was re-run with --force (both sides fresh, cache bypassed): it reproduces the body's table — Q_root 1.27%, omega 0.10%, gamma 42.17% (vs 42.16% stamped), layer-input lines identical. These movements are mechanism (the K^(2mu) conversion feeding the dispersion), not cache artifacts.

One earlier comment in this thread needs a precision note: the intermediate 8315b5a1 → 06be78884 comparison (analytic-derivative swap in isolation) quoted residuals of omega 1.9e-01 Hz / gamma 0.01%. That comparison diffed against the cached 06be78884 sample and its sub-percent residuals are therefore upper bounds contaminated by the cache artifact, not measured effects; the true effect of the derivative swap is at or below the ~1e-5 Hz sample-noise ceiling measured with forced runs. The conclusion stated there ("inert at defaults") stands, conservatively. The refresh comment above has been annotated accordingly.

[d-burg]

Closing summary: per-surface effect on the reference DIII-D-like equilibrium

Direct base-vs-head runs of examples/DIIID-like_SLAYER_example (n=1, uncoupled, AMR, midplane label, dc_type = none), base = develop 8f4be5fe, head = baf674659. The regression harness tracks only the 2/1–4/1 max-element diffs; this table is the full per-surface picture.

Surface	ψ_N	K	μ	K^(2μ)	Δ' (ψ_N ref)	Δ' (r_s ref)	γ: base → head	ω: base → head
2/1	0.518	0.913	0.540	0.906	+9.097	+8.240	+228.3 → +207.0 Hz (−9.3 %)	−10525.7 → −10526.8 Hz (+0.01 %)
3/1	0.770	1.225	0.529	1.239	−6.196	−7.678	−202.4 → −251.5 Hz (−24 %)	−9673.5 → −9674.5 Hz (+0.01 %)
4/1	0.893	1.360	0.558	1.409	−16.37	−23.07	−787.1 → −1118.8 Hz (−42 %)	−11086.6 → −11098.0 Hz (+0.10 %)
5/1	0.968	1.505	0.831	1.973	−2475	−4885	+30010 Hz → no valid root	−62703 Hz → —
6/1	0.993	1.550	0.521	1.579	−194.6	−307.2	−160885 Hz → no valid root	−52925 Hz → —
7/1	0.9985	1.554	0.504	1.560	+242.3	+378.0	+336507 → +471980 Hz (+40 %)	+34746 → +51226 Hz (+47 %)

What to read from it

• The physical headline is the 2/1 becoming less unstable (γ −9.3 %). Inside ψ_N ≈ 0.6 the flux label grows faster than the minor radius, so K < 1 and K^(2μ) < 1 reduces the tearing drive; everything outside has K^(2μ) > 1 and its already-stable Δ' is made more strongly stabilizing (3/1, 4/1 more damped).
• ω is untouched on the interior surfaces (≤ 0.1 %): the rotation frequency is set by ω_* and ω_E, not Δ', so the conversion affects growth, not rotation. The harness's "ω 0.10 %" line is the 4/1.
• The harness's "γ 42.16 % max element" is the 4/1 (−787 → −1119 Hz), a stable surface getting more stable — not the 2/1.
• A change the harness does not track: the 5/1 and 6/1 flip from "root found" to "no valid root." Before conversion the 5/1 reported γ = +30 kHz on a surface with Δ' = −2475 — an unphysical root the validity gate exists to reject (the example's own TOML comment states the 5/1 is expected to be dropped). With |Δ'_rs| roughly doubled (the 5/1 has the largest Mercier index on the table, μ = 0.83, hence the largest K^(2μ) = 1.97), the AMR search finds no admissible root on either surface. This is the documented intended behavior, but it is a visible change in the per-surface output that reviewers should know about.
• Surfaces at ψ_N > 0.99 (6/1, 7/1) are in the edge region where the radial label is ambiguous and the layer-overlap criterion applies; shown for completeness, not interpretation.
• Unchanged, by construction: every layer input (rs, shear, S, P_perp, tauk, iota_e), the ψ_N-referenced BVP matrix (SingularSurfaces/Delta_prime_matrix, retained verbatim), the whole outer-region solve (diiid_n1_riccati 17/17 bit-identical), and the GGJ path.

Other changes in this PR worth noting for review

1. PerSurface/k_ref, mu_mercier, delta_prime_conversion are new HDF5 datasets; PerSurface/Delta_prime_matrix is now the converted (r_s-referenced) matrix and its long_name says so.
2. The conversion is applied unconditionally on the SLAYER path (if !_is_ggj(model)). There is deliberately no off-switch: the unconverted matching was a convention mismatch, not an alternative mode. The pre-conversion matrix remains recoverable from SingularSurfaces/Delta_prime_matrix.
3. rs_method (radial label: :midplane default, :halfwidth, :fsa, :volume, :flux) is programmatic only — not exposed via TOML. The default label is unchanged, so the numbers above are the midplane-label results; label sensitivity is ±3 % on the 2/1 and O(1) at q ≳ 4 (InnerLayer.SLAYER - API - Choose one default radial label for the Fitzpatrick layer formalism #417).
4. The finite-difference da/dψ stencil is replaced by analytic interpolant derivatives for all labels; default results move at the stencil's own 1e-8–1e-9 truncation level (visible as the shear/S/D_norm lines in the harness table).

---

NO MERGE WITHOUT HUMAN REVIEW — this PR requires approval from its named reviewers (logan-nc, matt-pharr). Non-negotiable.

[logan-nc]

@d-burg the requested reviewers are not available → assigning @ebursch. Please post confirmation that this is ready to review and provide a synthesis/guide for him as this is possibly an overwhelming wall of text in the PR comments rn.

An artifact with workflows and figures, as per the claude.md instructions I shared earlier, would probably be quite helpful here.

@ebursch - you got this!

[d-burg]

Rebased onto develop d7e6e6fe; regression re-stamped on the current tip

The previous stamp was taken against 8f4be5fe, now 25 commits behind — including #419 (Keep the sign of crit), the #383 FourFitVars → MatrixSplines API break, and the #400 ForceFreeStates reorg. A !-tagged PR whose whole claim is "these six quantities moved and nothing else did" needs its evidence taken on top of those.

Rebase was textually clean — this PR touches only InnerLayer/ and Tearing/, which none of those 25 commits touch. Verified semantically too: the package loads, and the SLAYER/InnerLayer suites pass 273/273 (runtests_innerlayer, _slayer_params, _slayer_riccati, _slayer_inputs, _slayer_runner).

Regression, both sides forced (cache bypassed), pinned manifest 638da0a2, julia 1.11.6:

Quantity	previous stamp 8f4be5fe→baf67465	now d7e6e6fe→29df8700e
SLAYER r-based shear	6.009e-08	6.009e-08
SLAYER Lundquist S	6.967e-02	6.967e-02
SLAYER D_norm	3.127e-10	3.127e-10
SLAYER Q_root	8.892e-03 (1.27%)	8.889e-03 (1.27%)
SLAYER ω_Hz	1.143e+01 (0.10%)	1.143e+01 (0.10%)
SLAYER γ_Hz	3.318e+02 (42.16%)	3.317e+02 (42.14%)
diiid_n1_riccati	17/17 unchanged	17/17 unchanged

The signature is unchanged, which is the point: diiid_n1_riccati is still bit-identical across all 17 quantities including the BVP Δ' diagonal and the raw side-major matrix, so #419 and the MatrixSplines break do not interact with this PR. The residual γ wobble (42.16 → 42.14 %) is the ~1e-5 Hz root-search sample noise tracked in #418.

One caveat on the record: the manifest here (638da0a2) is a fresh resolve on this machine, not the 7e5c34ad of the original stamp. Both sides share it so the A/B is valid, but it is a different package set than the earlier runs in this thread.

One code change in the rebase (29df8700e): the k_ref guard tested only for an exact zero. A non-finite da/dψ — reachable where the interpolant's r² is clamped non-positive on extrapolated surfaces — would have silently zeroed the converted Δ' diagonal, and a negative one would raise DomainError from the non-integer exponent. It now falls back to the identity conversion with a warning. Provably inert on this deck: rs is bit-identical and the dispersion diffs above match the old stamp to three significant figures.

Also corrected the Mercier-unstable comment. Clamping μ to 0 leaves only the diagonal raw — the off-diagonals still carry K_i^(1/2)·K_j^(−1/2) — and since the exponents are 1/2 ± i|μ| there, it discards a genuine phase rather than taking a continuous limit.

Previous head baf674659 is preserved locally if anyone needs to diff against it.

[d-burg]

Jacobian audit of the radial labels: :fsa and :halfwidth are coordinate-dependent

Prompted by @logan-nc's question on #298 ("easily testable by running with Hamada and PEST"), I ran every label in radial_label through three jacobians on the Solovev deck. The jacobian fixes the poloidal-angle coordinate, not the flux label, so every quantity evaluated at fixed ψ_N must be invariant. Max relative spread across pest / hamada / equal_arc:

label	spread in r_s	spread in dr/dψ	verdict
:midplane (default)	≤ 1.3e-16	≤ 3.0e-14	exactly invariant
:volume	≤ 6.7e-10	≤ 9.8e-09	invariant (quadrature noise)
:flux	≤ 6.9e-10	≤ 9.9e-09	invariant (quadrature noise)
:halfwidth	1.5e-4 … 3.6e-4	1.4e-2 … 3.8e-2	coordinate-dependent
:fsa	3.8e-2 … 6.6e-2	5.9e-2 … 1.45e-1	coordinate-dependent

The cause is the same in both: they average over the coordinate's θ with a uniform measure. :fsa is a 128-point uniform mean, which is a true flux-surface average only if the jacobian is already uniform in θ — so it silently inherits whatever θ-packing the coordinate system has, and moves r_s by up to 6.6 % and dr/dψ by up to 14.5 % on nothing but a jac_type change. :halfwidth assumes θ = 0.5 is the inboard midplane, which is not true in every coordinate.

The three labels defined from flux alone — :midplane (a single θ = 0 chord), :volume, :flux — are clean, and the default is unaffected. But the name :fsa claims a flux-surface average it does not compute, and its docstring currently advertises it as "the closest geometric approximation to :flux at interior surfaces" with no warning that the approximation depends on jac_type.

Suggested resolution (not in this PR unless reviewers want it here): either weight the mean by the jacobian to make it a genuine FSA, or rename it and document both it and :halfwidth as coordinate-dependent conventions not to be used for Fitzpatrick-formalism quantities.

[d-burg]

@logan-nc ready to review — rebased onto the develop tip and re-stamped the regression on it (comment above); the previous evidence predated #419 and #383.

@ebursch — a reading order so you can skip the archaeology.

What the bug is. The dispersion relation Δ' − Δ_crit = S^(1/3)Δ̂(Q) was comparing two different reference lengths. The outer BVP Δ' is referenced to unit Δψ_N (the STRIDE convention); every term on the layer side — the Fitzpatrick slab Δ(Q), S = (μ₀r_s²/η)/τ_H, and the critical-Δ — is referenced to unit x̂ = (r−r_s)/r_s. They were subtracted directly.

The fix, in one line. Δ̂_ij = K_i^(1/2+μ_i) · Δ'_ij · K_j^(μ_j−1/2) with K = r_s·(dψ_N/dr)|_s and μ = √(−D_I); the diagonal is K^(2μ)·Δ'_kk, and at the slab point D_I = −1/4 it reduces to the textbook Δ̂ = r_s·Δ'_phys. It follows from Frobenius scaling: the outer solution is A_L|x|^(1/2−μ) + A_S|x|^(1/2+μ), so rescaling x_ψ = K·x̂ maps Â_L = A_L K^(1/2−μ) and Â_S = A_S K^(1/2+μ). Note any constant in GPEC's Δ' definition cancels in the ratio, so the result does not depend on that convention.

The four things worth your scrutiny, in order:

1. The derivation above — is the ψ_N ↔ r_s identification of K right?
2. The empirical check: K^(2μ) computed from TJ's own equilibrium files reproduces the measured Δ̂_TJ/Δ'_GPEC ratio to +1 % median over 23 pole-free points, and the factor is 0.84 at the 2/1 but 1.20 at the 3/1 while raw GPEC was +18 % high / −28 % low. One parameter-free formula correcting both directions is the part I would try hardest to break.
3. The per-surface table in the closing-summary comment — in particular that the 5/1 and 6/1 flip from "root found" to "no valid root". Intended (the 5/1 was reporting γ = +30 kHz on a surface with Δ' = −2475), but it is a visible output change the harness does not track.
4. Whether SingularSurfaces/Delta_prime_matrix (ψ_N) and Tearing/PerSurface/Delta_prime_matrix (r_s) sharing a dataset name is acceptable. long_name documents it; anyone keying a post-processing script on the name gets different numbers without warning.

What you can skip: everything about the radial-label options (:halfwidth, :volume, :flux). They are programmatic-only, not exposed via TOML, and the regression confirms they are inert at defaults.

Caveat you should know going in: the TJ benchmark that carries claim (2) lives in a private repo, so it is not reproducible from this PR. The in-repo tests verify the transform's algebra against its own formula and would pass identically if the exponents were physically wrong. If that is not good enough, say so and we will land a reproducible artifact.

[d-burg]

@ebursch — per @logan-nc's suggestion, a self-contained synthesis of this PR so you can skip the thread archaeology: the reference-length mismatch as a picture, the K^(2μ) derivation, the TJ benchmark before/after, the per-surface effect on the shipped deck, and the radial-label audit, ending in a checklist of what to verify independently.

https://claude.ai/code/artifact/af31c362-8129-473b-824e-94d3df982559