[ebursch]

Release note

• Audience: users
• Numerical impact: SLAYER layer parameters and growth rates move on every deck. P_⊥ is now built by default from the perpendicular particle diffusivity D_⊥, which on diiid_slayer_n1 is 0.22–0.33 of the old electron-χ value, so γ moves by +31 % (2/1), +18 % (3/1) and +24 % (4/1) while ω moves by ≤ 0.4 %. The new analytic χ profiles of the DIII-D-like kinetic file alone move γ by ≤ 1.5 %. The layer-overlap cap of diiid_truncation_n1 is unchanged (its widths move 2 %). (harness @ 4ce6bfe)
• Migration: in [SLAYER], rename chi_perp to chi_perp_e (also the chi_perp keyword of build_slayer_inputs and slayer_parameters). The new default P_perp_model = "D_perp" needs the ion diffusivity χ⊥,i: add a chi_i profile to the kinetic file or set chi_perp_i, or set P_perp_model = "chi_perp_e" to keep the previous P_⊥. tau_E is the whole-plasma confinement time.

SLAYER's perpendicular Prandtl number P_⊥, which enters the inner-layer Δ(Q) with P_φ, can now be built several ways through [SLAYER] P_perp_model, and by default comes from the perpendicular particle diffusivity of Fitzpatrick (2022), Eq. 16. The critical Δ keeps the electron χ⊥,e in every mode, so the two are no longer forced to be the same diffusivity. The kinetic HDF5 schema gains an optional chi_i profile.

