InnerLayer.SLAYER - BUGFIX - Restore the toroidal mode number in the diamagnetic frequencies

[d-burg]

Release note

• Audience: users
• Numerical impact: none on any shipped deck. Every SLAYER deck is nn = 1, where the restored factor is exactly 1.0 and the arithmetic is bit-identical. Results move for any run at n ≥ 2. (harness @ 47feada50)
• Migration: none.

build_slayer_inputs computed the electron and ion diamagnetic frequencies without the toroidal mode number, so Q_e and Q_i were low by a factor of n at every rational surface. Drift stabilization was therefore underestimated for any n ≥ 2 layer calculation.

The defect

_omega_star_at returns

ω_star_e = (2π / chi1) * (T_e * dn_e / n_e + dT_e)

which is (dp_e/dψ)/(e·n_e) — the n = 1 diamagnetic frequency. The physical quantity is ω_*e = k_θ·v_*e = (m/r_s)·(dp_e/dr)/(e·n_e·B). Writing that in flux coordinates, dψ/dr = r·B_φ/q absorbs the q in m = n·q and leaves an explicit

ω_*e = n · (dp_e/dψ) / (e · n_e)

Nothing downstream restored the n: Q_e = -tauk * omega_e in LayerParameters.jl, with no mode-number factor anywhere in the chain.

How it was found

Auditing whether the shear sign from #431 was load-bearing elsewhere in the layer physics. It is not — but the ExB rotation does pick up its mode number downstream (1im * mc.ntor * mc.rotation[k], Tearing/Dispersion/CoupledFullMatch.jl) while ω_* never did, and the omega kwarg to slayer_parameters is accepted but unused in the body. That asymmetry was the tell.

Confirmed against TJ

Fitzpatrick's TJ carries the mode number explicitly. Documentation/LayerParameters.tex:

$$\omega_{\ast k} = \frac{m_k,B_0,p_2'(\hat r_k)}{\mu_0,R_0^2,e,n_{ek},g_k,\hat r_k}$$

and TJ/Rational.cpp:180:

double wak = double(mres[k]) * B0 * ppk /mu0 /R0/R0 /e /nek /gk /rres[k];

Converting TJ's r-form to the flux form gives ratio TJ/GPEC = m/(q·g) = n.

The same cross-check validates the neighbouring quantities, which are not changed here: TJ's tHk = Lsk/m/VAk reduces to R₀/(n·s·V_A), matching GPEC's tau_h exactly; the W_d iteration and the :rfitzp critical-Δ match term for term.

Blast radius

quantity	effect
Q_e, Q_i	low by n — drift stabilization underestimated at n ≥ 2
iota_e = Q_e/(Q_e − Q_i)	unaffected — the factor cancels in the ratio
lu, tauk, D_norm, P_perp, delta_n, tau_r	unaffected

Entirely latent today: every SLAYER deck is nn = 1, and the only n > 1 example (Solovev_ideal_example_multi_n, nn_high = 2) is an ideal run that never builds a layer. The regression harness cannot catch this — there is no n > 1 SLAYER case.

Validation

• New test build_slayer_inputs: omega_star carries the toroidal mode number resolves the same surface at n = 1 and n = 2 and asserts Q/tauk doubles while iota_e does not move

• Verified non-vacuous: with the fix reverted, exactly the two scaling assertions fail and the iota_e assertion still passes — the predicted signature

• runtests_slayer_inputs 38/38

• Full suite julia -t 4 --project=. test/runtests.jl at 53fffc0d0: 62 testsets, 2348 assertions, 0 failures

• Post-review refactor (_omega_star_at now takes the mode number) re-verified: runtests_slayer_inputs 38/38, runtests_slayer_params 68/68, runtests_slayer_runner 73/73

Regression report

diiid_slayer_n1, develop (349a0c262) vs this head (47feada50), both sides --force, identical pinned environments.

Regression Report: diiid_slayer_n1
==========================================================================================
Ref 1: 349a0c262  @ 349a0c262 (2026-09-04)
       env: julia 1.11.6, arm64-apple-darwin24.0.0, manifest 7e5c34ad (pinned), 8 threads/8 BLAS
Ref 2: 47feada50  @ 47feada50 (2026-09-04)
       env: julia 1.11.6, arm64-apple-darwin24.0.0, manifest 7e5c34ad (pinned), 8 threads/8 BLAS
------------------------------------------------------------------------------------------
Quantity                            349a0c262  47feada50  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 Q_root [2/1,3/1,4/1]         [3 elem]   [3 elem]   5.1e-06            OK
SLAYER ω_Hz [2/1,3/1,4/1]           [3 elem]   [3 elem]   1.8e-06            OK
SLAYER γ_Hz [2/1,3/1,4/1]           [3 elem]   [3 elem]   1.276e-01 (0.02%)  ** 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)                         828.8s     453.5s                        --
==========================================================================================
Summary: 1 changed, 16 unchanged

Inert on the shipped deck, as predicted. Every layer quantity is bit-identical at 0.0e+00 — including Q_e and Q_i, now tracked by this case, which is the direct statement that the restored factor does nothing at n = 1.

