InnerLayer.SLAYER - FEATURE! - 🚨 Derive the toroidal critical-Δ geometric factor in the r_s reference

[d-burg]

Release note

• Audience: users
• Numerical impact: none at defaults (dc_type = "none"); with dc_type = "toroidal" the critical-Δ offset is now computed from the equilibrium rather than from a user-supplied dgeo_val, and its parallel-conduction closure uses the same toroidal geometry (Δ_crit +3.7 to +5.9 % on the DIII-D-like example). The HDF5 dataset Tearing/PerSurface/D_geo carries derived values instead of zeros. (harness @ 90c3201)
• Migration: runs at the default dc_type = "none" are unaffected. Under dc_type = "toroidal", drop dgeo_val to take the derived geometric factor; dc_type = "toroidal" no longer errors when it is absent. A prescribed dgeo_val now means the factor in the r_s reference, i.e. Eq. 59's own value × k_ref·v1/V_s (≈ 2 at large aspect ratio). A value that was passed as Eq. 59's own factor must be converted, or it gives about half the intended Δ_crit. Anything reading Tearing/PerSurface/D_geo must now expect derived values where it previously read zeros. :lar and :rfitzp are unchanged.

The toroidal critical-Δ of Connor, Ham, Hastie & Liu 2015 (PPCF 57 065001, Eq. 59) can now be used on any equilibrium: its geometric factor is derived per rational surface from the flux-surface averages GPEC already computes, expressed in the same r_s reference length as the rest of the SLAYER layer stack.

Left: the critical Δ at each rational surface of the DIII-D-like example under the cylindrical rfitzp model, under the toroidal geometric factor with a cylindrical parallel-conduction closure, and under the full toroidal model. The toroidal model is 22 to 60 % above rfitzp; the closure accounts for 3.7 to 5.9 % of that. Right: SLAYER growth rates with no critical Δ (the default) and under each model. The 2/1 goes from growing (192 s⁻¹) to damped: −153 s⁻¹ with rfitzp, −230 s⁻¹ with the toroidal model.

Regression report

regress --cases diiid_slayer_n1,diiid_slayer_n1_toroidal --refs 6fc89d021,90c320185 --force on feynman (SLURM, 8 threads; current develop vs this branch with develop merged; both fresh, manifest pinned).

Regression Report: diiid_slayer_n1
==============================================================================================
Ref 1: 6fc89d021  @ 6fc89d021 (2026-10-09)
       env: julia 1.11.9, x86_64-linux-gnu, manifest f826ee26 (pinned), 8 threads/8 BLAS, jit build
Ref 2: 90c320185  @ 90c320185 (2026-10-09)
       env: julia 1.11.9, x86_64-linux-gnu, manifest f826ee26 (pinned), 8 threads/8 BLAS, jit build
----------------------------------------------------------------------------------------------
Quantity                                 6fc89d021  90c320185  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_e                               [6 elem]   [6 elem]   0.0e+00           OK
SLAYER Q_i                               [6 elem]   [6 elem]   0.0e+00           OK
SLAYER toroidal critical-Δ factor D_geo  [6 elem]   [6 elem]   8.151e+00 (Inf%)  ** CHANGED **
SLAYER Q_root [2/1,3/1,4/1]              [3 elem]   [3 elem]   1.7e-06           OK
SLAYER ω_Hz [2/1,3/1,4/1]                [3 elem]   [3 elem]   6.1e-02           OK
SLAYER γ_Hz [2/1,3/1,4/1]                [3 elem]   [3 elem]   8.8e-03           OK
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)                              280.9s     288.3s                       --
==============================================================================================
Summary: 1 changed, 17 unchanged