P_perp_model	P_⊥	Reference
D_perp (default)	τ_R D_⊥/r_s², with D_⊥ = c_β² η/μ₀ + (2/3)(1 − c_β²)(n T_e′ χ⊥,e + n T_i′ χ⊥,i)/p′	Fitzpatrick, Phys. Plasmas 29, 032507 (2022), Eqs. 4–6, 16
chi_perp_e	τ_R χ⊥,e/r_s² (the only option before; Fitzpatrick's TJ convention)	—
chi_perp_i	τ_R χ⊥,i/r_s²	—
P_phi	P_⊥ = P_φ	Fitzpatrick, Phys. Plasmas 30, 092512 (2023), Eq. 143
c_beta	P_⊥ = c_β², the χ⊥ → 0 limit of D_perp	J.-K. Park, Phys. Plasmas 29, 072506 (2022), Eq. 9 with K = 0
tau_E	P_⊥ = P_φ = τ_R/τ_E(r_s), with τ_E(r) = r²/χ_E and a uniform χ_E = (3/16)a²/τ̄_E	Fitzpatrick (2023), Eqs. 37–38, 81–83, 131

What changed

• P_⊥ models in slayer_parameters, selected by [SLAYER] P_perp_model and passed through TearingProblem and the layer-overlap scan, so the scan and the tearing stage see one plasma.
• D_perp default. A run that needs χ⊥,i and has neither a chi_i profile nor chi_perp_i stops with a message naming both. slayer_parameters keeps chi_perp_e as its own default, since it has no profile gradients.
• tau_E takes the whole-plasma τ̄_E, converted to a uniform χ_E with κ_e = 3/16 (the T_e ∝ 1 − r̂² value of Eq. 38) and the boundary minor radius from the flux label √(2ψ_t/B₀).
• Kinetic file. TkMkr_D3Dlike_Hmode_kinetic.h5 carries analytic χ profiles (below) in place of profiles derived from a measured discharge, and gains chi_i; n, T and ω_E are byte-identical to develop's. The forms and their sources are in the file's provenance attribute and in the stability docs.
• Truncation deck pins P_perp_model = "chi_perp_e": under D_perp the DIII-D-like layers are about 30 % narrower and do not overlap inside ψ = 1, so the pin keeps the layer-overlap cap exercised.
• Harness case diiid_slayer_n1_chi_perp_e keeps a baseline for the previous model.
• Docs: the P_⊥ models and the analytic χ in stability.md, chi_i in the kinetic schema, the Fitzpatrick 2022 and 2023 papers in docs/resources, and the Park citation in the SLAYER file headers corrected to Phys. Plasmas 29, 072506 (2022).
• An old deck's chi_perp is rejected with a hint that it is now chi_perp_e.

chi_e   = 2.75 + ψ + 3 exp(−25 (ψ − 0.85)²)        [m²/s, ψ = ψ_N]
chi_i   = 1/2 + (3/2) ψ²⁰ + (5/4) exp(−25 (ψ − 0.7)²)
chi_phi = 3/2 − (5/4) (1 − ψ)⁴

These are not fits to a discharge. Their magnitudes are of the same order as the TRANSP power-balance χ_e and χ_i of DIII-D H-modes in Hong et al., Phys. Plasmas 30, 072512 (2023), Fig. 5; χ_e keeps the factor-2 raise that makes the default deck exercise the overlap cut; χ_φ tracks χ_i (Prandtl number of order one, rising outward) as reported for DIII-D by Yoshida et al., Nucl. Fusion 52, 123005 (2012). They do not resolve the pedestal transport barrier, where power-balance χ falls to about 0.1–0.3 m²/s (Callen et al., Nucl. Fusion 50, 064004 (2010)).

Regression report

regress --refs 9fb32c9a6,4ce6bfea5 --force on feynman, all 20 cases (develop against this branch with develop merged). Tests and the docs build pass at the head.

case	result
diiid_slayer_n1	P_⊥ ×0.22–0.33 at all six surfaces; γ 222.6 → 291.0 (2/1), −288.5 → −339.5 (3/1), −1124 → −1395 s⁻¹ (4/1); ω ≤ 0.4 %; no root flag changes
diiid_slayer_n1_chi_perp_e (new)	no develop run (develop rejects the new key); at this head its γ is within 0.8 %, 0.9 % and 1.5 % of develop's diiid_slayer_n1, the effect of the analytic χ alone
diiid_truncation_n1	cap, psilim 0.99498, qlim 6.2 and surfaces unchanged; Eq. 100 widths 2.05 %
every other case	unchanged

Full regression report

Regression Report: diiid_slayer_n1
===========================================================================================
Ref 1: 9fb32c9a6  @ 9fb32c9a6 (2026-10-20)
       env: julia 1.11.9, x86_64-linux-gnu, manifest 8e24ca66 (pinned), 8 threads/8 BLAS, jit build
Ref 2: 4ce6bfea5  @ 4ce6bfea5 (2026-10-20)
       env: julia 1.11.9, x86_64-linux-gnu, manifest 8e24ca66 (pinned), 8 threads/8 BLAS, jit build
-------------------------------------------------------------------------------------------
Quantity                            9fb32c9a6  4ce6bfea5  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]   5.589e+01 (77.66%)  ** CHANGED **
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]   6.787e-03 (1.09%)   ** CHANGED **
SLAYER ω_Hz [2/1,3/1,4/1]           [3 elem]   [3 elem]   4.435e+01 (0.40%)   ** CHANGED **
SLAYER γ_Hz [2/1,3/1,4/1]           [3 elem]   [3 elem]   2.713e+02 (24.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)                         298.2s     262.5s                         --
===========================================================================================
Summary: 4 changed, 13 unchanged

Regression Report: diiid_truncation_n1
=====================================================================================================
Ref 1: 9fb32c9a6  @ 9fb32c9a6 (2026-10-20)
       env: julia 1.11.9, x86_64-linux-gnu, manifest 8e24ca66 (pinned), 8 threads/8 BLAS, jit build
Ref 2: 4ce6bfea5  @ 4ce6bfea5 (2026-10-20)
       env: julia 1.11.9, x86_64-linux-gnu, manifest 8e24ca66 (pinned), 8 threads/8 BLAS, jit build