The lone γ line is the threaded root search's own reproducibility, not this change. Its measured same-source floor is 0.076/0.122/0.145 Hz per surface (#418), and #431 saw 1.276e-01 on this same case; here it is 1.276e-01 Hz (0.02%). With every layer input bit-identical, a 0.13 Hz shift in the extracted root can only come from the root search. An earlier run of this same pair at 25d285eb8 returned 1.9e-05 on this line — the spread between those two runs is itself the floor being described.

The n = 2 demonstration

There is no n > 1 SLAYER case in the harness, so this defect was invisible to it. While preparing this PR I built one — diiid_slayer_n2, the n = 1 deck re-run via [overrides] at nn = 2 — to produce the measurement below. It is not part of this PR: on review (@jhalpern30) it was dropped as regression bloat, since it cost 273s/ref to guard something a unit test covers for free. Its two novel quantities, Q_e and Q_i, were folded into diiid_slayer_n1 instead, and the n-scaling is guarded by test/runtests_slayer_inputs.jl.

The measurement stands as validation even though the case is gone. develop (e09795119) vs 21fd7240e:

SLAYER minor radius rs              [13 elem]  [13 elem]  0.0e+00              OK
SLAYER Lundquist S                  [13 elem]  [13 elem]  0.0e+00              OK
SLAYER D_norm                       [13 elem]  [13 elem]  0.0e+00              OK
SLAYER tauk                         [13 elem]  [13 elem]  0.0e+00              OK
SLAYER iota_e                       [13 elem]  [13 elem]  0.0e+00              OK
SLAYER Q_e                          [13 elem]  [13 elem]  1.158e+00 (100.00%)  ** CHANGED **
SLAYER Q_i                          [13 elem]  [13 elem]  2.122e+00 (100.00%)  ** CHANGED **
SLAYER Q_root [3/2,4/2,5/2]         [3 elem]   [3 elem]   8.407e-01 (100.02%)  ** CHANGED **
SLAYER ω_Hz [3/2,4/2,5/2]           [3 elem]   [3 elem]   1.094e+04 (100.02%)  ** CHANGED **
SLAYER γ_Hz [3/2,4/2,5/2]           [3 elem]   [3 elem]   8.441e-01 (0.13%)    ** CHANGED **
SLAYER no_root flags [3/2,4/2,5/2]  [3 elem]   [3 elem]   0.0e+00              OK

With every geometric and resistive input bit-identical, Q_e and Q_i move by exactly 100% — the restored factor n = 2 — and the mode rotation ω_Hz doubles with them (drift rotation ∝ ω_*), while γ_Hz moves 0.13%. The root inventory is unchanged. On develop, an n = 2 run underestimates the diamagnetic drive by exactly half; any 3/2 analysis needs this PR first.

The harness fix that made [overrides] usable here (8ae06cd23, plus the 4c70d4232 gitignore) is kept — it is a general capability for any deck with relative file references, not specific to the removed case.

Notes for reviewers

The mode number is applied inside _omega_star_at(ψ, n_tor), which returns the physical ω_. (It was originally applied at the call site; @jhalpern30 correctly pushed back that a helper named for ω_ should return ω_*.) The parameter is n_tor rather than n because n_e and n_i are densities in the same scope.

The compute_omega_star=false branch reads profiles.omega_e/omega_i, which are literal zeros in production: run_slayer.jl:44 hardcodes zeros(npsi), compute_omega_star defaults to true, and no caller in src/ passes it. The kinetic file schema has no ω_* column either, so nothing can supply one. That branch is a library-only path.

Note

Found alongside a second, independent discrepancy in D_norm — see the companion PR. The two are in different files and do not conflict.

[jhalpern30]

I am fully trusting you on the physics here - all of my comments are more general suggestions. Leaving this as approved for simplicity whenever you're done looking at the comments since this definitely doesn't need a second runthrough, but if you disagree/decide to ignore any of the below please comment why for posterity

[jhalpern30]

Personally, I think the something = an entire if-else block is kinda confusing code, but it seems like that's in a lot of this file so its probably ok.

But if I am understanding this PR correctly - the actual omega* has an included n, so its a little bit deceptive to have _omega_star_at return something that isn't actually omega* and needs to be multiplied by n. I think its more clear to build it into the helper even though Claude argues against it _omega_star_at(psi, n). Then you can remove this 3 line comment up here and just leave it as is

[jhalpern30]

When I see override, I think of something that is already set; in this case its just another branch of how to set this value.

I think more true to form would be

ω_e_use, ω_i_use = (prof.omega_e, prof.omega_i)
# Override ω_*e, ω_*i with spline-derivative values when requested
compute_omega_star && (ω_e_use, ω_i_use = _omega_star_at(psi, n))

or something like that. The PR comment does suggest that prof.omega_e/i are never actually used, but I didn't factor that in here since I am unfamiliar with the rest of the relevant code

[jhalpern30]

So the entire purpose of this regression case is to guard the n-dependence of the layer module... that seems like overkill for a regression test that takes 5 minutes to run. I speak selfishly because my (and I'm assuming other's) regress running is automatically chosen by Claude. So if I change anything related to this module's code, Claude will queue up both of these cases to test something relatively small.

So personally, I am against regression case bloat and don't want them to take the place of unit tests. So my question here for your consideration - if the n1 and n2 case here actually needed and worth the additional regression harness runtime, which can already get quite long? Or is this something where we could keep a single regression case (perhaps n2 to try to catch n-dependent things like this in the future) and supplement with unit tests that check for specifics to slayer more cheaply?