• D_geo: changes from zeros to the derived factor, which is the intended change.
• Roots: unchanged within the harness floors in this run. Neither side seeds the AMR triangulation (that fix is in Tearing - BUGFIX! - 🚨 Correct the coupled SLAYER determinant and Doppler-shift it by the kinetic-file E×B rotation #463), so γ can differ by about 0.1 % from run to run.
• Against the previous head (e84e1afa7), after merging develop 6fc89d021: in both cases every layer input, D_geo, D_R and Δ_crit are bit-identical. The γ values differ by at most 0.06 %, which is the same triangulation noise.
• Commits since 2762880: the Connor et al. 2015 entry in docs/development/references.md and class lines in the case files. An error on an unusable reference length was added and then taken back out, so that function is unchanged; the hazard is tracked in InnerLayer.SLAYER - BUGFIX - Fail clearly when the minor-radius derivative at a rational surface is unusable #505. None of this changes a computed value: regress --cases diiid_slayer_n1_toroidal --refs d008c30e1,e84e1afa7 --force on feynman gives 5 unchanged (D_R, Δ_crit, Q_root, γ, no-root flags).
• diiid_slayer_n1_toroidal: has no develop baseline, because develop errors on dc_type = "toroidal" without a dgeo_val. The χ∥-closure commit is therefore compared against the previous head, regress --refs d0e8a0261,276288065 (same job):

Regression Report: diiid_slayer_n1_toroidal
==========================================================================================
Ref 1: d0e8a0261  @ d0e8a0261 (2026-09-25)
       env: julia 1.11.9, x86_64-linux-gnu, manifest f826ee26 (pinned), 8 threads/8 BLAS
Ref 2: 276288065  @ 276288065 (2026-09-26)
       env: julia 1.11.9, x86_64-linux-gnu, manifest f826ee26 (pinned), 8 threads/8 BLAS
------------------------------------------------------------------------------------------
Quantity                            d0e8a0261  276288065  Diff               Status       
------------------------------------------------------------------------------------------
SLAYER resistive interchange D_R    [6 elem]   [6 elem]   0.0e+00            OK           
SLAYER critical-Δ offset D_c        [6 elem]   [6 elem]   1.368e+01 (4.95%)  ** CHANGED **
SLAYER Q_root [2/1,3/1,4/1]         [3 elem]   [3 elem]   2.959e-03 (0.42%)  ** CHANGED **
SLAYER γ_Hz [2/1,3/1,4/1]           [3 elem]   [3 elem]   1.105e+02 (3.49%)  ** CHANGED **
SLAYER no_root flags [2/1,3/1,4/1]  [3 elem]   [3 elem]   0.0e+00            OK           
Runtime (s)                         257.2s     255.1s                        --           
==========================================================================================
Summary: 3 changed, 2 unchanged

The per-surface numbers are in the closure section below. D_R and every layer input are identical.

Review package

A self-contained page with the derivation, the implementation data path (where the coordinate transformation enters), the symbolic and numerical verification ladder, and all figures embedded: Toroidal Critical-Δ Review Package.

Base branch

Originally based on and targeted at bugfix/slayer-dprime-reference-length (#403), which defines the r_s reference length K = r_s·dψ_N/dr used below. #403 has since merged into develop (f36ecef86), so this branch was rebased onto #403's final head 12f64da8c and the PR now targets develop with a single commit.

Derivation

Connor's Eq. 59, in Hamada coordinates with ψ(V) the toroidal flux, χ(V) the poloidal flux, primes d/dV, ι/2π = χ'/ψ' = 1/q, Λ = ψ'² ι'/2π, α = 2πn/χ':

$$\Delta'_{crit} = \frac{\pi^{3/2}}{2}\left(\frac{\chi_\parallel}{\chi_\perp}\right)^{1/4} V_s \left(\frac{\alpha^2\Lambda^2}{\langle B^2\rangle\langle|\nabla V|^2\rangle}\right)^{1/4}(-D_R)$$

In GPEC quantities (ψ_N normalized poloidal flux, chi1 = 2π·psio = dχ/dψ_N with χ the full poloidal flux, v1 = dV/dψ_N, q1 = dq/dψ_N):

$$\iota' = 2\pi\frac{d(1/q)}{dV} = -\frac{2\pi}{q^2},\frac{q_1}{v_1},\qquad \psi_t' = \frac{d\psi_t}{d\psi_N}\frac{d\psi_N}{dV} = \frac{q,\texttt{chi1}}{v_1},\qquad \Lambda = \left(\frac{q,\texttt{chi1}}{v_1}\right)^2\left(-\frac{q_1}{q^2 v_1}\right),\qquad \alpha = \frac{2\pi n}{\texttt{chi1}/v_1}$$

$$\langle B^2\rangle = \oint B^2,\frac{\texttt{jac}}{V'},d\theta,\qquad \langle|\nabla V|^2\rangle = v_1^2\oint|\nabla\psi_N|^2,\frac{\texttt{jac}}{V'},d\theta$$

Two points where the paper's notation has to be read carefully, both verified against the paper's own large-aspect-ratio limit (Eq. 61):

1. Connor writes ⟨A⟩ = (1/2π)∮A dθ, but with Hamada θ of period 1 this must be the normalized flux-surface average (⟨1⟩ = 1), i.e. the jac/V' weighting above. A literal 1/2π leaves a (2π)^{1/2} mismatch with Eq. 61.
2. Eq. 59 is the jump in dΨ/dY with Y = (V−V_s)/V_s. GPEC's layer stack, the rfitzp critical-Δ, and the Tearing - BUGFIX! - Convert Δ' to the r_s reference length before slab-layer matching #403-converted Δ' are all referenced to x̂ = (r−r_s)/r_s. Because Eq. 59 is derived at H = 0 (Frobenius exponent μ = 1/2), the conversion is linear: multiply by r_s(dV/dr)/V_s = K·v1/V_s. V_s cancels, so no volume integral is needed. Without this conversion Eq. 59 gives exactly half of rfitzp at large aspect ratio (V ∝ r²).

With the conversion the geometric factor reduces at large aspect ratio to √(n s r_s/R₀) with s = (r_s/q)dq/dr, which is exactly the factor inside rfitzp, so the two branches coincide there and the paper's Eq. 61 is recovered after dividing by r_s.

Implementation

• ResistGeometry (src/ForceFreeStates/Surfaces/ResistEval.jl) gains avg_dpsisq = ⟨|∇ψ_N|²⟩, an eighth column of the existing θ-loop.

• toroidal_dgeo (src/InnerLayer/SLAYER/LayerInputs.jl, exported) evaluates

  dgeo = k_ref · v1 · (α²Λ²/(⟨B²⟩ · v1² · ⟨|∇ψ_N|²⟩))^{1/4}
  

  with α and Λ as above. k_ref · v1 is the Y → x̂ reference conversion; this is the only place the coordinate transformation enters, and k_ref is the same K that delta_prime_to_rs_reference applies to Δ'.

• build_slayer_inputs derives dgeo_val whenever the surface carries a ResistGeometry (sing.restype, populated by resist_eval_all!), for every dc_type; only :toroidal consumes it, via _solve_dc_tmp's 0.5·(−D_R)·π^{3/2}·(χ∥/χ⊥)^{1/4}·dgeo, whose χ∥ closure now uses the same toroidal geometry (next section). It throws only if dc_type=:toroidal is requested on a surface without a ResistGeometry. An explicit dgeo_val still overrides.

Toroidal χ∥ closure

Eq. 59 needs χ∥/χ⊥. GPEC gets χ∥ by harmonically combining a collisional (Spitzer–Härm) value with a free-streaming value. That combination is solved self-consistently with the critical width W_d.

Until 276288065, :toroidal used Eq. 59's toroidal geometry only in the final prefactor. The W_d loop and the free-streaming χ∥ inside it still used the cylindrical field-line pitch n·s/R₀ (Fitzpatrick 1995, Phys. Plasmas 2 825, Eq. 132 and Sec. VII). Both now come from the balance that produces Eq. 59.

The balance. Connor's Eq. 18 with the Eq. 20 coefficients, written in x = V − V_s, reads

$$\chi_\perp\langle|\nabla V|^2\rangle,\partial_x^2T = \chi_\parallel,\frac{\alpha^2\Lambda^2}{\langle B^2\rangle},x^2T .$$

That is χ⊥k⊥² = χ∥k∥² with k∥(x) = |αΛ|·x/√⟨B²⟩ [1/m].

• No hidden factor in k∥. The choice α = 2πn/χ′ together with the phase exp(2πinu/χ′) fixes Hamada angles of period 1, so B·∇u = Λx exactly. The "Λx/φ" printed on Connor's p. 3 is dimensionally inconsistent. The paper's own Eq. 59 → 61 reduction confirms this.
• Transport width: x_T = (χ⊥/χ∥)^{1/4}/D_V, with D_V = (α²Λ²/(⟨B²⟩⟨|∇V|²⟩))^{1/4}.
• Eq. 59 in terms of it: (π^{3/2}/2)·D_R·V_s/x_T.

Critical width. Fitzpatrick's w_D (Connor Eq. 65) is 2√2 times the cylindrical x_T. In the x̂ reference, with the same conversion K·v1 = r_s·dV/dr as D_geo:

$$W_d = \frac{2\sqrt2,x_T}{K v_1} = \frac{\sqrt8,(\chi_\perp/\chi_\parallel)^{1/4}}{D_{geo}},\qquad \Delta'_{crit} = \frac{\sqrt2,\pi^{3/2}(-D_R)}{W_d},$$

The second relation is the same identity that links :lar and :rfitzp.

Free-streaming χ∥. The form and the 2/√π flux-limit factor are unchanged; only k∥ changes, to its value from the balance at x = W_d·K·v1:

$$\chi_{fs} = \frac{2v_{te}}{\sqrt\pi,K_\parallel W_d},\qquad K_\parallel = \frac{|\alpha\Lambda|,K,v_1}{\sqrt{\langle B^2\rangle}}.$$

K∥ is computed by the new toroidal_kpar.

In a circular cylinder, αΛ = −n s B₀/(4π²R₀²r²) and K·v1 = 4π²R₀r², so K∥ = n·s/R₀. That is exactly the previous code.

K∥ is not D_geo²/r_s. Algebraically, K∥ = (D_geo²/r_s)·g, with g = r_s√⟨|∇ψ_N|²⟩/K. The metric factor g is 1 in a cylinder, but 0.57–0.76 on the DIII-D-like surfaces. Taking K∥ from D_geo would put it 24–43 % off and make χ∥ depend on the radial label. So K∥ is computed separately.

Limits.

• Collisional (χ_SH ≪ χ_fs): K∥ drops out, and the result is unchanged from before.
• Free-streaming: the fixed point solves in closed form, χ∥^{3/4} = 2v_te·D_geo/(√(8π)·K∥·χ⊥^{1/4}), so Δ′_crit ∝ D_geo^{4/3}·K∥^{−1/3}. Relative to the previous closure, the factor is (D_geo·K_cyl/(D_cyl·K∥))^{1/3}.

That factor bounds the change on the DIII-D-like equilibrium (midplane label):

surface	ψ_N	D_geo/D_cyl	K∥/K_cyl	metric g	free-streaming Δ_crit change
2/1	0.518	1.172	1.044	0.761	+3.9 %
3/1	0.770	1.248	1.093	0.702	+4.5 %
4/1	0.893	1.316	1.134	0.655	+5.1 %
5/1	0.968	1.391	1.180	0.610	+5.6 %
6/1	0.993	1.457	1.210	0.570	+6.4 %

Measured with the diiid_slayer_n1_toroidal harness case, which uses the same equilibrium with psihigh = 0.9995 and so also reaches the 7/1. It compares d0e8a0261 (before) with 276288065 (after):

surface	Δ_crit before	Δ_crit after	change	free-streaming bound	γ before [s⁻¹]	γ after [s⁻¹]
2/1	16.89	17.52	+3.7 %	+3.9 %	−214.2	−229.8
3/1	16.22	16.91	+4.3 %	+4.5 %	−796.4	−820.3
4/1	44.93	47.10	+4.8 %	+5.1 %	−3164	−3274
5/1	276.1	289.8	+5.0 %	+5.6 %
6/1	49.80	52.27	+5.0 %	+6.4 %
7/1	26.56	28.13	+5.9 %

Every surface stays inside its free-streaming bound, and close to it: at these temperatures χ∥ is mostly free-streaming. D_R, D_geo and every layer input are bit-identical. The three tracked modes were already stable and become 3.0–7.3 % more stable.

Checks. All but the last are in the unit tests:

• Large-aspect-ratio limit (runtests_tj_analytic.jl): on TJ circular equilibria at ε = 0.1, 0.05 and 0.025, K∥/(n s/R₀) − 1 = (0.37–0.47)·ε². That is second order; the test requires |dev| ≤ ε² and a quartering from ε = 0.05 to 0.025. D_geo's first-order deviation is entirely the metric factor, g − 1 ≈ −2·(D_geo/D_cyl − 1), which checks the K∥–D_geo identity numerically.
• Closure algebra (runtests_slayer_params.jl, Test 1e):
  • Given cylindrical (D_geo, K∥), :toroidal reproduces :lar to 1e-12 and :rfitzp to 1e-8.
  • Both limits match their closed forms to 1e-8 for three (D_geo, K∥) pairs.
  • dgeo_val = 0 still gives no offset.
• Radial-label invariance (runtests_slayer_inputs.jl):
  • The toroidal Δ_crit/K is identical across :midplane, :flux and :volume to 1e-8. The previous cylindrical closure does not have this property, because its r-based shear and r_s change with the label.
  • K∥ and D_geo scale with K exactly.
  • A prescribed dgeo_val keeps the derived K∥, and without a ResistGeometry the closure falls back to the cylindrical K∥ with a warning.
• Scaling: K∥·a is invariant under B₀ → B₀/2 and (a, R₀) → 2(a, R₀).
• Independent re-derivation: a separate AI reviewer, given only the two papers and the diff, re-derived each step and confirmed it: k∥, W_d, the cylinder limit, the K∥–D_geo identity, the closed forms, label invariance, and the code. It also confirmed that v1_local, avg_bsq and avg_dpsisq are the normalized-average quantities the formulas need.

Still a modelling choice (also in the code comments and docstrings):

• Averaging: k∥ is taken in Connor's flux-surface-averaged sense, |αΛ|/√⟨B²⟩, so its variation along a field line is not resolved.
• Closure form: the closure itself is unchanged: Spitzer–Härm and free streaming, combined harmonically, with the 2/√π factor, and k∥ evaluated at the full W_d.
• Trapped electrons: there is no trapped-electron reduction of the free-streaming χ∥. Connor et al. treat χ∥ as an input, so that would be new modelling.
• Cylindrical models: :lar and :rfitzp keep the cylindrical geometry by definition.

Verification

• Symbolic (SymPy), 12/12: Λ ≡ ψ'χ'' − χ'ψ'' = ψ'²(ι/2π)'; the chain-rule images of ψ_t', Λ, α equal the code expressions; the code formula equals Eq. 59 × r_s(dV/dr)/V_s with V_s cancelling; the circular limit gives ½√(n s r/R) before and √(n s r/R) after conversion for an arbitrary q(r); Eq. 61·r_s equals the rfitzp formula through Connor's Eq. 65; the factor is invariant under B → λB and lengths → μ·lengths, while the Fortran one-power form scales as (λ/μ)^{-1/2}. The consequences that can be checked numerically are pinned by this PR's unit tests rather than by the script: the large-aspect-ratio limit and its linear ε-scaling (runtests_tj_analytic.jl), scale invariance, which the one-power form fails by √2 (runtests_tj_analytic.jl), and the exact :toroidal/:rfitzp prefactor identity (runtests_slayer_inputs.jl).
• Analytic chain Eq. 59 → Eq. 61 → Lutjens offset √2π^{3/2}(−D_R)/w_D (Eq. 65) → the rfitzp code formula: all consistent, pinning the prefactor.
• Large-aspect-ratio limit on TJ circular equilibria: dgeo/√(n s r_s/R₀) = 1.001–1.012 at ε = 0.05, ≤ 1.09 at ε = 0.3 near the edge The residual is a genuine O(ε) toroidal correction: it equals (0.144–0.185)·r_s/R₀ at every ε from 0.1 to 0.0125 and halves with ε to within 1–2 %. The unit test bounds it by the local r_s/R₀ at ε = 0.05 and requires it to halve from ε = 0.05 to 0.025 within 10 %, which a wrong O(1) factor (e.g. the Y-vs-r_s factor of 2, or a first-power ψ′) cannot satisfy.
• Scale invariance: dgeo unchanged to 6 digits under B₀ → B₀/2 and (a, R₀) → 2(a, R₀) (unit test). The Fortran STRIDE form on the slayer_growthrate branch of the Fortran GPEC, which carries one power of ψ_t' in Λ instead of two, fails this by √2 and is off by ψ_t'^{-1/2} (2.3× at ε = 0.1, 3.2× at ε = 0.05). That branch also feeds D_I rather than D_R and writes log(q) as shear; none of that is used here.
• Independent metric: ⟨B²⟩, ⟨|∇ψ_N|²⟩ and V' recomputed from R(ψ,θ), Z(ψ,θ), F(ψ), psio alone agree with ResistGeometry to 1e-7–1e-4 on interior surfaces of the circular and DIII-D-like cases.
• End-to-end prefactor: given cylindrical geometry, :toroidal reproduces :lar to 1e-12 and :rfitzp to 1e-8 (unit test). This checks the shared χ-matching prefactor. With toroidal geometry the closure checks are listed in the closure section above.
• Sign: an interchange-stable surface (D_R < 0) gets a positive, stabilizing offset (unit test).
• Wiring checks, not physics: the derivation reaches toroidal_dgeo with the surface's inputs, and the radial label reaches it only through k_ref.
• Tests at 276288065: runtests_slayer_params.jl 81/81, runtests_slayer_inputs.jl 84/84, runtests_tj_analytic.jl 28/28, runtests_slayer_runner.jl 79/79.

Radial-label sensitivity

The numbers below come from one-off runs, not from code in this PR, so this is how to regenerate them from the public API. On examples/DIIID-like_SLAYER_example, form the equilibrium and a fresh Riccati Δ′, then call build_slayer_inputs(equil, sings, profiles; rs_method=lab, dc_type=:toroidal) and the same with dc_type=:rfitzp, for lab in (:midplane, :flux, :volume). The ratio of the two dc_tmp values, and Δ′_rs/dc_tmp, give the rows. The TJ ε scan is runtests_tj_analytic.jl's _lar_equilibrium(ε) swept over ε.

DIII-D-like SLAYER example, n = 1 (Δ' diagonal reproduces #403's table: 9.097, −6.196, −16.374):

	2/1	3/1	4/1
Δ_crit toroidal / rfitzp, :midplane (current default)	1.17	1.25	1.32
Δ_crit toroidal / rfitzp, :flux	0.96	0.94	0.91
Δ_crit toroidal / rfitzp, :volume	1.00	1.01	1.01
label spread of the margin Δ'_rs/Δ_crit, rfitzp	31 %	32 %	41 %
label spread of the margin Δ'_rs/Δ_crit, toroidal	7.7 %	10 %	15 %

Because dgeo ∝ K exactly and the converted Δ' ∝ K^{2μ} with μ ≈ 1/2, the toroidal threshold margin is nearly label-invariant; the cylindrical rfitzp formula (∝ r_s/√(da/dψ)) is not. The residual toroidal spread in this table is ~1 % from K^{2μ−1} and ~6 % from χ∥ through the W_d loop. The table predates 276288065, which makes that loop toroidal and removes the second source exactly: Δ_crit/K is now identical across labels to 1e-8 (unit test). With either branch the 2/1 sits at about half its threshold under every label, whereas with dc_type = "none" it is unstable.

TJ circular ε scan (fresh Δ' per point, surfaces 2/1 and 3/1): the label spread of the margin and the toroidal/rfitzp deviation both vanish roughly linearly in ε (rfitzp spread 1.3–1.7 % at ε = 0.05, toroidal 0.2 %; toroidal/rfitzp 1.006–1.012 midplane, 0.995–0.998 flux). No residual at small ε, so no implementation artefact; finite-ε differences are label physics.

Notes for reviewers

1. Reference-conversion exponent — decision: linear. The critical-Δ is converted from Connor's Y = (V−V_s)/V_s reference to x̂ = (r−r_s)/r_s with the linear rule (Frobenius exponent ½, exact in the paper's H = 0 ordering). The alternative would be this PR's base's K^(2α) rule with the equilibrium's Mercier exponent α = √(−D_I). Linear is the defensible choice: Δ_crit is a layer-side quantity like the slab Δ̂(Q), W_d and rfitzp, all defined with the slab exponents ½ ± ½; only the outer Δ' genuinely carries ½ ± α, which is why Tearing - BUGFIX! - Convert Δ' to the r_s reference length before slab-layer matching #403 converts it with K^(2α). Applying c^(2α) to the offset would put an outer-region exponent on a quantity derived without one; the consistent generalization would be a finite-D_I layer theory with thermal conduction, which changes the structure of the D_R term and does not exist. Recorded in the toroidal_dgeo docstring. Cost of the alternative, c^(2α−1) on Δ_crit (always stabilizing), and the residual label-covariance mismatch of the threshold under the linear rule, K^(2α−1), on the DIII-D-like equilibrium (midplane label):

   surface	ψ_N	α	c	c^(2α−1)	K^(2α−1)
   3/2 (q = 1.5)	0.273	0.584	2.11	1.13	0.90
   2/1	0.518	0.540	2.21	1.07	0.99
   3/1	0.770	0.529	2.38	1.05	1.01
   4/1	0.893	0.558	2.60	1.12	1.04

   At the 2/1 and 3/1, where accurate γ matters most, the ambiguity is 5–7 % on Δ_crit and the label mismatch is ≤ 1 %. The 10 % mismatch at the 3/2 comes from the midplane label deep in the Shafranov-shifted core (K = 0.52), an argument for the flux label (InnerLayer.SLAYER - API - Choose one default radial label for the Fitzpatrick layer formalism #417), not for changing the exponent. The 5/1 (α = 0.83, D_R = −0.67) is outside the small-D_R ordering under either rule.

2. Default radial label. This PR inherits :midplane from Tearing - BUGFIX! - Convert Δ' to the r_s reference length before slab-layer matching #403. The evidence in InnerLayer.SLAYER - API - Choose one default radial label for the Fitzpatrick layer formalism #417 and the tables above point to :flux (Fitzpatrick's toroidal-flux label, the one the layer formulas are derived in) as the default for shaped plasmas; :volume only coincides with rfitzp by accident of near-circularity. Flipping the default moves rs, S, shear and every SLAYER result and is deliberately left to InnerLayer.SLAYER - API - Choose one default radial label for the Fitzpatrick layer formalism #417 as its own !-tagged PR.

3. D_geo is written for every dc_type once a ResistGeometry exists; it is a diagnostic and does not feed the dispersion unless dc_type = "toroidal".

---

🤖 Generated with Claude Code

[d-burg]

Review package

Toroidal Critical-Δ Review Package — a single self-contained page with the derivation (including the two notation readings that had to be checked against Eq. 61), the implementation data path showing the one place the coordinate transformation enters, the symbolic (SymPy, 12/12) and numerical verification ladder, the radial-label sensitivity tables and the TJ ε scan, with all figures embedded.

Decision recorded in the PR body and the toroidal_dgeo docstring: the Y → x̂ reference conversion of the critical-Δ stays linear (exponent ½, the paper's H = 0 ordering); the K^(2α) rule remains reserved for the outer Δ′, which genuinely carries the Mercier exponent. Cost of the alternative, c^(2α−1) on Δ_crit: 7 % at the 2/1, 5 % at the 3/1, 13 % at the 3/2, all in the stabilizing direction.

🤖 Generated with Claude Code

[Copilot]

Copilot review overview

🟡 Changes recommended

The new toroidal closure API has inconsistent K∥ defaults and lacks validation for invalid geometry values.

Review effort: Balanced
Findings: 2

Open (2)

• Derive SLAYERParameters default K∥ from preceding fields · New
• Validate kpar_val and geometric-factor overrides before iteration · New

What changed in this PR

Derives toroidal critical-Δ geometry from equilibrium data and propagates it through SLAYER’s thermal-conduction closure.

Changes:

• Adds D_geo, K∥, and ⟨|∇ψ_N|²⟩ calculations.
• Integrates derived geometry into SLAYER and HDF5 output.
• Adds analytic, unit, and regression coverage.

THIRD-PARTY HUMAN REVIEW IS REQUIRED BEFORE MERGE.

File	Description
src/​ForceFreeStates/​Surfaces/​ResistEval.jl	Adds the gradient flux-surface average.
src/​InnerLayer/​SLAYER/​LayerInputs.jl	Derives toroidal geometry per surface.
src/​InnerLayer/​SLAYER/​LayerParameters.jl	Applies geometry to the closure.
src/​InnerLayer/​SLAYER/​SLAYER.jl	Exports geometry helpers.
src/​InnerLayer/​InnerLayer.jl	Re-exports geometry helpers.
src/​Tearing/​Runner/​Control.jl	Documents control semantics.
src/​Tearing/​Runner/​HDF5Output.jl	Updates critical-Δ metadata.
test/​runtests_tj_analytic.jl	Tests asymptotic and scale limits.
test/​runtests_slayer_params.jl	Tests closure algebra.
test/​runtests_slayer_inputs.jl	Tests geometry derivation and wiring.
regression-harness/​cases/​diiid_slayer_n1.toml	Tracks derived D_geo.
regression-harness/​cases/​diiid_slayer_n1_toroidal.toml	Adds toroidal regression coverage.
docs/​development/​references.md	Adds the Connor et al. reference.

---

💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.