-----------------------------------------------------------------------------------------------------
Quantity                               9fb32c9a6      4ce6bfea5      Diff               Status
-----------------------------------------------------------------------------------------------------
layer-overlap cap applied              1              1              0.0e+00            OK
layer-overlap psilim                   9.985081e-01   9.985081e-01   0.0e+00            OK
first overlapping surface index        7              7              0.0e+00            OK
diffusive-resistive layer widths (Δψ)  [7 elem]       [7 elem]       1.320e-04 (2.05%)  ** CHANGED **
qlim                                   6.200000e+00   6.200000e+00   0.0e+00            OK
psilim                                 9.949767e-01   9.949767e-01   0.0e+00            OK
mpert                                  37             37             0.0e+00            OK
# singular surfaces                    5              5              0.0e+00            OK
singular q values                      [5 elem]       [5 elem]       0.0e+00            OK
singular psi locations                 [5 elem]       [5 elem]       0.0e+00            OK
plasma energy Re(ep[1])                -1.372382e+00  -1.372382e+00  0.0e+00            OK
total energy Re(et[1])                 8.005423e-01   8.005423e-01   0.0e+00            OK
Runtime (s)                            227.1s         199.9s                            --
=====================================================================================================
Summary: 1 changed, 11 unchanged

diiid_slayer_n1_chi_perp_e fails on develop, which does not know P_perp_model, so the harness exits with one failed run; gal_resistive_pe extracts nothing on either ref, as on develop. Every other case is unchanged.

Notes for reviewers

• Taken over from @ebursch while he is on leave. His commits are kept; the develop merges, the review changes, the D_perp default, the analytic χ and the docs are on top.
• Why D_⊥ for both the layer and Eq. 100. Fitzpatrick (2025) derives the diffusive-resistive width (Eq. 100) from the layer equations of Fitzpatrick (2022), where P_⊥ = τ_R/τ_⊥ with τ_⊥ = r²/D_⊥ multiplies the diffusion of the density variable (2022, Eq. 29). The 2025 paper writes χ_E in that slot and sets it to 1 m²/s, relying on the 2022 finding that D_⊥ is comparable to the energy diffusivities. Eq. 16 makes that precise; on these profiles it gives D_⊥ = 0.22–0.32 χ_e at the outer rational surfaces, because the weights n T_s′/p′ sum to less than one where density gradients carry part of p′.
• Every P_⊥ model is μ₀χ/η, so r_s cancels except through η's trapped fraction. Only tau_E carries a radial label, through a; it uses the flux label whatever rs_method is, so the pending radial-label default change (InnerLayer.SLAYER - API - Choose one default radial label for the Fitzpatrick layer formalism #417) does not touch these models.
• κ_e is fixed at 3/16 rather than computed from the profile, because Eq. 38 needs dT_e/dr at the edge, which a pedestal makes ill-conditioned (0.045–0.12 on this profile depending on how the edge slope is taken). The docstring states the order-2–3 shape uncertainty, and that in H-mode the χ_E implied by τ̄_E falls below the core power-balance χ.
• n and T in the example file are unchanged: they are the long-standing DIII-D-like H-mode profiles, rough fits that every pinned harness quantity depends on.
• Overlap with InnerLayer.SLAYER - DOCS - 🌱 Describe chi_tor as the toroidal momentum diffusivity, not a heat diffusivity #471, which edits the chi_tor docstrings in the same files: whichever merges second needs a small conflict resolution.

---

NO MERGE WITHOUT HUMAN REVIEW. This PR needs a third-party human reviewer's approval of the head that lands.

🤖 Generated with Claude Code

[ebursch]

Effect of each P_perp_model on the DIII-D-like SLAYER example

Setup: examples/DIIID-like_SLAYER_example at 27d0744, run once per model with only P_perp_model changed. The run is uncoupled, uses AMR and is validity-gated. χ⊥,e, χ⊥,i and χ_φ come from the chi_e, chi_i and chi_phi profiles in TkMkr_D3Dlike_Hmode_kinetic.h5. For tau_E I used τ_E = 0.1 s, a typical DIII-D H-mode value, which is not derived from this equilibrium. All three surfaces found a root in every run.

P_perp_model	P_⊥ (2/1, 3/1, 4/1)	P_φ (2/1, 3/1, 4/1)	γ [Hz] (2/1, 3/1, 4/1)	Δγ vs default (2/1, 3/1, 4/1)	ω [Hz] 2/1
chi_perp_e	50.3, 44, 35.7	34.7, 24.7, 19	192, -251, -1012	—	-10527
chi_perp_i	29.9, 26.4, 12.5	34.7, 24.7, 19	213, -274, -1217	+11 %, -9 %, -20 %	-10523
P_phi	34.7, 24.7, 19	34.7, 24.7, 19	206, -277, -1123	+8 %, -10 %, -11 %	-10525
D_perp	16.4, 17.5, 10.9	34.7, 24.7, 19	242, -296, -1249	+26 %, -18 %, -23 %	-10516
c_beta	0.0133, 0.0088, 0.00645	34.7, 24.7, 19	313, -588, -1669	+63 %, -134 %, -65 %	-10133
tau_E	34.5, 30.8, 27.6	34.5, 30.8, 27.6	206, -264, -1034	+8 %, -5 %, -2 %	-10525

The Δγ column is (γ − γ_default)/|γ_default|, so + means more unstable. For scale, the run-to-run AMR noise on γ is about 0.01 %.

Takeaways

• Only c_beta changes the answer qualitatively. Its P_⊥ = c_β² ≈ 0.01, three to four orders of magnitude below every other model. The 2/1 drive rises about 60 %, the stable 3/1 and 4/1 roots damp 65–135 % harder, and ω on the 2/1 shifts about 4 %. Every other model shifts ω by less than 0.1 %.
• The χ-based models move γ by roughly 10–25 %.
  • chi_perp_i uses the example's χ⊥,i, which is about 0.4–0.6 of χ⊥,e at these surfaces. That halves P_⊥ and moves γ by about 10 % at the 2/1 and 3/1 and 20 % at the 4/1.
  • D_perp gives the smallest P_⊥ of the transport-based models, so it moves γ the furthest, by 18–26 %.
  • Across these models, a lower P_⊥ makes the 2/1 more unstable and the 3/1 and 4/1 more stable.
• P_phi and tau_E land close together. With τ_E = 0.1 s, τ_R/τ_E happens to be close to the χ_φ-based P_φ on the 2/1. tau_E also overrides P_φ, so its 3/1 and 4/1 sit nearer the default. γ scales with the chosen τ_E, so these numbers only illustrate the model.
• chi_perp_e (the default) reproduces develop exactly; see the regression report in the PR description.

The unit tests (Pkg.test() at 27d0744 on feynman) pass.

---

🤖 Generated with Claude Code

[ebursch]

Example kinetic-file χ profiles vs the DIII-D 147131 AOT source

Each panel compares TkMkr_D3Dlike_Hmode_kinetic.h5 (red dashed) with the AOT profile it was built from (blue): DIII-D 147131 at 2300 ms, mapped ρ_N → ψ_N with EFIT01 RHOVN. The grey lines mark the 2/1, 3/1 and 4/1 surfaces used in the SLAYER example. The bottom-right panel overlays the three .h5 profiles on one axis.

Dataset	Status	Source
chi_e (χ⊥,e)	existing, unchanged	AOT CHI_E, smoothed
chi_phi (χ_φ)	existing, unchanged	AOT ANG_MOM_DIFF, smoothed
chi_i (χ⊥,i)	added in this PR	AOT CHI_I; spurious zero at ψ_N = 1 dropped; Gaussian-smoothed with σ = 12 of 201 ρ_N points

At the SLAYER surfaces (new chi_i):

Surface	ψ_N	chi_i (.h5)	AOT CHI_I	.h5 vs AOT	chi_i / chi_e
2/1	0.518	1.13	1.06	+7 %	0.59
3/1	0.770	1.80	2.10	−15 %	0.60
4/1	0.893	1.15	0.82	+39 %	0.35

Notes

• Smoothing width: σ = 12 is the width that best rebuilds the existing chi_e from AOT CHI_E, about 10 % RMS in log over the profile.
• Edge treatment differs: the existing chi_e and chi_phi also smooth over the pedestal well, the edge spike in CHI_E and the negative ANG_MOM_DIFF beyond ψ_N ≈ 0.94, so they stay positive and smooth to the edge. The file doesn't record how those edge treatments were done. The new chi_i uses only the Gaussian filter, so it follows the AOT well near ψ_N ≈ 0.93 and the rise to about 2 m²/s at the edge more closely.
• 4/1 surface: the 4/1 sits on the inner flank of that well, where the smoothing lifts chi_i 39 % above AOT. A narrower σ would track AOT more closely there, at the cost of keeping more of the AOT structure elsewhere.

---

🤖 Generated with Claude Code

[ebursch]

@logan-nc @d-burg This adds a variety of options for perpendicular diffusion inputs. It doesn't change the default yet, so selecting the chi_perp_e is the same as what currently exists. You can see in the comparison that the only really different one is the c_beta which is good news (and encouragement to have this increased fidelity vs the Fortran slayer implementation). I think that we should move to full D_perp as default since there is not a literature justification for chi_perp_e ~ D_perp and if people already need chi_phi and chi_perp_e inputs than its not too crazy to also have them input chi_perp_i. It's encouraging that the tau_E choice is right order of magnitude (using tau_E = chi_phi and P_perp = P_phi) so that is also a helpful, user friendly, and control friendly option. After this gets merged, then either @d-burg or I should make a PR changing the default to D_perp and removing the chi_perp_e and chi_perp_i options. Let me know if you have any questions on this. Thanks!

[logan-nc]

At a high level, this looks good. I agree that if someone has chi_phi and chi_perp_e, they likely have chi_perp_i on hand. If they don't have things, or are lazy then letting them assert a "reasonable" tau_E for their device is a great way to get results quick.

[d-burg]

Thanks -- the switch is cleanly done, the default path is the same arithmetic as before, and the rename is complete. I checked the new models against the papers: D_perp matches Fitzpatrick (2022) Eqs. 4–6 and 16, c_beta matches Park (2022) and is the exact χ⊥ → 0 limit of the coded D_perp, and P_phi matches Fitzpatrick (2023) Eq. 143. tau_E does not yet match Eq. 131 (inline at line 400). I also compared the old and new example file: only chi_i and the provenance text were added, and every pre-existing dataset is byte-identical.

A mix of top level feedback from myself and claude:

Blocking

1. Docs CI. Two docstring lines parse as Markdown links (suggestion inline at LayerParameters.jl:220).
2. Rebase onto develop. This now conflicts in 7 files since #298 landed. Note for the binary file: #298 raised chi_e in TkMkr_D3Dlike_Hmode_kinetic.h5, and this branch's copy predates that, so chi_i needs to be added to develop's file rather than taking this branch's version, otherwise the chi_e change is undone. The chi_i/chi_e ratios in your table will change accordingly.

On making D_perp the default and removing chi_perp_e / chi_perp_i

I agree D_perp is the better-founded model, but I'd settle a few things before switching:

• D_perp is Fitzpatrick's derived perpendicular particle diffusivity: c_β²η⊥/μ₀ plus (2/3)(1 − c_β²)(n T_e′χ⊥,e + n T_i′χ⊥,i)/p′, built from the electron and ion energy diffusivities. By comparison chi_perp_e is an electron-only proxy, so I'd describe it in the docstring and release note as the backward-compatible choice rather than "the" P⊥. One caveat for the default: where the temperature profiles are flat, Eq. 16 reduces to the c_beta floor, which your table shows moves γ by 60–130 %, so flat-temperature regions of a real profile would change answers a lot.
• tau_E needs fixing before it is advertised as the quick option: Eq. 131's τ_E is the local r²/χ_E, not the whole-plasma confinement time, so as documented this mode is off by r̂_s²/κ_e (details inline).
• χ⊥,e stays a required input in every mode because it sets the critical-Δ width, so it can't go away even if the model option does. This PR is already a breaking rename of a key every SLAYER deck sets; I'd rather agree the end state now so users migrate once.
• Please add one harness case with a non-default P_perp_model (your table already has the reference numbers), so a later default change has a baseline.

Data question. The example file's provenance attribute now names DIII-D 147131 at 2300 ms as the source of chi_i; the previous text named no shot. Has releasing this AOT-derived profile in a public repository been cleared under the DIII-D data policy? Since the provenance is being rewritten anyway, it should also use the dataset names chi_e and chi_phi, not "chi_perp" and "chi_tor".

Small, outside the diff

• docs/src/kinetic_forces.md documents the kinetic HDF5 schema without chi_i and calls chi_e "the" perpendicular heat diffusivity; the run_slayer docstring likewise mentions only chi_e/chi_phi.
• An old deck with chi_perp = ... is rejected by the generic unknown-key error, which is the important part, but a one-line hint when chi_perp is among the unknown keys ("renamed to chi_perp_e") would save users a search through several dozen names.
• Fitzpatrick (2022) and (2023) are cited for three of the six models but aren't in docs/resources; adding them would let later reviewers check the equations.

[d-burg]

These two lines are what fails the Documentation job: [m²/s] (used by ...) and [s] (required by ...) parse as Markdown links (the parser skips the space between ] and (), and this docstring is spliced into inner_layer.md. A comma after the unit bracket, as the build_slayer_inputs docstring already does, removes both. The suggestion also adds :D_perp to the models that use chi_perp_i.

Suggested change

	- `chi_perp_i` -- ion perpendicular heat diffusivity [m²/s] (used by `P_perp_model=:chi_perp_i`)
	- `dlnn`, `dlnTe`, `dlnTi` -- radial log-gradients of n, T_e, T_i in any common radial variable
	(only their ratios enter; required by `P_perp_model=:D_perp`)
	- `tau_E` -- whole-plasma energy confinement time [s] (required by `P_perp_model=:tau_E`)
	- `chi_perp_i` -- ion perpendicular heat diffusivity [m²/s], used by `P_perp_model=:chi_perp_i` and `:D_perp`
	- `dlnn`, `dlnTe`, `dlnTi` -- radial log-gradients of n, T_e, T_i in any common radial variable
	(only their ratios enter; required by `P_perp_model=:D_perp`)
	- `tau_E` -- whole-plasma energy confinement time [s], required by `P_perp_model=:tau_E`

[d-burg]

I checked this against Fitzpatrick (2022) Eqs. 4–6 and 16: the expression and the docstring formula match the paper. Two small things.

The intermediate ratios have a singularity the formula itself does not: with a vanishing ion (or electron) pressure gradient, frac_i and tau_F both go infinite and the product gives Inf·0 = NaN, so this throws although Eq. 16 is finite there. Writing the weights over one common denominator is the same arithmetic (η_e τ_e/(1+η_e) = T_e·dlnTe/dp and likewise for ions), leaves only p′ = 0 singular, and makes the physical content visible: the bracket is (n T_e′ χ⊥,e + n T_i′ χ⊥,i)/p′.

Suggested change

	# D_⊥ from Fitzpatrick, Phys. Plasmas 29, 032507 (2022), Eqs. 4-6 and 16 (η_⊥ = η_∥), with
	# η_s/(1+η_s) = dlnT_s/(dlnn+dlnT_s) and τ = (T_e/T_i)(1+η_e)/(1+η_i) written flat-density safe.
	frac_e = dlnTe / (dlnn + dlnTe)
	frac_i = dlnTi / (dlnn + dlnTi)
	tau_F = (t_e / t_i) * (dlnn + dlnTe) / (dlnn + dlnTi)
	D_perp = c_beta^2 * eta / MU_0 +
	(2 / 3) * (1 - c_beta^2) * (frac_e * tau_F / (1 + tau_F) * chi_perp_e + frac_i / (1 + tau_F) * chi_perp_i)
	# D_⊥ from Fitzpatrick, Phys. Plasmas 29, 032507 (2022), Eqs. 4-6 and 16 (η_⊥ = η_∥); the weights η_s τ_s/(1+η_s) equal n T_s'/p'.
	dp = t_e * (dlnn + dlnTe) + t_i * (dlnn + dlnTi)
	D_perp = c_beta^2 * eta / MU_0 +
	(2 / 3) * (1 - c_beta^2) * (t_e * dlnTe * chi_perp_e + t_i * dlnTi * chi_perp_i) / dp

It would also help to say in the docstring what Eq. 16 implies at the edges, since users may be surprised: where T_e and T_i are flat, the χ terms vanish and P⊥ falls to the c_beta value (≈ 0.01 here), and where the species pressure gradients have opposite signs (hollow density) the weights can exceed 1. A test at a flat-temperature point would pin the first.

[d-burg]

tau_E and the log-gradients default to NaN so that a model which needs them fails loudly, but chi_perp_i defaults to 1.0 m²/s, so :chi_perp_i or :D_perp called without it returns a plausible-looking number. I'd default it to NaN here and in build_slayer_inputs, and move the isfinite(...) && ... > 0 check out of the :D_perp branch to just after the if block so it applies to P_perp for every model: today a zero, negative or NaN χ⊥,i (a partly-zero profile, or spline undershoot) gives P⊥ ≤ 0 in :chi_perp_i mode with no error.

[d-burg]

Small citation fix: in Park (2022) the C² term is in Eq. (9) and is defined in the text just below Eq. (10); Eq. (10) itself carries P, which is our P_tor. I'd cite "Eq. 9, with C² = c_β² + (1 − c_β²)K and K = 0" here and at line 423. On line 230, please also define τ_e = τ/(1+τ), τ_i = 1/(1+τ), η_s = dlnT_s/dlnn, and note that this τ is (T_e/T_i)(1+η_e)/(1+η_i), the opposite orientation to the tau = T_i/T_e field of SLAYERParameters.

[d-burg]

In tau_E mode P_tor is replaced by τ_R/τ_E, so the file's chi_phi and the scalar chi_tor are ignored; that is only stated in the slayer_parameters docstring. Please say it in the SLAYERControl docstring and the TOML annotation too (suggestion on the example deck), and log one @info when the mode is active. This warning should also depend on the model: as written it reports a chi_tor fallback in tau_E mode, and a chi_perp_e fallback in c_beta, P_phi and tau_E mode with dc_type = "none", where those values do not enter the result.

[d-burg]

Suggested change

	P_perp_model = "chi_perp_e" # Perpendicular Prandtl number model: "chi_perp_e", "chi_perp_i", "P_phi", "D_perp", "c_beta", or "tau_E" (needs tau_E [s])
	P_perp_model = "chi_perp_e" # Perpendicular Prandtl number model: "chi_perp_e", "chi_perp_i", "P_phi", "D_perp", "c_beta", or "tau_E" (needs tau_E [s]; also replaces the χ_φ-based P_tor)

[d-burg]

With T_e = T_i this check passes whichever way round the temperature ratio is written, and the struct's own tau is T_i/T_e, so I'd add a T_e ≠ T_i point (hand numbers: τ = 3/2, τ_e = 3/5, τ_i = 2/5; not run). One build_slayer_inputs test with P_perp_model=:D_perp would also cover the gradient forwarding from the profile splines, which is otherwise untested. I haven't run this test myself so this is just a suggestion, you should confirm on your end whether it checks out.

Suggested change

	@test d1.P_perp ≈ d1.tau_r * D_hand / d1.rs^2 rtol = 1e-12
	@test d1.P_perp ≈ d1.tau_r * D_hand / d1.rs^2 rtol = 1e-12
	# T_e = 2T_i ⇒ τ = 3/2, τ_e = 3/5, τ_i = 2/5: fixes the orientation of the temperature ratio
	d2 = slayer_parameters(; base..., t_e=2000.0, P_perp_model=:D_perp, grads..., chi_perp_i=2.0)
	D_hand2 = d2.c_beta^2 * d2.eta / MU_0 + (2 / 3) * (1 - d2.c_beta^2) * ((2 / 3) * (3 / 5) * 1.0 + (3 / 4) * (2 / 5) * 2.0)
	@test d2.P_perp ≈ d2.tau_r * D_hand2 / d2.rs^2 rtol = 1e-12

[d-burg]

Missing, non-positive and NaN tau_E are caught here, but inf gets through and only fails later inside slayer_parameters, after the equilibrium and stability stages. Adding isfinite makes the two checks agree.

Suggested change

	ctrl.P_perp_model !== :tau_E || (ctrl.tau_E !== nothing && ctrl.tau_E > 0) ||
	ctrl.P_perp_model !== :tau_E || (ctrl.tau_E !== nothing && isfinite(ctrl.tau_E) && ctrl.tau_E > 0) ||

[d-burg]

I checked this against Fitzpatrick (2023), Phys. Plasmas 30, 092512. Eq. 131 is P⊥ = P_φ = τ_R/τ_E, but the τ_E there is the local energy confinement time τ_E(r) = r²/χ_E with χ_E uniform (Eqs. 81–83), not the whole-plasma confinement time. Eqs. 140–142 make this explicit: P⊥(r̂) = (τ_R0/τ_E0)·ĵ_oh(r̂) with τ_E0 = a²/χ_E, so the r² cancels and P⊥ = μ₀χ_E/η(r_s). The whole-plasma time from W_e/P_e or a scaling law is τ̄_E = κ_e a²/χ_E (Eqs. 37–38), with κ_e = (3/2)∫₀¹ T̂_e r̂ dr̂ / [−dT̂_e(1)/dr̂].

With tau_E documented as the whole-plasma confinement time, this code's P⊥ differs from the paper's by r̂_s²/κ_e. For a T_e profile like (1 − r̂²), κ_e ≈ 0.19, so the 3/1 and 4/1 surfaces come out about 2–4× larger than Eq. 131 gives, and the error grows outward.

Two ways to make it match: take the whole-plasma τ̄_E and use P⊥ = τ_R(r_s)·κ_e a²/(r_s² τ̄_E), which needs κ_e from the T_e profile and the minor radius; or take a uniform χ_E (= κ_e a²/τ̄_E) as the input, which makes this mode the same as the χ-based P⊥ with one number for every surface. Either way, the docstring should say which confinement time is meant and cite Eqs. 81–83 alongside Eq. 131.

[logan-nc]

Thanks for the thorough review @d-burg
I think @ebursch will tackle these changes post paternity leave.

I'd just like to reinforce that DIII-D data cannot be included in this open repo directly. The approach for the "DIIID-like" example should be to cite a paper if available and make a roughly matching functional form of the profile needed. Ideally, note the functional form somewhere in annotations / docs so it's clear its not a direct profile fit to data (eg "y = y0 (1 -x^2) where y0=3.0 roughly matches [Name Journal 2026]" or some such).

Addressed in 4ce6bfe: develop merged, the review changes folded in, D_perp made the default and the example chi profiles replaced by analytic forms. See the updated description.