diff --git a/docs/development/references.md b/docs/development/references.md
index 43f3daf74..978e6826f 100644
--- a/docs/development/references.md
+++ b/docs/development/references.md
@@ -76,6 +76,16 @@ The SLAYER slab inner-layer solver and the resistive-layer width diagnostics are
   - Location: `docs/resources/2026-Burgess-Tearing Stability Prediction Combining Toroidal Calculations With a Two-Fluid Slab Layer.pdf`
   - Describes: The two-fluid slab layer, its four regimes, and the `S^(1/3)` outer-inner matching
 
+- **Fitzpatrick (2022)**: "Influence of anomalous perpendicular transport on linear tearing mode dynamics in tokamak plasmas"
+  - Location: `docs/resources/2022-Fitzpatrick-Influence of anomalous perpendicular transport on linear tearing mode dynamics in tokamak plasmas.pdf`
+  - Published: Physics of Plasmas **29**, 032507 (2022)
+  - Describes: The perpendicular particle diffusivity D_⊥ built from the electron and ion energy diffusivities (Eqs. 4-6, 16), used by `P_perp_model = :D_perp`
+
+- **Fitzpatrick (2023)**: "Error-field penetration thresholds in ohmically heated ITER and SPARC plasmas"
+  - Location: `docs/resources/2023-Fitzpatrick-Error-field penetration thresholds in ohmically heated ITER and SPARC plasmas.pdf`
+  - Published: Physics of Plasmas **30**, 092512 (2023)
+  - Describes: P_⊥ = P_φ = τ_R/τ_E with the local τ_E = r²/χ_E (Eqs. 81-83, 131, 143) and the whole-plasma τ̄_E = κ_e a²/χ_E (Eqs. 37-38), used by `P_perp_model = :P_phi` and `:tau_E`
+
 - **Fitzpatrick (2025)**: "Response of a magnetically diverted tokamak plasma to a resonant magnetic perturbation"
   - Location: `docs/resources/2025-Fitzpatrick-Response of a magnetically diverted tokamak plasma to a resonant magnetic perturbation.pdf`
   - Published: Nuclear Fusion (2025), doi 10.1088/1741-4326/ae4fdd (open access: arXiv:2511.07666)
diff --git a/docs/resources/2022-Fitzpatrick-Influence of anomalous perpendicular transport on linear tearing mode dynamics in tokamak plasmas.pdf b/docs/resources/2022-Fitzpatrick-Influence of anomalous perpendicular transport on linear tearing mode dynamics in tokamak plasmas.pdf
new file mode 100644
index 000000000..184539270
Binary files /dev/null and b/docs/resources/2022-Fitzpatrick-Influence of anomalous perpendicular transport on linear tearing mode dynamics in tokamak plasmas.pdf differ
diff --git a/docs/resources/2023-Fitzpatrick-Error-field penetration thresholds in ohmically heated ITER and SPARC plasmas.pdf b/docs/resources/2023-Fitzpatrick-Error-field penetration thresholds in ohmically heated ITER and SPARC plasmas.pdf
new file mode 100644
index 000000000..014a7212c
Binary files /dev/null and b/docs/resources/2023-Fitzpatrick-Error-field penetration thresholds in ohmically heated ITER and SPARC plasmas.pdf differ
diff --git a/docs/src/kinetic_forces.md b/docs/src/kinetic_forces.md
index ced2a590d..e8e616fb8 100644
--- a/docs/src/kinetic_forces.md
+++ b/docs/src/kinetic_forces.md
@@ -12,8 +12,9 @@ extension:
   - **HDF5 (`.h5`/`.hdf5`)** — the GPEC kinetic schema (recommended). Datasets
     at the file root: required `psi` (normalized poloidal flux), `n_e`, `T_i`,
     `T_e`, `omega_E` (and `n_i`, defaulting to `n_e` if omitted); optional
-    `omega_tor`, `chi_e` (perpendicular heat diffusivity ``\chi_\perp``),
-    `chi_phi` (toroidal momentum diffusivity ``\chi_\phi``), and `zeff`
+    `omega_tor`, `chi_e` and `chi_i` (electron and ion perpendicular heat diffusivities
+    ``\chi_{\perp,e}``, ``\chi_{\perp,i}``), `chi_phi` (toroidal momentum diffusivity
+    ``\chi_\phi``), and `zeff`
     (effective charge ``Z_\mathrm{eff}``, dimensionless). Each dataset
     carries a `units` attribute; the root carries `schema_version` and
     `provenance`. Densities are m⁻³, temperatures eV, frequencies rad/s,
@@ -25,7 +26,7 @@ extension:
     `psi_n n_i n_e T_i[eV] T_e[eV] omega_E`, retained for backward
     compatibility. Header rows are skipped.
 
-The NTV calculation consumes `n_i, n_e, T_i, T_e, omega_E`; `chi_e`/`chi_phi`/`zeff`,
+The NTV calculation consumes `n_i, n_e, T_i, T_e, omega_E`; `chi_e`/`chi_i`/`chi_phi`/`zeff`,
 when present, are carried for the resistive-layer (SLAYER) analysis and ignored
 here. That analysis takes ``Z_\mathrm{eff}(\psi)`` from the `zeff` dataset when the file
 has one, otherwise from `n_i/n_e` with the `[SLAYER]` `impurity_charge`, otherwise
diff --git a/docs/src/stability.md b/docs/src/stability.md
index 01f84bdba..8d5e19eaa 100644
--- a/docs/src/stability.md
+++ b/docs/src/stability.md
@@ -333,28 +333,49 @@ resolve it from, so the cap cannot be applied there and setting the flag warns i
 `psilim_from_layer_overlap = true`. It then caps `qlim` before the `rational_offset_fraction` step, so `rational_offset_fraction`
 still selects the final surface from inside the cap. A cap beyond `psihigh` has no effect.
 
-The flag is off by default, so every shipped deck is unchanged by it. With the shipped
-DIII-D-like profiles the overlap point is ψ ≈ 0.9985: outside the ideal and Riccati decks'
-`psihigh = 0.995`, where the cap could not bind even if enabled, but inside the SLAYER deck's
-0.9995, where enabling it would exclude the q = 7 surface. `examples/DIIID-like_truncation_example`
-raises `psihigh` to 0.9999 with the flag on: the cap excludes q = 7 and q = 8, and the resulting
-domain (`psilim ≈ 0.99498`, last surface q = 6) reproduces the hand-chosen `psihigh = 0.995` of the
-ideal deck. Whenever `psilim` lands past the overlap point, with the flag on or off, GPEC warns;
-on the shipped decks that includes the SLAYER deck.
-
-**Which perpendicular diffusivity.** The two formulas that take `chi_perp` want different
-quantities. The layer width goes through `P_perp = τ_R/τ_E` with `τ_E = r²/χ_E`, and χ_E is the
-one-fluid perpendicular *energy* diffusivity, the gradient-weighted combination of the electron and
-ion channels, `χ_E = (n_e χ_e ∇T_e + n_i χ_i ∇T_i)/(n_e ∇T_e + n_i ∇T_i)`. The critical-Δ island
-width `W_d ∝ (χ_⊥/χ_∥)^{1/4}` describes electron-temperature flattening against an electron parallel
-conductivity, so its χ_⊥ is the electron channel. GPEC currently feeds both from one profile, the
-kinetic file's `chi_e`, which in the shipped DIII-D-like data is the electron power-balance
-diffusivity: the right channel for the critical Δ, and a proxy for χ_E in the layer width. This is
-the same single-χ convention as Fitzpatrick's TJ code, which also assumes `T_e = T_i`.
-Separating the two is deferred: it needs an ion `chi_i` in the kinetic-profile schema, and the
-power-balance ion diffusivity available upstream is not yet usable (it comes back negative across
-most of the radius on the cases checked). Since both widths scale as `χ^{1/4}`, a factor-2 error in
-χ moves them by about 19%, but that can be enough to move the overlap cut by a surface.
+The flag is off by default, so every shipped deck is unchanged by it. With the default
+`P_perp_model = "D_perp"` the layers on the shipped DIII-D-like profiles do not overlap inside
+ψ = 1. With the electron-channel `"chi_perp_e"` they are about 40% wider and the overlap point is
+ψ ≈ 0.9985: outside the ideal and Riccati decks' `psihigh = 0.995`, but inside the SLAYER deck's
+0.9995, where the cap would exclude the q = 7 surface. `examples/DIIID-like_truncation_example`
+pins `"chi_perp_e"` and raises `psihigh` to 0.9999 with the flag on: the cap excludes q = 7 and
+q = 8, and the resulting domain (`psilim ≈ 0.99498`, last surface q = 6) reproduces the hand-chosen
+`psihigh = 0.995` of the ideal deck. Whenever `psilim` lands past the overlap point, with the flag
+on or off, GPEC warns.
+
+**Which perpendicular diffusivity.** The two formulas that take a perpendicular diffusivity want
+different quantities. The layer width goes through `P_perp`, which multiplies the diffusion of the
+total pressure; under the closure `T_s = T_s(n)` it is set by the perpendicular particle diffusivity
+`D_⊥` of Fitzpatrick, Phys. Plasmas **29**, 032507 (2022), Eq. 16, which weights the electron and ion energy diffusivities by
+`n T_s′/p′`. The critical-Δ island width `W_d ∝ (χ_⊥/χ_∥)^{1/4}` describes electron-temperature
+flattening against an electron parallel conductivity, so its χ_⊥ is the electron channel `chi_e`
+in every case. `[SLAYER] P_perp_model` selects how `P_perp` is built: from `D_⊥` (`"D_perp"`,
+the default), from `chi_e` alone (the single-χ convention of Fitzpatrick's TJ code), from
+`chi_i`, from the toroidal Prandtl number, from its `c_β²` floor, or from a whole-plasma
+confinement time (`"tau_E"`); see `slayer_parameters`. `"tau_E"` is for a plasma with no χ
+profiles: its uniform χ_E sets both `P_perp` and the toroidal `P_tor`, so the file's `chi_phi`
+and `[SLAYER] chi_tor` are not used in that mode. On the DIII-D-like
+profiles `D_⊥` is 0.22–0.32 of `chi_e` at the outer rational surfaces, because the weights sum to
+less than one where density gradients carry part of `p′`. Since the widths scale as `χ^{1/4}`, a
+factor-2 change in χ moves them by about 19%, which can be enough to move the overlap cut by a
+surface.
+
+**The shipped diffusivities.** The DIII-D-like kinetic file carries analytic χ profiles in the
+normalized poloidal flux ψ, in m²/s:
+
+```
+chi_e   = 2.75 + ψ + 3 exp(−25 (ψ − 0.85)²)
+chi_i   = 1/2 + (3/2) ψ²⁰ + (5/4) exp(−25 (ψ − 0.7)²)
+chi_phi = 3/2 − (5/4) (1 − ψ)⁴
+```
+
+They are not fits to a measured discharge. Their magnitudes are of the same order as the TRANSP
+power-balance χ_e and χ_i reported for DIII-D H-modes by Hong *et al.*, Phys. Plasmas **30**,
+072512 (2023), Fig. 5; `chi_e` includes a factor-2 raise so that the default deck exercises the
+resistive-layer overlap cut. `chi_phi` is set to track `chi_i` (Prandtl number of order one, rising
+outward), as reported for DIII-D by Yoshida *et al.*, Nucl. Fusion **52**, 123005 (2012). The forms
+do not resolve the H-mode pedestal transport barrier, where power-balance χ falls to about
+0.1–0.3 m²/s (Callen *et al.*, Nucl. Fusion **50**, 064004 (2010)).
 
 | Dataset | Contents |
 |---|---|
diff --git a/examples/DIIID-like_SLAYER_example/gpec.toml b/examples/DIIID-like_SLAYER_example/gpec.toml
index c32ca67c4..d5f74e8da 100644
--- a/examples/DIIID-like_SLAYER_example/gpec.toml
+++ b/examples/DIIID-like_SLAYER_example/gpec.toml
@@ -82,10 +82,11 @@ scan_mode               = "amr"      # Q-plane scan strategy: "amr" (adaptive re
 coupling_mode           = "uncoupled" # "uncoupled" (per-surface) or "coupled" (multi-surface determinant)
 D_c_type                = "none"     # Critical-Δ offset selector: "none", "lar", "fitzpatrick", or "toroidal"
 ion_mass                = 2.0        # Ion mass in proton-mass units (2.0 = deuterium)
-zeff                    = 1.0        # fallback only; Z_eff(ψ) comes from the kinetic file (its zeff profile, else n_i/n_e with zimp)
+zeff                    = 1.0        # fallback only; Z_eff(ψ) comes from the kinetic file (its zeff profile, else n_i/n_e with impurity_charge)
 impurity_charge         = 6.0        # Impurity charge closing quasi-neutrality when Z_eff(ψ) is derived from n_i/n_e (6.0 = carbon)
-chi_perp                = 1.0       # fallback only; the kinetic file supplies χ⊥(ψ)
+chi_perp_e              = 1.0       # fallback only; the kinetic file supplies χ⊥,e(ψ)
 chi_tor                 = 1.0       # fallback only; the kinetic file supplies χ_φ(ψ)
+P_perp_model            = "D_perp"   # Perpendicular Prandtl number model: "D_perp" (default), "chi_perp_e", "chi_perp_i", "P_phi", "c_beta", or "tau_E" (needs tau_E [s]; also replaces the χ_φ-based P_tor)
 # omega_E_kHz           = {"2/1" = 0.0, "3/1" = 3.0} # override of the E×B rotation Ω_E/2π per unit n [kHz], keyed by m/n; unlisted surfaces use the kinetic file
 pole_threshold_adaptive = true       # SLAYER |Δ| spans many orders; adapt to 10·median(|Δ|)
 filter_above_poles      = true       # Discard roots above the highest pole γ
@@ -93,7 +94,7 @@ filter_outside_re       = true       # Condition the above-pole filter on the +
 polish_roots            = true       # refine each root to the true zero; enables the validity gate
 store_scan              = false      # per-surface scans omitted from HDF5 (uncoupled has one per surface)
 # Standardized kinetic-profile file (GPEC HDF5 kinetic schema), shared with the
-# sibling ideal example; carries n_e, T_e, T_i, omega_E, chi_e (χ⊥), chi_phi (χ_φ).
+# sibling ideal example; carries n_e, T_e, T_i, omega_E, chi_e (χ⊥,e), chi_i (χ⊥,i), chi_phi (χ_φ).
 profile_file = "../DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode_kinetic.h5" # Kinetic-profile file (relative to the run dir), required when SLAYER is enabled
 
 [SLAYER.scan_grid]
diff --git a/examples/DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode_kinetic.h5 b/examples/DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode_kinetic.h5
index dab58dc3f..185329e82 100644
Binary files a/examples/DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode_kinetic.h5 and b/examples/DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode_kinetic.h5 differ
diff --git a/examples/DIIID-like_ideal_example/gpec.toml b/examples/DIIID-like_ideal_example/gpec.toml
index d706bbed7..2f1345eca 100644
--- a/examples/DIIID-like_ideal_example/gpec.toml
+++ b/examples/DIIID-like_ideal_example/gpec.toml
@@ -88,7 +88,7 @@ write_outputs_to_HDF5 = true            # Write ErrorFields/CoilSensitivities/ t
 verbose = true                          # Log per-coil-set progress and linearity diagnostics
 
 [KineticForces]
-kinetic_file = "TkMkr_D3Dlike_Hmode_kinetic.h5"  # GPEC HDF5 kinetic schema (psi,n_i,n_e,T_i,T_e,omega_E,+chi_e,chi_phi). Legacy .gpeckf/.kin still readable.
+kinetic_file = "TkMkr_D3Dlike_Hmode_kinetic.h5"  # GPEC HDF5 kinetic schema (psi,n_i,n_e,T_i,T_e,omega_E,+chi_e,chi_i,chi_phi). Legacy .gpeckf/.kin still readable.
 nn = 1                                  # Toroidal mode number (matches nn_low/nn_high)
 bounce_harmonic_max = 6                 # Bounce harmonics retained
 ion_charge = 1                          # Main-ion charge (deuterium)
diff --git a/examples/DIIID-like_truncation_example/gpec.toml b/examples/DIIID-like_truncation_example/gpec.toml
index 341c01263..22001feb1 100644
--- a/examples/DIIID-like_truncation_example/gpec.toml
+++ b/examples/DIIID-like_truncation_example/gpec.toml
@@ -44,6 +44,8 @@ psilim_from_layer_overlap = true   # Cap psilim where adjacent resistive layers
 
 [SLAYER]
 # The growth-rate analysis is not run here; the section is present because the layer-overlap scan
-# takes its kinetic profiles, and its chi_perp(psi), from the SLAYER profile_file.
+# takes its kinetic profiles, and its chi_perp_e(psi), from the SLAYER profile_file. P_perp is pinned to the
+# electron channel so the cap binds: with the default D_perp the layers are ~30% narrower and do not overlap inside psi = 1.
 enabled = false                # Run the SLAYER tearing-mode analysis
+P_perp_model = "chi_perp_e"    # Perpendicular Prandtl number model (pinned; see above)
 profile_file = "../DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode_kinetic.h5" # Kinetic-profile file (relative to the run dir), required when SLAYER is enabled
diff --git a/regression-harness/cases/diiid_slayer_n1_chi_perp_e.toml b/regression-harness/cases/diiid_slayer_n1_chi_perp_e.toml
new file mode 100644
index 000000000..2968a1750
--- /dev/null
+++ b/regression-harness/cases/diiid_slayer_n1_chi_perp_e.toml
@@ -0,0 +1,10 @@
+# Regression case: the diiid_slayer_n1 deck with P_perp built from the electron diffusivity chi_e alone,
+# the model GPEC used before the D_perp default. It keeps a baseline for the non-default P_perp models
+# and for the electron-channel convention of Fitzpatrick's TJ code.
+[case]
+name = "diiid_slayer_n1_chi_perp_e"
+base = "diiid_slayer_n1"   # Inherit example_dir, [quantities.*] and [overrides] from this case
+description = "DIII-D-like H-mode equilibrium, n=1, SLAYER tearing-mode analysis with P_perp_model = chi_perp_e"
+
+[overrides]
+"SLAYER.P_perp_model" = "chi_perp_e"   # P_perp from chi_e alone, as before the D_perp default
diff --git a/src/Equilibrium/KineticProfiles.jl b/src/Equilibrium/KineticProfiles.jl
index 3d4132d72..a73bb5d35 100644
--- a/src/Equilibrium/KineticProfiles.jl
+++ b/src/Equilibrium/KineticProfiles.jl
@@ -66,7 +66,8 @@ temperatures eV, frequencies rad/s, diffusivities m²/s.
 | `T_i`/`T_e`         | ion / electron temperature                                                                                        |
 | `omega_E`           | ExB rotation                                                                                                      |
 | `omega_tor`         | toroidal rotation (optional)                                                                                      |
-| `chi_e`             | perpendicular heat diffusivity χ⊥ (shipped data: electron power-balance χ_e; see the stability docs)              |
+| `chi_e`             | electron perpendicular heat diffusivity χ⊥,e                                                                      |
+| `chi_i`             | ion perpendicular heat diffusivity χ⊥,i                                                                           |
 | `chi_phi`           | toroidal momentum diffusivity χ_φ                                                                                 |
 | `zeff`              | effective charge Z_eff(ψ) supplied directly (optional); the resistive-layer analysis prefers it to its closure    |
 | `species_densities` | named per-species densities (e.g. `"n_D"`, `"n_T"`) for explicit multi-ion input; `nothing` for ASCII / no extras |
@@ -81,6 +82,7 @@ struct KineticProfileData
     omega_E::Union{Nothing,Vector{Float64}}
     omega_tor::Union{Nothing,Vector{Float64}}
     chi_e::Union{Nothing,Vector{Float64}}
+    chi_i::Union{Nothing,Vector{Float64}}
     chi_phi::Union{Nothing,Vector{Float64}}
     zeff::Union{Nothing,Vector{Float64}}
     species_densities::Union{Nothing,Dict{String,Vector{Float64}}}
@@ -88,11 +90,11 @@ struct KineticProfileData
 end
 
 function KineticProfileData(; psi, n_i=nothing, n_e=nothing, T_i=nothing, T_e=nothing,
-    omega_E=nothing, omega_tor=nothing, chi_e=nothing, chi_phi=nothing, zeff=nothing,
+    omega_E=nothing, omega_tor=nothing, chi_e=nothing, chi_i=nothing, chi_phi=nothing, zeff=nothing,
     species_densities=nothing, provenance="")
     _v(x) = x === nothing ? nothing : Float64.(collect(x))
     return KineticProfileData(Float64.(collect(psi)), _v(n_i), _v(n_e), _v(T_i), _v(T_e),
-        _v(omega_E), _v(omega_tor), _v(chi_e), _v(chi_phi), _v(zeff), species_densities, String(provenance))
+        _v(omega_E), _v(omega_tor), _v(chi_e), _v(chi_i), _v(chi_phi), _v(zeff), species_densities, String(provenance))
 end
 
 const _KINETIC_H5_EXTS = (".h5", ".hdf5", ".he5")
@@ -100,7 +102,7 @@ const _KINETIC_H5_EXTS = (".h5", ".hdf5", ".he5")
 # Dataset units written into / expected from the HDF5 kinetic schema.
 const _KINETIC_H5_UNITS = Dict("psi" => "normalized poloidal flux", "n_i" => "m^-3",
     "n_e" => "m^-3", "T_i" => "eV", "T_e" => "eV", "omega_E" => "rad/s",
-    "omega_tor" => "rad/s", "chi_e" => "m^2/s", "chi_phi" => "m^2/s", "zeff" => "dimensionless")
+    "omega_tor" => "rad/s", "chi_e" => "m^2/s", "chi_i" => "m^2/s", "chi_phi" => "m^2/s", "zeff" => "dimensionless")
 
 """
     read_kinetic_file(path; group="/") -> KineticProfileData
@@ -139,7 +141,7 @@ function _read_kinetic_h5(path::AbstractString; group::AbstractString="/")
         haskey(ats, "provenance") && (prov = string(read(ats["provenance"])))
         # Non-standard datasets named `n_*` (e.g. "n_D", "n_T") are named per-species density
         # profiles for explicit multi-ion input; other names are reserved for future schema fields.
-        standard = ("psi", "n_i", "n_e", "T_i", "T_e", "omega_E", "omega_tor", "chi_e", "chi_phi", "zeff")
+        standard = ("psi", "n_i", "n_e", "T_i", "T_e", "omega_E", "omega_tor", "chi_e", "chi_i", "chi_phi", "zeff")
         extras = Dict{String,Vector{Float64}}()
         for k in keys(g)
             (k in standard || !startswith(k, "n_")) && continue
@@ -149,7 +151,7 @@ function _read_kinetic_h5(path::AbstractString; group::AbstractString="/")
         return KineticProfileData(; psi=Float64.(vec(read(g["psi"]))),
             n_i=rd("n_i"), n_e=rd("n_e"), T_i=rd("T_i"), T_e=rd("T_e"),
             omega_E=rd("omega_E"), omega_tor=rd("omega_tor"),
-            chi_e=rd("chi_e"), chi_phi=rd("chi_phi"), zeff=rd("zeff"),
+            chi_e=rd("chi_e"), chi_i=rd("chi_i"), chi_phi=rd("chi_phi"), zeff=rd("zeff"),
             species_densities=(isempty(extras) ? nothing : extras), provenance=prov)
     end
 end
@@ -180,6 +182,7 @@ function write_kinetic_h5(path::AbstractString, data::KineticProfileData;
         put("omega_E", data.omega_E)
         put("omega_tor", data.omega_tor)
         put("chi_e", data.chi_e)
+        put("chi_i", data.chi_i)
         put("chi_phi", data.chi_phi)
         put("zeff", data.zeff)
         # Named per-species densities (multi-ion input); dataset names must be `n_*` to round-trip.
diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl
index 886d551c8..1da5ed652 100755
--- a/src/GeneralizedPerturbedEquilibrium.jl
+++ b/src/GeneralizedPerturbedEquilibrium.jl
@@ -564,20 +564,21 @@ function _overlap_slayer_inputs(inputs, dir_path::AbstractString)
     (inputs isa AbstractDict && haskey(inputs, "SLAYER")) || return nothing
     ctrl = Runner.slayer_control_from_toml(inputs["SLAYER"])
     isempty(ctrl.profile_file) && !ctrl.enabled && return nothing
-    loaded = try
-        Runner._load_profiles(ctrl, dir_path)
-    catch err
-        ctrl.enabled && rethrow()
-        @warn "Layer-overlap scan skipped: could not read the [SLAYER] profile_file." exception = err
-        return nothing
-    end
     # chi(psi) and Z_eff(psi) from the file where it carries them, else the control's scalar
     # fallbacks -- the same precedence run_slayer uses. Every other plasma parameter the layer
     # widths depend on is the run's own [SLAYER] setting, so the scan and the SLAYER stage see
     # one plasma.
+    loaded, chi = try
+        l = Runner._load_profiles(ctrl, dir_path)
+        l, Runner._slayer_chi_inputs(l, ctrl; warn=false)
+    catch err
+        ctrl.enabled && rethrow()
+        @warn "Layer-overlap scan skipped: could not read the [SLAYER] profile_file or the diffusivities its P_perp_model needs." exception = err
+        return nothing
+    end
     return (; profiles=loaded.profiles,
-        chi_perp=loaded.chi_perp === nothing ? ctrl.chi_perp : loaded.chi_perp,
-        chi_tor=loaded.chi_tor === nothing ? ctrl.chi_tor : loaded.chi_tor,
+        chi.chi_perp_e, chi.chi_perp_i, chi.chi_tor,
+        tau_E=ctrl.tau_E, P_perp_model=ctrl.P_perp_model,
         zeff=loaded.zeff === nothing ? ctrl.zeff : loaded.zeff,
         bt=ctrl.bt, mu_i=ctrl.ion_mass, theta=ctrl.theta_sample,
         resistivity_model=Runner._build_resistivity_model(ctrl.resistivity_model),
@@ -604,7 +605,8 @@ function _layer_overlap_scan(overlap_inputs, intr::ForceFreeStatesInternal, equi
     # Eq. (30), so the scan is driven in that label regardless of the SLAYER default.
     scans = [
         Tearing.resistive_layer_overlap(equil, overlap_inputs.profiles; n_tor=n, rs_method=:flux,
-            chi_perp=overlap_inputs.chi_perp, chi_tor=overlap_inputs.chi_tor,
+            chi_perp_e=overlap_inputs.chi_perp_e, chi_perp_i=overlap_inputs.chi_perp_i, chi_tor=overlap_inputs.chi_tor,
+            tau_E=overlap_inputs.tau_E, P_perp_model=overlap_inputs.P_perp_model,
             bt=overlap_inputs.bt, mu_i=overlap_inputs.mu_i, zeff=overlap_inputs.zeff, theta=overlap_inputs.theta,
             resistivity_model=overlap_inputs.resistivity_model, lnLambda_form=overlap_inputs.lnLambda_form)
         for n in intr.nlow:intr.nhigh
diff --git a/src/InnerLayer/SLAYER/LayerInputs.jl b/src/InnerLayer/SLAYER/LayerInputs.jl
index 80ebf533a..ed9bc0aa1 100644
--- a/src/InnerLayer/SLAYER/LayerInputs.jl
+++ b/src/InnerLayer/SLAYER/LayerInputs.jl
@@ -186,10 +186,15 @@ profiles, without an intermediate file round-trip.
 
   - `zeff`      -- effective charge. Scalar or a callable of `psi`, evaluated per
     surface (default `1.0`).
-
-  - `chi_perp`  -- perpendicular heat diffusivity [m²/s]. Scalar or a
+  - `chi_perp_e`  -- electron perpendicular heat diffusivity [m²/s]. Scalar or a
     callable of `psi` (default `1.0`).
-
+  - `chi_perp_i` -- ion perpendicular heat diffusivity [m²/s], scalar or callable of `psi`;
+    required by `P_perp_model = :chi_perp_i` and `:D_perp` (default `nothing`).
+  - `tau_E` -- whole-plasma energy confinement time τ̄_E [s], required by `P_perp_model = :tau_E`
+    (default `nothing`). That mode takes the minor radius `a` of the ψ = 1 surface from the
+    `:flux` label, whatever `rs_method` is.
+  - `P_perp_model` -- how `P_perp` is built; see [`slayer_parameters`](@ref)
+    (default `:D_perp`, as in the `[SLAYER]` deck).
   - `chi_tor`   -- toroidal heat diffusivity [m²/s]. Scalar or a callable
     of `psi` (default `1.0`).
 
@@ -247,7 +252,10 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles;
     mu_i::Real=2.0,
     zeff=1.0,
     z_i::Real=1.0,
-    chi_perp=1.0,
+    chi_perp_e=1.0,
+    chi_perp_i=nothing,
+    tau_E::Union{Real,Nothing}=nothing,
+    P_perp_model::Symbol=:D_perp,
     chi_tor=1.0,
     dr_val=nothing,
     dgeo_val=nothing,
@@ -256,7 +264,7 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles;
     resistivity_model::NeoResistivityModel=SauterNeoModel(),
     lnLambda_form::Symbol=:nrl)
     R0_use = R0 === nothing ? equil.ro : Float64(R0)
-    _eval(x, ψ) = x isa Real ? Float64(x) : Float64(x(ψ))
+    _eval(x, ψ) = x === nothing ? nothing : x isa Real ? Float64(x) : Float64(x(ψ))
 
     # Compute physical B_T = F(ψ) / (2π·R₀) per surface from the F spline
     # when `bt` is not explicitly supplied.
@@ -270,19 +278,19 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles;
         end
 
     _rs_at, _da_dpsi_at = radial_label(equil; rs_method=rs_method, theta=theta, q_beyond=q_beyond)
+    # Boundary minor radius on the flux label, needed only to turn a whole-plasma τ̄_E into a uniform χ_E.
+    a_minor = P_perp_model === :tau_E ? first(radial_label(equil; rs_method=:flux, theta=theta, q_beyond=q_beyond))(1.0) : nothing
 
     # Per-surface ω_*e, ω_*i (diamagnetic frequencies) from spline derivatives.
     # In flux coordinates ω_* = n·(dp/dψ)/(e·n_e), so the surface's toroidal mode
     # number enters here. Main-ion density is taken equal to the electron density
     # (quasi-neutrality, matching the staging step).
     chi1 = 2π * equil.psio
+    _kinetic_at(ψ) = (n_e=Float64(profiles.n_e(ψ)), dn_e=Float64(profiles.n_e(ψ; deriv=DerivOp(1))),
+        T_e=Float64(profiles.T_e(ψ)), dT_e=Float64(profiles.T_e(ψ; deriv=DerivOp(1))),
+        T_i=Float64(profiles.T_i(ψ)), dT_i=Float64(profiles.T_i(ψ; deriv=DerivOp(1))))
     _omega_star_at(ψ, n_tor) = begin
-        n_e = Float64(profiles.n_e(ψ))
-        dn_e = Float64(profiles.n_e(ψ; deriv=DerivOp(1)))
-        T_e = Float64(profiles.T_e(ψ))
-        dT_e = Float64(profiles.T_e(ψ; deriv=DerivOp(1)))
-        T_i = Float64(profiles.T_i(ψ))
-        dT_i = Float64(profiles.T_i(ψ; deriv=DerivOp(1)))
+        (; n_e, dn_e, T_e, dT_e, T_i, dT_i) = _kinetic_at(ψ)
         ω_star_e = n_tor * (2π / chi1) * (T_e * dn_e / n_e + dT_e)
         ω_star_i = -n_tor * (2π / (Float64(z_i) * chi1)) * (T_i * dn_e / n_e + dT_i)
         return (ω_star_e, ω_star_i)
@@ -305,6 +313,7 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles;
         n_res = sing.n[1]
 
         prof = profiles(psi)
+        kin = _kinetic_at(psi)
         # ω_*e, ω_*i from the density and temperature spline derivatives, at this surface's n.
         ω_e_use, ω_i_use = _omega_star_at(psi, n_res)
 
@@ -392,8 +401,12 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles;
             omega_e=ω_e_use, omega_i=ω_i_use,
             qval=q, sval_r=sval_r, bt=_bt_at(psi),
             rs=rs, R0=R0_use, mu_i=mu_i, zeff=_eval(zeff, psi),
-            chi_perp=_eval(chi_perp, psi),
+            chi_perp_e=_eval(chi_perp_e, psi),
             chi_tor=_eval(chi_tor, psi),
+            chi_perp_i=_eval(chi_perp_i, psi),
+            dlnn=kin.dn_e / kin.n_e, dlnTe=kin.dT_e / kin.T_e, dlnTi=kin.dT_i / kin.T_i,
+            tau_E=tau_E, a_minor=a_minor,
+            P_perp_model=P_perp_model,
             m=m_res, n=n_res,
             dr_val=dr_val_k,
             dgeo_val=dgeo_val_k,
diff --git a/src/InnerLayer/SLAYER/LayerParameters.jl b/src/InnerLayer/SLAYER/LayerParameters.jl
index 2b0aa066d..a8d34110b 100644
--- a/src/InnerLayer/SLAYER/LayerParameters.jl
+++ b/src/InnerLayer/SLAYER/LayerParameters.jl
@@ -16,7 +16,7 @@
 
 Dimensionless layer-physics parameters at one rational surface for the
 Fitzpatrick two-fluid drift-MHD SLAYER inner-layer model (Fitzpatrick
-2023; Park et al. 2022), plus dimensional auxiliaries required for
+2023; Park 2022), plus dimensional auxiliaries required for
 de-normalization. The parametrization uses `P_perp`, `P_tor`, and
 `D_norm` (not the older `pr`/`pe`/`ds` set).
 
@@ -28,7 +28,7 @@ de-normalization. The parametrization uses `P_perp`, `P_tor`, and
 | `lu`       | Lundquist number S = τ_R / τ_H                                    |
 | `c_beta`   | Compressibility √(β_local / (1 + β_local))                        |
 | `D_norm`   | (d_β/r_s) · S^(1/3) · √ι_e  (Fitzpatrick normalized scale)        |
-| `P_perp`   | Perpendicular Prandtl number τ_R / τ_⊥                            |
+| `P_perp`   | Perpendicular Prandtl number τ_R / τ_⊥ (set by `P_perp_model`)    |
 | `P_tor`    | Toroidal-direction Prandtl number τ_R / τ_‖tor                    |
 | `Q_e`      | Normalized electron diamagnetic: −tauk · ω_*e                     |
 | `Q_i`      | Normalized ion diamagnetic:      −tauk · ω_*i                     |
@@ -104,6 +104,9 @@ end
 # `dc_tmp = 0` branch.
 const ALLOWED_DC_TYPES = (:none, :lar, :fitzpatrick, :toroidal)
 
+# Allowed P_perp_model values selecting how P_perp is built. `:chi_perp_e` is the default.
+const ALLOWED_P_PERP_MODELS = (:chi_perp_e, :chi_perp_i, :P_phi, :D_perp, :c_beta, :tau_E)
+
 """
     r_based_shear(rs, q, dq_dpsi, da_dpsi) -> Float64
 
@@ -134,7 +137,7 @@ end
 # large-aspect-ratio offsets are Connor, Ham, Hastie & Liu, PPCF 57, 065001 (2015) Eqs. (59) and (61),
 # after Lütjens, Luciani & Garbet, Phys. Plasmas 8, 4267 (2001). Returns dc_tmp as a Float64.
 function _solve_dc_tmp(; dc_type::Symbol, dr_val::Real, dgeo_val::Real,
-    chi_perp::Real, t_e::Real, zeff::Real, tau_ee::Real,
+    chi_perp_e::Real, t_e::Real, zeff::Real, tau_ee::Real,
     rs::Real, R0::Real, sval_r::Real, n_tor::Integer,
     max_iter::Integer=100, tol::Real=1e-10)
     dc_type in ALLOWED_DC_TYPES ||
@@ -151,7 +154,7 @@ function _solve_dc_tmp(; dc_type::Symbol, dr_val::Real, dgeo_val::Real,
         chi_par_lmfp = (2.0 * R0 * vte) / (sqrt(π) * n_tor * abs(sval_r) * Wd)
         chi_par = (chi_par_smfp * chi_par_lmfp) /
                   (chi_par_smfp + chi_par_lmfp)
-        Wd_new = sqrt(8.0) * (chi_perp / chi_par)^0.25 *
+        Wd_new = sqrt(8.0) * (chi_perp_e / chi_par)^0.25 *
                  (1.0 / sqrt((rs / R0) * abs(sval_r) * n_tor))
         if abs(Wd_new - Wd) / max(abs(Wd), 1e-30) < tol
             Wd = Wd_new
@@ -167,13 +170,13 @@ function _solve_dc_tmp(; dc_type::Symbol, dr_val::Real, dgeo_val::Real,
 
     if dc_type === :lar
         return 0.5 * (-dr_val) * π^1.5 *
-               (chi_par / chi_perp)^0.25 *
+               (chi_par / chi_perp_e)^0.25 *
                sqrt((n_tor * abs(sval_r)) / (R0 * rs))
     elseif dc_type === :fitzpatrick
         return -(sqrt(2.0) * π^1.5 * dr_val) / Wd
     elseif dc_type === :toroidal
         return 0.5 * (-dr_val) * π^1.5 *
-               (chi_par / chi_perp)^0.25 * dgeo_val
+               (chi_par / chi_perp_e)^0.25 * dgeo_val
     end
     return 0.0
 end
@@ -181,8 +184,11 @@ end
 """
     slayer_parameters(; n_e, t_e, t_i, omega_e, omega_i,
                         qval, sval_r, bt, rs, R0, mu_i, zeff,
-                        chi_perp, chi_tor,
+                        chi_perp_e, chi_tor,
                         m, n,
+                        chi_perp_i=nothing, tau_E=nothing, a_minor=nothing,
+                        dlnn=nothing, dlnTe=nothing, dlnTi=nothing,
+                        P_perp_model=:D_perp,
                         dr_val=0.0, dgeo_val=0.0,
                         dc_type=:none, ising=0,
                         resistivity_model=SauterNeoModel(),
@@ -211,9 +217,29 @@ parametrization (P_perp/P_tor/D_norm; the older magnetic/electron Prandtl
   - `R0`      -- major radius [m]
   - `mu_i`    -- ion mass in proton-mass units (e.g. 2.0 for D)
   - `zeff`    -- effective charge
-  - `chi_perp`, `chi_tor` -- perpendicular / toroidal heat diffusivity [m²/s]. `chi_perp` feeds both
-    `P_perp` (which wants the one-fluid energy diffusivity χ_E) and the critical-Δ `W_d` (which wants
-    the electron channel); one profile serves both for now, see the stability docs
+  - `chi_perp_e`, `chi_tor` -- electron perpendicular heat / toroidal momentum diffusivity [m²/s]
+  - `chi_perp_i` -- ion perpendicular heat diffusivity [m²/s], required 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 τ̄_E [s], required by `P_perp_model=:tau_E`
+  - `a_minor` -- plasma minor radius `a` [m], required by `P_perp_model=:tau_E` (`build_slayer_inputs` takes it from the flux label)
+
+  An input a model does not use may be left as `nothing`; one it uses must be given.
+  - `P_perp_model` -- how `P_perp` is built (default `:D_perp`, as in the `[SLAYER]` deck):
+      + `:D_perp` -- τ_R/τ_⊥ with τ_⊥ = r_s²/D_⊥, D_⊥ the perpendicular particle diffusivity of
+        Fitzpatrick, Phys. Plasmas 29, 032507 (2022), Eq. 16 (with Eqs. 4-6 and η_⊥ = η_∥). Its
+        τ = (T_e/T_i)(1+η_e)/(1+η_i) is the inverse orientation of the `tau` = T_i/T_e field
+        of `SLAYERParameters`. With flat temperatures it falls to the `:c_beta` value
+      + `:chi_perp_e` -- τ_R·χ⊥,e/r_s²
+      + `:chi_perp_i` -- τ_R·χ⊥,i/r_s²
+      + `:P_phi`  -- P_perp = P_tor (Fitzpatrick, Phys. Plasmas 30, 092512 (2023), Eq. 143)
+      + `:c_beta` -- P_perp = c_β² (J.-K. Park, Phys. Plasmas 29, 072506 (2022), Eq. 9 with K = 0),
+        the χ⊥ → 0 limit of `:D_perp`
+      + `:tau_E`  -- P_perp = P_tor = τ_R/τ_E(r_s) (Fitzpatrick 2023, Eqs. 37-38, 81-83 and 131), from
+        a uniform χ_E = κ_e a²/τ̄_E with κ_e = 3/16, the value for T_e ∝ 1 − r̂² (Eq. 38). Other
+        profile shapes change κ_e by a factor of order 2-3. Both Prandtl numbers come from χ_E, so
+        `chi_tor` is ignored in this mode
+    `chi_perp_e` still sets the critical-Δ `dc_tmp` in every mode.
   - `m`, `n`  -- poloidal / toroidal mode numbers at the surface
   - `dr_val`, `dgeo_val` -- inputs for the critical-Δ formula
   - `dc_type` -- one of `:none`, `:lar`, `:fitzpatrick`, `:toroidal`
@@ -266,8 +292,11 @@ function slayer_parameters(;
     omega_e::Real, omega_i::Real,
     qval::Real, sval_r::Real, bt::Real,
     rs::Real, R0::Real, mu_i::Real, zeff::Real,
-    chi_perp::Real, chi_tor::Real,
+    chi_perp_e::Real, chi_tor::Real,
     m::Integer, n::Integer,
+    chi_perp_i::Union{Real,Nothing}=nothing, tau_E::Union{Real,Nothing}=nothing, a_minor::Union{Real,Nothing}=nothing,
+    dlnn::Union{Real,Nothing}=nothing, dlnTe::Union{Real,Nothing}=nothing, dlnTi::Union{Real,Nothing}=nothing,
+    P_perp_model::Symbol=:D_perp,
     dr_val::Real=0.0, dgeo_val::Real=0.0,
     dc_type::Symbol=:none, ising::Integer=0,
     resistivity_model::NeoResistivityModel=SauterNeoModel(),
@@ -368,11 +397,46 @@ function slayer_parameters(;
     lbeta = (5.0 / 3.0) * MU_0 * n_e * E_CHG * (t_e + t_i) / bt^2
     c_beta = sqrt(lbeta / (1.0 + lbeta))
 
-    # Effective Prandtl-like transport ratios
-    tau_perp = rs^2 / chi_perp
-    P_perp = tau_r / tau_perp
+    # Effective Prandtl-like transport ratios. A diffusivity the model uses must be given, finite and positive.
+    _need(name, v) = (v !== nothing && isfinite(v) && v > 0) ||
+                     throw(ArgumentError("slayer_parameters: P_perp_model=:$P_perp_model needs a positive, finite $name; got $v"))
+    _need("chi_perp_e", chi_perp_e)
+    P_perp_model === :tau_E || _need("chi_tor", chi_tor)
+    P_perp_model in (:chi_perp_i, :D_perp) && _need("chi_perp_i", chi_perp_i)
     tau_tor = rs^2 / chi_tor
     P_tor = tau_r / tau_tor
+    if P_perp_model === :tau_E
+        _need("tau_E", tau_E)
+        _need("a_minor", a_minor)
+        # Fitzpatrick 2023, Eqs. 37-38, 81-83 and 131.
+        chi_E = (3 / 16) * a_minor^2 / tau_E
+        P_tor = tau_r / (rs^2 / chi_E)
+    end
+    P_perp = if P_perp_model === :chi_perp_e
+        tau_r / (rs^2 / chi_perp_e)
+    elseif P_perp_model === :chi_perp_i
+        tau_r / (rs^2 / chi_perp_i)
+    elseif P_perp_model === :P_phi || P_perp_model === :tau_E
+        P_tor
+    elseif P_perp_model === :D_perp
+        # Fitzpatrick, Phys. Plasmas 29, 032507 (2022), Eqs. 4-6 and 16 with η_⊥ = η_∥, written over the pressure gradient p′/n.
+        all(x -> x !== nothing && isfinite(x), (dlnn, dlnTe, dlnTi)) ||
+            throw(ArgumentError("slayer_parameters: P_perp_model=:D_perp needs finite dlnn, dlnTe and dlnTi; got $dlnn, $dlnTe, $dlnTi"))
+        dp = t_e * (dlnn + dlnTe) + t_i * (dlnn + dlnTi)
+        dp != 0 || throw(ArgumentError("slayer_parameters: P_perp_model=:D_perp is undefined where the pressure gradient vanishes " *
+                                       "(t_e(dlnn+dlnTe) + t_i(dlnn+dlnTi) = 0); choose another P_perp_model for this surface"))
+        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
+        tau_r / (rs^2 / D_perp)
+    elseif P_perp_model === :c_beta
+        c_beta^2
+    else
+        throw(ArgumentError("slayer_parameters: unknown P_perp_model=$P_perp_model. " *
+                            "Allowed: $(ALLOWED_P_PERP_MODELS)"))
+    end
+    isfinite(P_perp) && P_perp > 0 ||
+        throw(ArgumentError("slayer_parameters: P_perp_model=:$P_perp_model gave P_perp=$P_perp; check its inputs " *
+                            "(chi_perp_e=$chi_perp_e, chi_perp_i=$chi_perp_i, dlnn=$dlnn, dlnTe=$dlnTe, dlnTi=$dlnTi)"))
 
     # Normalized beta-related width and Δ-normalization
     d_beta = c_beta * d_i
@@ -381,7 +445,7 @@ function slayer_parameters(;
 
     # Critical-Δ offset from chi_parallel matching
     dc_tmp = _solve_dc_tmp(; dc_type=dc_type, dr_val=dr_val, dgeo_val=dgeo_val,
-        chi_perp=chi_perp, t_e=t_e, zeff=zeff,
+        chi_perp_e=chi_perp_e, t_e=t_e, zeff=zeff,
         tau_ee=tau_ee, rs=rs, R0=R0, sval_r=sval_r,
         n_tor=n)
 
diff --git a/src/InnerLayer/SLAYER/LayerThickness.jl b/src/InnerLayer/SLAYER/LayerThickness.jl
index a07d4ea48..9e55b178d 100644
--- a/src/InnerLayer/SLAYER/LayerThickness.jl
+++ b/src/InnerLayer/SLAYER/LayerThickness.jl
@@ -2,7 +2,7 @@
 #
 # Resistive inner-layer *width* (del_s) diagnostic for the SLAYER two-fluid
 # drift-MHD slab layer (Fitzpatrick, Tearing Mode Dynamics in Tokamak
-# Plasmas, IOP 2023; Park et al. 2022, Phys. Plasmas 29, 122505; Burgess
+# Plasmas, IOP 2023; Park 2022, Phys. Plasmas 29, 072506; Burgess
 # et al. 2026). This is a distinct quantity from the matching index
 # Δ = π/W' returned by the dispersion solver (Riccati.jl): it integrates a
 # separate reduced Riccati ODE evaluated at the electron diamagnetic
diff --git a/src/InnerLayer/SLAYER/Riccati.jl b/src/InnerLayer/SLAYER/Riccati.jl
index 8f4e518c6..22c57e51c 100644
--- a/src/InnerLayer/SLAYER/Riccati.jl
+++ b/src/InnerLayer/SLAYER/Riccati.jl
@@ -2,7 +2,7 @@
 #
 # Inner-layer Δ via the Fitzpatrick two-fluid drift-MHD layer Riccati ODE
 # (Fitzpatrick, Tearing Mode Dynamics in Tokamak Plasmas, IOP 2023; Park
-# et al. 2022, Phys. Plasmas 29, 122505; Burgess et al. 2026), for the
+# 2022, Phys. Plasmas 29, 072506; Burgess et al. 2026), for the
 # pressureless tearing channel (no parallel flow, pe = 0). The A, B, C
 # coefficients, the dW/dp Riccati equation, the large-p asymptotic boundary
 # conditions, the regime test D² > ι_e P_⊥/P_tor^(2/3), the small-p
@@ -169,7 +169,7 @@ pressureless layer produces only the tearing channel.
 # Algorithm
 
 Implements the Fitzpatrick two-fluid drift-MHD layer formulation
-(Fitzpatrick 2023; Park et al. 2022) for the pressureless tearing channel
+(Fitzpatrick 2023; Park 2022) for the pressureless tearing channel
 (no parallel flow, pe = 0). Integrates
 `dW/dp = -(fA'/p)·W − W²/p + (fB/(fA·fC))·p³` from a large `p_start`
 (selected by `_riccati_f_initial` according to whether
diff --git a/src/InnerLayer/SLAYER/SLAYER.jl b/src/InnerLayer/SLAYER/SLAYER.jl
index 7f6c95c00..efc81cd31 100644
--- a/src/InnerLayer/SLAYER/SLAYER.jl
+++ b/src/InnerLayer/SLAYER/SLAYER.jl
@@ -1,8 +1,8 @@
 # SLAYER.jl
 #
 # SLAYER (Slab Layer) two-fluid drift-MHD inner-layer model (Fitzpatrick,
-# Tearing Mode Dynamics in Tokamak Plasmas, IOP 2023; Park et al. 2022,
-# Phys. Plasmas 29, 122505; Burgess et al. 2026). The dispersion path uses
+# Tearing Mode Dynamics in Tokamak Plasmas, IOP 2023; Park 2022,
+# Phys. Plasmas 29, 072506; Burgess et al. 2026). The dispersion path uses
 # the Fitzpatrick layer formulation — P_perp / P_tor transport, c_beta
 # compressibility, D_norm normalized ion-skin scale, two-fluid drift
 # coupling via Q_e, Q_i, iota_e (see `Riccati.jl`).
diff --git a/src/Tearing/LayerOverlap.jl b/src/Tearing/LayerOverlap.jl
index 64c187d2d..4fc96d21e 100644
--- a/src/Tearing/LayerOverlap.jl
+++ b/src/Tearing/LayerOverlap.jl
@@ -155,7 +155,7 @@ Riccati solve; `delta_s` is solved per surface as a reported comparison only.
   - `rs_method` -- radial label the widths and their ψ-conversion Jacobian are taken in
     (default `:midplane`; see `radial_label`). The caller wiring this scan into the domain
     truncation passes `:flux`, the label Fitzpatrick Eq. (100) is anchored to
-  - `bt`, `mu_i`, `zeff`, `chi_perp`, `chi_tor`, `resistivity_model`, `lnLambda_form`
+  - `bt`, `mu_i`, `zeff`, `chi_perp_e`, `chi_perp_i`, `chi_tor`, `tau_E`, `P_perp_model`, `resistivity_model`, `lnLambda_form`
     -- passed through to `build_slayer_inputs`. The defaults are placeholders: a caller
     scanning a real run must pass that run's own values, as the domain-truncation wiring does
 
@@ -176,8 +176,11 @@ function resistive_layer_overlap(equil, profiles::KineticProfiles;
     bt=nothing,
     mu_i::Real=2.0,
     zeff=1.0,
-    chi_perp=1.0,
+    chi_perp_e=1.0,
+    chi_perp_i=nothing,
     chi_tor=1.0,
+    tau_E::Union{Real,Nothing}=nothing,
+    P_perp_model::Symbol=:D_perp,
     resistivity_model::NeoResistivityModel=SauterNeoModel(),
     lnLambda_form::Symbol=:nrl)
 
@@ -316,7 +319,8 @@ function resistive_layer_overlap(equil, profiles::KineticProfiles;
         # which makes iota_e singular) is recorded and skipped rather than aborting the scan.
         params = try
             build_slayer_inputs(equil, [sing], profiles; rs_method=rs_method, q_beyond=q_beyond,
-                bt=bt, mu_i=mu_i, zeff=zeff, chi_perp=chi_perp, chi_tor=chi_tor,
+                bt=bt, mu_i=mu_i, zeff=zeff, chi_perp_e=chi_perp_e, chi_perp_i=chi_perp_i, chi_tor=chi_tor,
+                tau_E=tau_E, P_perp_model=P_perp_model,
                 dr_val=0.0, dc_type=:none, theta=theta,
                 resistivity_model=resistivity_model, lnLambda_form=lnLambda_form)
         catch err
diff --git a/src/Tearing/Runner/Control.jl b/src/Tearing/Runner/Control.jl
index 1e9c2ba2e..493383660 100644
--- a/src/Tearing/Runner/Control.jl
+++ b/src/Tearing/Runner/Control.jl
@@ -46,10 +46,17 @@ constructor.
   - `impurity_charge` -- charge of the single impurity species that closes quasi-neutrality when
     `Z_eff(ψ)` is derived from `n_i/n_e` (default 6.0, carbon); the main ion has charge 1.
     Unused when the file supplies a `zeff` profile
-  - `chi_perp`, `chi_tor` -- fallback perpendicular / toroidal heat
-    diffusivity [m²/s], used only when the kinetic file carries no usable
-    `chi_e`/`chi_phi` profile (dataset absent or all-zero); otherwise the
-    file's χ⊥(ψ)/χ_φ(ψ) take precedence
+  - `chi_perp_e`, `chi_perp_i`, `chi_tor` -- fallback electron / ion perpendicular heat and
+    toroidal momentum diffusivity [m²/s], used only when the kinetic file carries no usable
+    `chi_e`/`chi_i`/`chi_phi` profile (dataset absent or all-zero); otherwise the file's
+    χ⊥,e(ψ)/χ⊥,i(ψ)/χ_φ(ψ) take precedence. `chi_perp_i` defaults to `nothing`, so a model that
+    needs χ⊥,i fails loudly when neither the file nor the deck supplies it. Each must be positive
+  - `tau_E`    -- whole-plasma energy confinement time τ̄_E [s], required by
+    `P_perp_model = :tau_E` (default `nothing`)
+  - `P_perp_model` -- how the perpendicular Prandtl number is built: `:D_perp` (default; needs
+    χ⊥,i), `:chi_perp_e`, `:chi_perp_i`, `:P_phi`, `:c_beta` or `:tau_E`; see
+    `InnerLayer.SLAYER.slayer_parameters`. `:tau_E` sets both P_perp and P_tor from τ̄_E, so it
+    ignores `chi_phi`/`chi_tor`. `chi_perp_e` still sets the critical-Δ width in every mode
   - `D_R`, `D_geo`  -- critical-Δ formula inputs. `nothing` (default)
     auto-derives them from the equilibrium: `D_R` from the resistive
     interchange index `D_R = E + F + H²` at each surface, `D_geo` from the
@@ -96,9 +103,9 @@ there is one consistent interface for resistive and kinetic profiles.
 
   - `profile_file`   -- path to a kinetic-profile file (relative to the run
     dir), required when SLAYER is enabled. HDF5 (`.h5`) files use the GPEC
-    kinetic schema and may carry `chi_e` (χ⊥) and `chi_phi` (χ_φ); ASCII
+    kinetic schema and may carry `chi_e` (χ⊥,e), `chi_i` (χ⊥,i) and `chi_phi` (χ_φ); ASCII
     tables are also accepted but carry no χ. When a χ dataset is absent or
-    all-zero, the scalar `chi_perp`/`chi_tor` fallbacks below are used (so a
+    all-zero, the scalar `chi_perp_e`/`chi_perp_i`/`chi_tor` fallbacks below are used (so a
     file can keep the χ keys set to 0 to defer to the scalars).
   - `profile_group`  -- group within the HDF5 file (default `"/"`)
 
@@ -120,8 +127,11 @@ there is one consistent interface for resistive and kinetic profiles.
     ion_mass::Float64 = 2.0
     zeff::Float64 = 1.0
     impurity_charge::Float64 = 6.0
-    chi_perp::Float64 = 1.0
+    chi_perp_e::Float64 = 1.0
+    chi_perp_i::Union{Float64,Nothing} = nothing
     chi_tor::Float64 = 1.0
+    tau_E::Union{Float64,Nothing} = nothing
+    P_perp_model::Symbol = :D_perp
     D_R::Union{Float64,Nothing} = nothing
     D_geo::Union{Float64,Nothing} = nothing
     theta_sample::Float64 = 0.0
@@ -177,6 +187,7 @@ const _UNSUPPORTED_INNER_MODELS = (:ggj_ray, :ggj_shooting, :ggj_galerkin)
 const _VALID_SCAN_MODES = (:amr, :brute_force)
 const _VALID_COUPLING_MODES = (:uncoupled, :coupled)
 const _VALID_D_C_TYPES = (:none, :lar, :fitzpatrick, :toroidal)
+const _VALID_P_PERP_MODELS = (:chi_perp_e, :chi_perp_i, :P_phi, :D_perp, :c_beta, :tau_E)
 const _VALID_RESISTIVITY_MODELS = (:sauter, :redl, :spitzer, :spitzer_harm)
 const _VALID_LNLAMBDA_FORMS = (:nrl, :sauter, :wesson)
 
@@ -197,6 +208,15 @@ function validate(ctrl::SLAYERControl)
     ctrl.D_c_type in _VALID_D_C_TYPES ||
         throw(ArgumentError("SLAYERControl: D_c_type=$(ctrl.D_c_type) " *
                             "not in $(_VALID_D_C_TYPES)"))
+    ctrl.P_perp_model in _VALID_P_PERP_MODELS ||
+        throw(ArgumentError("SLAYERControl: P_perp_model=$(ctrl.P_perp_model) " *
+                            "not in $(_VALID_P_PERP_MODELS)"))
+    for (name, v) in (("chi_perp_e", ctrl.chi_perp_e), ("chi_perp_i", ctrl.chi_perp_i), ("chi_tor", ctrl.chi_tor))
+        v === nothing || (isfinite(v) && v > 0) ||
+            throw(ArgumentError("SLAYERControl: $name must be positive and finite [m²/s]; got $v"))
+    end
+    ctrl.P_perp_model !== :tau_E || (ctrl.tau_E !== nothing && isfinite(ctrl.tau_E) && ctrl.tau_E > 0) ||
+        throw(ArgumentError("SLAYERControl: P_perp_model=tau_E requires a positive, finite tau_E [s]"))
     ctrl.resistivity_model in _VALID_RESISTIVITY_MODELS ||
         throw(ArgumentError("SLAYERControl: resistivity_model=$(ctrl.resistivity_model) " *
                             "not in $(_VALID_RESISTIVITY_MODELS)"))
@@ -257,9 +277,10 @@ function slayer_control_from_toml(section::AbstractDict)
     # Validate keys against the struct fields
     field_names = Set(String.(fieldnames(SLAYERControl)))
     unknown = [k for k in keys(flat) if !(k in field_names)]
+    hint = "chi_perp" in unknown ? " (`chi_perp` is now `chi_perp_e`.)" : ""
     isempty(unknown) ||
         throw(ArgumentError("slayer_control_from_toml: unknown keys " *
-                            "$(unknown) in [SLAYER] section. Known: " *
+                            "$(unknown) in [SLAYER] section.$hint Known: " *
                             "$(sort(collect(field_names)))."))
 
     # Coerce types where needed
@@ -267,11 +288,11 @@ function slayer_control_from_toml(section::AbstractDict)
     for (k, v) in flat
         sym = Symbol(k)
         if sym in (:inner_model, :scan_mode, :coupling_mode, :D_c_type,
-            :resistivity_model, :lnLambda_form)
+            :resistivity_model, :lnLambda_form, :P_perp_model)
             kwargs[sym] = v isa Symbol ? v : Symbol(String(v))
         elseif sym in (:Q_re_range, :Q_im_range)
             kwargs[sym] = _as_range(v)
-        elseif sym in (:bt, :D_R, :D_geo)
+        elseif sym in (:bt, :D_R, :D_geo, :tau_E, :chi_perp_i)
             # Allow explicit nothing (auto-derive) or a number (override)
             kwargs[sym] = v === nothing ? nothing : Float64(v)
         elseif sym === :omega_E_kHz
diff --git a/src/Tearing/Runner/TearingProblem.jl b/src/Tearing/Runner/TearingProblem.jl
index 583e13453..a5c9aeb10 100644
--- a/src/Tearing/Runner/TearingProblem.jl
+++ b/src/Tearing/Runner/TearingProblem.jl
@@ -33,12 +33,8 @@ function _tearing_params(::SLAYERModel, equil, surfaces, loaded, control)
     # physical B_T = F(psi)/(2*pi*R_0) per surface from the equilibrium's F-spline, which is
     # what its docstring already prescribes.
     bt = control.bt
-    # χ⊥/χ_φ from the kinetic file when present, else the scalar fallbacks.
-    chi_perp = loaded.chi_perp === nothing ? control.chi_perp : loaded.chi_perp
-    chi_tor = loaded.chi_tor === nothing ? control.chi_tor : loaded.chi_tor
-    (loaded.chi_perp === nothing || loaded.chi_tor === nothing) && @warn(
-        "SLAYER: no usable chi_e/chi_phi profile(s) (dataset absent or all-zero); " *
-        "using the scalar chi_perp/chi_tor fallback for the missing one(s).")
+    # χ⊥,e/χ⊥,i/χ_φ from the kinetic file when present, else the scalar fallbacks.
+    chi = _slayer_chi_inputs(loaded, control)
     # Z_eff(ψ) from the kinetic file when it carries n_i, else the scalar fallback.
     zeff = loaded.zeff === nothing ? control.zeff : loaded.zeff
     if loaded.zeff_source === :scalar
@@ -52,8 +48,11 @@ function _tearing_params(::SLAYERModel, equil, surfaces, loaded, control)
         bt=bt,
         mu_i=control.ion_mass,
         zeff=zeff,
-        chi_perp=chi_perp,
-        chi_tor=chi_tor,
+        chi_perp_e=chi.chi_perp_e,
+        chi_perp_i=chi.chi_perp_i,
+        chi_tor=chi.chi_tor,
+        tau_E=control.tau_E,
+        P_perp_model=control.P_perp_model,
         dr_val=control.D_R,
         dgeo_val=control.D_geo,
         dc_type=control.D_c_type,
diff --git a/src/Tearing/Runner/run_slayer.jl b/src/Tearing/Runner/run_slayer.jl
index 01f421173..72eed83f4 100644
--- a/src/Tearing/Runner/run_slayer.jl
+++ b/src/Tearing/Runner/run_slayer.jl
@@ -8,9 +8,9 @@
 # ---------------------------------------------------------------------
 # Read kinetic profiles through the shared `Equilibrium.read_kinetic_file`
 # reader and adapt them to the SLAYER inner layer. Returns the spline-based
-# `KineticProfiles` plus χ⊥(ψ)/χ_φ(ψ) callables built from the file's
-# `chi_e`/`chi_phi` (or `nothing` when the file carries no usable χ, in which
-# case the caller falls back to the scalar `control.chi_perp`/`chi_tor`), and
+# `KineticProfiles` plus χ⊥,e(ψ)/χ⊥,i(ψ)/χ_φ(ψ) callables built from the file's
+# `chi_e`/`chi_i`/`chi_phi` (or `nothing` when the file carries no usable χ, in which
+# case the caller falls back to the scalar `control.chi_perp_e`/`chi_perp_i`/`chi_tor`), and
 # a Z_eff(ψ) callable: the file's own `zeff` profile when it carries one, else
 # derived from its `n_i/n_e` (or `nothing` when the file carries neither, in which
 # case the caller falls back to the scalar `control.zeff`), and whether the file
@@ -38,13 +38,19 @@ function _load_profiles(control::SLAYERControl, dir_path::AbstractString)
     # χ⊥(ψ)/χ_φ(ψ) splines from the file. A χ array that is absent OR all-zero
     # is treated as "not provided" — χ must be positive (χ=0 ⇒ τ_⊥→∞), and the
     # all-zero sentinel lets a file keep the chi_e/chi_phi keys while deferring
-    # to the scalar control.chi_perp/chi_tor fallback. Build the spline only
+    # to the scalar control.chi_perp_e/chi_tor fallback. Build the spline only
     # from a usable (present, not all-zero) array.
     psi_xs = collect(Float64, data.psi)
-    _chi_spline(v) = (v === nothing || all(iszero, v)) ? nothing :
-                     cubic_interp(psi_xs, collect(Float64, v))
-    chi_perp = _chi_spline(data.chi_e)
-    chi_tor = _chi_spline(data.chi_phi)
+    function _chi_spline(name, v)
+        (v === nothing || all(iszero, v)) && return nothing
+        minimum(v) > 0 || error("run_slayer: kinetic file '$path' has $name with minimum $(minimum(v)) ≤ 0; diffusivities must be positive")
+        return cubic_interp(psi_xs, collect(Float64, v))
+    end
+    model = control.P_perp_model
+    uses_e = model in (:chi_perp_e, :D_perp) || control.D_c_type !== :none
+    chi_perp_e = uses_e ? _chi_spline("chi_e", data.chi_e) : nothing
+    chi_perp_i = model in (:chi_perp_i, :D_perp) ? _chi_spline("chi_i", data.chi_i) : nothing
+    chi_tor = model !== :tau_E ? _chi_spline("chi_phi", data.chi_phi) : nothing
     # Z_eff(ψ): a profile supplied in the file wins, since a real impurity mix is not recoverable
     # from n_i/n_e. Otherwise quasi-neutrality with one impurity species of charge `impurity_charge` and a
     # charge-1 main ion, the single-ion closure the kinetic spline set uses:
@@ -63,8 +69,39 @@ function _load_profiles(control::SLAYERControl, dir_path::AbstractString)
         cubic_interp(psi_xs, z)
     end
     zeff_source = _usable(data.zeff) ? :file : zeff === nothing ? :scalar : :closure
-    return (profiles=profiles, chi_perp=chi_perp, chi_tor=chi_tor, zeff=zeff, zeff_source=zeff_source,
-        has_omega_E=data.omega_E !== nothing)
+    return (profiles=profiles, chi_perp_e=chi_perp_e, chi_perp_i=chi_perp_i, chi_tor=chi_tor, zeff=zeff, zeff_source=zeff_source,
+        has_omega_E=(data.omega_E !== nothing))
+end
+
+# χ⊥,e/χ⊥,i/χ_φ for the layer: the kinetic file's profiles where present, else the control's
+# scalar fallbacks. With `warn`, a fallback is reported only for a χ the selected P_perp model
+# (and, for χ⊥,e, the critical Δ) actually uses.
+function _slayer_chi_inputs(loaded, control::SLAYERControl; warn::Bool=true)
+    model = control.P_perp_model
+    pick(f, s) = f === nothing ? s : f
+    model in (:chi_perp_i, :D_perp) && loaded.chi_perp_i === nothing && control.chi_perp_i === nothing &&
+        throw(
+            ArgumentError(
+                "SLAYER: P_perp_model = $model needs the ion diffusivity χ⊥,i: add a chi_i profile to the " *
+                "kinetic file, set [SLAYER] chi_perp_i, or choose P_perp_model = \"chi_perp_e\"."
+            )
+        )
+    if warn
+        uses_e = model in (:chi_perp_e, :D_perp) || control.D_c_type !== :none
+        uses_i = model in (:chi_perp_i, :D_perp)
+        uses_tor = model !== :tau_E
+        fallbacks = String[]
+        uses_e && loaded.chi_perp_e === nothing && push!(fallbacks, "chi_e → chi_perp_e = $(control.chi_perp_e)")
+        uses_i && loaded.chi_perp_i === nothing && push!(fallbacks, "chi_i → chi_perp_i = $(control.chi_perp_i)")
+        uses_tor && loaded.chi_tor === nothing && push!(fallbacks, "chi_phi → chi_tor = $(control.chi_tor)")
+        isempty(fallbacks) || @warn "SLAYER: no usable kinetic-file profile (dataset absent or all-zero) for " *
+                                    "$(join(fallbacks, ", ")); using the scalar fallback."
+        model === :tau_E && @info "SLAYER: P_perp_model = tau_E sets P_perp and P_tor from a uniform " *
+              "χ_E = (3/16)a²/τ̄_E with τ̄_E = $(control.tau_E) s; the χ_φ profile is not used."
+    end
+    return (; chi_perp_e=pick(loaded.chi_perp_e, control.chi_perp_e),
+        chi_perp_i=pick(loaded.chi_perp_i, control.chi_perp_i),
+        chi_tor=pick(loaded.chi_tor, control.chi_tor))
 end
 
 # ---------------------------------------------------------------------
@@ -337,7 +374,7 @@ function run_slayer_from_inputs(model::InnerLayer.InnerLayerModel,
     # A surface with no known rotation applies no shift.
     Ω_shift = [isnan(Ω) ? 0.0 : Ω for Ω in Ω_E]
     coupled && _is_ggj(model) && any(!iszero, Ω_shift) && @warn(
-            "SLAYER: E×B rotation is not applied to GGJ surfaces, which carry no time normalization.")
+        "SLAYER: E×B rotation is not applied to GGJ surfaces, which carry no time normalization.")
     q_shifts = coupled ? Float64[_q_shift(params[k], Ω_shift[k]) for k in 1:n] : zeros(Float64, n)
     scs = [_build_surface_coupling(model, params[k], dp[k, k], q_shifts[k]) for k in 1:n]
     # Coupled Q is scanned in the reference surface's normalization.
diff --git a/test/runtests.jl b/test/runtests.jl
index bc068435a..9f5fc3aef 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -62,6 +62,7 @@ else
     include("./runtests_dispersion_scan.jl")
     include("./runtests_dispersion_amr.jl")
     include("./runtests_dispersion_polish.jl")
+    include("./runtests_slayer_profile_transport.jl")
     include("./runtests_slayer_runner.jl")
     include("./runtests_kinetic.jl")
     include("./runtests_kinetic_validity.jl")
diff --git a/test/runtests_dispersion_coupled.jl b/test/runtests_dispersion_coupled.jl
index c25857969..456f49624 100644
--- a/test/runtests_dispersion_coupled.jl
+++ b/test/runtests_dispersion_coupled.jl
@@ -20,12 +20,12 @@
         InnerLayerResponse(m.a + m.b * ComplexF64(Q), zero(ComplexF64))
 
     function _slayer_ref()
-        return slayer_parameters(;
+        return slayer_parameters(; P_perp_model=:chi_perp_e,
             n_e=5.0e19, t_e=1000.0, t_i=1000.0,
             omega_e=-1.0e4, omega_i=5.0e3,
             qval=2.0, sval_r=1.0, bt=2.0,
             rs=0.5, R0=1.7, mu_i=2.0, zeff=1.0,
-            chi_perp=1.0, chi_tor=1.0, m=2, n=1)
+            chi_perp_e=1.0, chi_tor=1.0, m=2, n=1)
     end
 
     @testset "Constructor validation" begin
@@ -199,12 +199,12 @@
         # exactly.
         # A second surface with different n_e and T_e, so tauk differs and the rescale is exercised.
         p_a = _slayer_ref()
-        p_b = slayer_parameters(;
+        p_b = slayer_parameters(; P_perp_model=:chi_perp_e,
             n_e=2.0e19, t_e=500.0, t_i=500.0,
             omega_e=-1.0e4, omega_i=5.0e3,
             qval=3.0, sval_r=1.5, bt=2.0,
             rs=0.6, R0=1.7, mu_i=2.0, zeff=1.0,
-            chi_perp=1.0, chi_tor=1.0, m=3, n=1)
+            chi_perp_e=1.0, chi_tor=1.0, m=3, n=1)
         m = SLAYERModel()
         sc1 = surface_coupling(m, p_a, 0.0 + 0im)
         sc2 = surface_coupling(m, p_b, 0.0 + 0im)
diff --git a/test/runtests_dispersion_residual.jl b/test/runtests_dispersion_residual.jl
index 28a8d52f4..4758f725b 100644
--- a/test/runtests_dispersion_residual.jl
+++ b/test/runtests_dispersion_residual.jl
@@ -25,7 +25,7 @@
             omega_e=-1.0e4, omega_i=5.0e3,
             qval=2.0, sval_r=1.0, bt=2.0,
             rs=0.5, R0=1.7, mu_i=2.0, zeff=1.0,
-            chi_perp=1.0, chi_tor=1.0, m=2, n=1)
+            chi_perp_e=1.0, chi_tor=1.0, m=2, n=1; P_perp_model=:chi_perp_e)
     end
 
     @testset "Constructor scale defaults" begin
diff --git a/test/runtests_dispersion_rotation.jl b/test/runtests_dispersion_rotation.jl
index 80094c173..c25efa4ef 100644
--- a/test/runtests_dispersion_rotation.jl
+++ b/test/runtests_dispersion_rotation.jl
@@ -48,10 +48,10 @@
         # Q_layer = τ_k·ω_lab − τ_k·n·Ω_E, so the rotating determinant at Q + τ_ref·n·Ω_E equals the
         # static one at Q: every root keeps its γ and moves by ω_lab(Ω_E) − ω_lab(0) = +n·Ω_E. This pins
         # how `_q_shift` composes with the determinant; the next testset ties the sign to the physical axis.
-        mk(; qval, rs, m) = slayer_parameters(; n_e=5.0e19, t_e=1000.0, t_i=1000.0,
+        mk(; qval, rs, m) = slayer_parameters(; P_perp_model=:chi_perp_e, n_e=5.0e19, t_e=1000.0, t_i=1000.0,
             omega_e=-1.0e4, omega_i=5.0e3,
             qval=qval, sval_r=1.0, bt=2.0, rs=rs, R0=1.7,
-            mu_i=2.0, zeff=1.0, chi_perp=1.0, chi_tor=1.0, m=m, n=2)
+            mu_i=2.0, zeff=1.0, chi_perp_e=1.0, chi_tor=1.0, m=m, n=2)
         params = [mk(; qval=1.5, rs=0.5, m=3), mk(; qval=2.0, rs=0.6, m=4)]
         slayer = SLAYERModel(; variant=:fitzpatrick)
         dp = ComplexF64[-2.0 0.5; 0.5 -3.0]
@@ -75,8 +75,8 @@
         # the `i·conj(Q)` bridge in `solve_inner` moves it to −τ_k·ω_*e.
         slayer = SLAYERModel(; variant=:fitzpatrick)
         for (ω_e, ω_i) in ((-1.0e4, 5.0e3), (1.0e4, -5.0e3), (-1.0e4, 0.0))
-            p = slayer_parameters(; n_e=5.0e19, t_e=1000.0, t_i=1000.0, omega_e=ω_e, omega_i=ω_i,
-                qval=2.0, sval_r=1.0, bt=2.0, rs=0.5, R0=1.7, mu_i=2.0, zeff=1.0, chi_perp=1.0, chi_tor=1.0, m=2, n=1)
+            p = slayer_parameters(; P_perp_model=:chi_perp_e, n_e=5.0e19, t_e=1000.0, t_i=1000.0, omega_e=ω_e, omega_i=ω_i,
+                qval=2.0, sval_r=1.0, bt=2.0, rs=0.5, R0=1.7, mu_i=2.0, zeff=1.0, chi_perp_e=1.0, chi_tor=1.0, m=2, n=1)
             Q_electron = p.tauk * ω_e + 0.01im
             @test abs(solve_inner(slayer, p, Q_electron).tearing) < 0.1 * abs(solve_inner(slayer, p, -conj(Q_electron)).tearing)
         end
diff --git a/test/runtests_kinetic.jl b/test/runtests_kinetic.jl
index 96e2c965b..5c99e003f 100644
--- a/test/runtests_kinetic.jl
+++ b/test/runtests_kinetic.jl
@@ -737,8 +737,9 @@
             da = E.read_kinetic_file(gpeckf)
             db = E.read_kinetic_file(h5)
             # legacy ASCII carries no optional fields; the example .h5 carries chi
-            @test da.chi_e === nothing && da.chi_phi === nothing
-            @test db.chi_e !== nothing && db.chi_phi !== nothing
+            @test da.chi_e === nothing && da.chi_i === nothing && da.chi_phi === nothing
+            @test db.chi_e !== nothing && db.chi_i !== nothing && db.chi_phi !== nothing
+            @test all(>(0), db.chi_i)
             # core columns identical (the .h5 was generated from the .gpeckf)
             @test da.psi == db.psi
             @test da.n_e == db.n_e
@@ -766,6 +767,9 @@
             @test dc.chi_phi ≈ db.chi_phi
             @test dc.n_e ≈ db.n_e
             @test dc.provenance == "roundtrip"
+            @test dc.chi_i ≈ db.chi_i
+            E.write_kinetic_h5(tmp, E.KineticProfileData(; psi=db.psi, chi_e=db.chi_e))
+            @test E.read_kinetic_file(tmp).chi_i === nothing
             rm(tmp; force=true)
         end
 
diff --git a/test/runtests_layer_overlap_fitzpatrick.jl b/test/runtests_layer_overlap_fitzpatrick.jl
index 5de7412f5..2d91b59fc 100644
--- a/test/runtests_layer_overlap_fitzpatrick.jl
+++ b/test/runtests_layer_overlap_fitzpatrick.jl
@@ -70,11 +70,11 @@
         P = Psi_of(r)
         te = te_of(P)
         ne = ne_of(P)
-        p = slayer_parameters(; n_e=ne, t_e=te, t_i=te,
+        p = slayer_parameters(; P_perp_model=:chi_perp_e, n_e=ne, t_e=te, t_i=te,
             omega_e=1.0e4, omega_i=5.0e3,
             qval=q_of(r), sval_r=s_of(r), bt=B0,
             rs=r * A_MIN, R0=R0, mu_i=MNUM, zeff=ZEFF,
-            chi_perp=CHI, chi_tor=CHI, m=1, n=n,
+            chi_perp_e=CHI, chi_tor=CHI, m=1, n=n,
             resistivity_model=SpitzerHarmModel(), lnLambda_form=:nrl)
         return p
     end
@@ -259,6 +259,9 @@ end
         @test overlap_inputs(Dict("profile_file" => "own_zeff.h5", "impurity_charge" => 10.0, "zeff" => 2.5); dir=tmp).zeff(kin.psi[j]) ≈ 3.2
         Equilibrium.write_kinetic_h5(joinpath(tmp, "no_ni.h5"), Equilibrium.KineticProfileData(; base...))
         @test overlap_inputs(Dict("profile_file" => "no_ni.h5", "zeff" => 2.5); dir=tmp).zeff == 2.5
+        # A P_perp_model whose χ⊥,i the file lacks skips the scan with a warning rather than aborting the run.
+        @test (@test_logs (:warn, r"Layer-overlap scan skipped") match_mode = :any overlap_inputs(
+            Dict("profile_file" => "no_ni.h5", "P_perp_model" => "D_perp"); dir=tmp)) === nothing
     end
 
     widths(over) = GeneralizedPerturbedEquilibrium._layer_overlap_scan(overlap_inputs(over), _fresh(), equil).width_dr
@@ -283,7 +286,7 @@ end
         profs = Utilities.KineticProfiles(; psi=psi_grid,
             n_e=fill(1.0e19, npsi), T_e=fill(500.0, npsi), T_i=fill(500.0, npsi),
             omega_E=zeros(npsi))
-        scan = Tearing.resistive_layer_overlap(equil, profs; n_tor=1, m_max=0, extrapolate=false)
+        scan = Tearing.resistive_layer_overlap(equil, profs; P_perp_model=:chi_perp_e, n_tor=1, m_max=0, extrapolate=false)
         @test scan isa Tearing.LayerOverlapScan
         @test isempty(scan.psi)
         @test scan.psihigh === nothing && scan.first_overlap === nothing
diff --git a/test/runtests_slayer_inputs.jl b/test/runtests_slayer_inputs.jl
index bbd80bacc..070a57217 100644
--- a/test/runtests_slayer_inputs.jl
+++ b/test/runtests_slayer_inputs.jl
@@ -2,6 +2,7 @@
     using GeneralizedPerturbedEquilibrium
     using GeneralizedPerturbedEquilibrium.Equilibrium
     using GeneralizedPerturbedEquilibrium.Utilities
+    using GeneralizedPerturbedEquilibrium.Utilities: MU_0
     using GeneralizedPerturbedEquilibrium.InnerLayer
     using GeneralizedPerturbedEquilibrium.ForceFreeStates: SingType
     using TOML
@@ -75,7 +76,7 @@
         # dr_val=0.0 bypasses the build_slayer_inputs requirement that sing.restype be
         # pre-populated by ForceFreeStates.resist_eval_all! — the test sings here are
         # minimal stubs without restype, so we supply dr_val explicitly.
-        sl = build_slayer_inputs(equil, sings, profiles; bt=2.0, dr_val=0.0)
+        sl = build_slayer_inputs(equil, sings, profiles; P_perp_model=:chi_perp_e, bt=2.0, dr_val=0.0)
 
         @test length(sl) == 2
         @test sl[1] isa SLAYERParameters
@@ -128,9 +129,9 @@
         sings = [_mk_sing(; psi=0.3, q=2.0, q1=1.5, m=2, n=1)]
         bt_phys = Float64(equil.profiles.F_spline(0.3)) / (2π * equil.ro)
 
-        sl_default = build_slayer_inputs(equil, sings, profiles; dr_val=0.0)
-        sl_at_phys = build_slayer_inputs(equil, sings, profiles; bt=bt_phys, dr_val=0.0)
-        sl_double = build_slayer_inputs(equil, sings, profiles; bt=2 * bt_phys, dr_val=0.0)
+        sl_default = build_slayer_inputs(equil, sings, profiles; P_perp_model=:chi_perp_e, dr_val=0.0)
+        sl_at_phys = build_slayer_inputs(equil, sings, profiles; P_perp_model=:chi_perp_e, bt=bt_phys, dr_val=0.0)
+        sl_double = build_slayer_inputs(equil, sings, profiles; P_perp_model=:chi_perp_e, bt=2 * bt_phys, dr_val=0.0)
 
         # Leaving bt unset must be identical to passing the physical field explicitly.
         @test sl_default[1].lu ≈ sl_at_phys[1].lu rtol = 1e-12
@@ -141,35 +142,53 @@
         @test sl_double[1].lu != sl_default[1].lu
     end
 
-    @testset "build_slayer_inputs: chi_perp/chi_tor as scalars and callables" begin
+    @testset "build_slayer_inputs: chi_perp_e/chi_tor as scalars and callables" begin
         sings = [_mk_sing(; psi=0.5, q=2.4, q1=1.2, m=2, n=1)]
 
         # Scalar (dr_val=0.0 bypasses the sing.restype requirement; see comment above)
-        sl_s = build_slayer_inputs(equil, sings, profiles;
-            bt=2.0, chi_perp=2.0, chi_tor=1.5, dr_val=0.0)
+        sl_s = build_slayer_inputs(equil, sings, profiles; P_perp_model=:chi_perp_e,
+            bt=2.0, chi_perp_e=2.0, chi_tor=1.5, dr_val=0.0)
         # Callable with matching value
         chi_p(psi) = 2.0 + 0.0 * psi
         chi_t(psi) = 1.5 + 0.0 * psi
-        sl_c = build_slayer_inputs(equil, sings, profiles;
-            bt=2.0, chi_perp=chi_p, chi_tor=chi_t, dr_val=0.0)
+        sl_c = build_slayer_inputs(equil, sings, profiles; P_perp_model=:chi_perp_e,
+            bt=2.0, chi_perp_e=chi_p, chi_tor=chi_t, dr_val=0.0)
         @test sl_s[1].P_perp ≈ sl_c[1].P_perp
         @test sl_s[1].P_tor ≈ sl_c[1].P_tor
 
         # Callable with ψ-dependence changes the result
         chi_p_var(psi) = 1.0 + 10.0 * psi                     # χ⊥(0.5) = 6.0 > 2.0
-        sl_var = build_slayer_inputs(equil, sings, profiles;
-            bt=2.0, chi_perp=chi_p_var, chi_tor=1.5, dr_val=0.0)
+        sl_var = build_slayer_inputs(equil, sings, profiles; P_perp_model=:chi_perp_e,
+            bt=2.0, chi_perp_e=chi_p_var, chi_tor=1.5, dr_val=0.0)
         # P_perp = τ_r · χ⊥ / r² grows with χ⊥, so the varying-χ case at
         # ψ=0.5 (χ⊥=6) gives a *larger* P_perp than the scalar χ⊥=2.
         @test sl_var[1].P_perp > sl_s[1].P_perp
         @test sl_var[1].P_perp ≈ sl_s[1].P_perp * 6.0 / 2.0 rtol = 1e-10
     end
 
+    @testset "build_slayer_inputs: P_perp_model forwards the profile gradients and minor radius" begin
+        sings = [_mk_sing(psi=0.5, q=2.4, q1=1.2, m=2, n=1)]
+        kw = (; bt=2.0, chi_perp_e=2.0, chi_perp_i=1.0, chi_tor=1.5, dr_val=0.0)
+        ref = build_slayer_inputs(equil, sings, profiles; P_perp_model=:chi_perp_e, kw...)[1]
+        # Flat n_e and linear T_e, T_i: the D_⊥ weights are T_s′/(T_e′ + T_i′) = 7/13 and 6/13
+        dp = build_slayer_inputs(equil, sings, profiles; kw..., P_perp_model=:D_perp)[1]
+        D_hand = dp.c_beta^2 * dp.eta / MU_0 + (2 / 3) * (1 - dp.c_beta^2) * (2.0 * 7 / 13 + 1.0 * 6 / 13)
+        @test dp.P_perp ≈ ref.P_perp * D_hand / 2.0 rtol = 1e-8
+        # τ̄_E mode takes a from the boundary surface on the flux label, whatever rs_method is
+        te = build_slayer_inputs(equil, sings, profiles; kw..., P_perp_model=:tau_E, tau_E=0.1)[1]
+        a = first(radial_label(equil; rs_method=:flux))(1.0)
+        ref_fsa = build_slayer_inputs(equil, sings, profiles; P_perp_model=:chi_perp_e, kw..., rs_method=:fsa)[1]
+        te_fsa = build_slayer_inputs(equil, sings, profiles; kw..., P_perp_model=:tau_E, tau_E=0.1, rs_method=:fsa)[1]
+        @test te_fsa.P_perp ≈ ref_fsa.P_perp * ((3 / 16) * a^2 / 0.1) / 2.0 rtol = 1e-12
+        @test te.P_perp ≈ ref.P_perp * ((3 / 16) * a^2 / 0.1) / 2.0 rtol = 1e-12
+        @test te.P_tor == te.P_perp
+    end
+
     @testset "build_slayer_inputs: rs_method radial labels are self-consistent" begin
         sings = [_mk_sing(; psi=0.5, q=2.4, q1=1.2, m=2, n=1)]
         got = Dict{Symbol,Any}()
         for rsm in (:midplane, :halfwidth, :fsa, :volume, :flux)
-            sl = build_slayer_inputs(equil, sings, profiles; bt=2.0, dr_val=0.0, rs_method=rsm)
+            sl = build_slayer_inputs(equil, sings, profiles; P_perp_model=:chi_perp_e, bt=2.0, dr_val=0.0, rs_method=rsm)
             got[rsm] = sl[1]
             @test isfinite(sl[1].rs) && sl[1].rs > 0
             @test isfinite(sl[1].k_ref) && sl[1].k_ref > 0
@@ -188,24 +207,24 @@
         sings = [_mk_sing(; psi=0.5, q=2.4, q1=1.2, m=2, n=1)]
 
         # dc_type=:none and dr_val=0.0 → dc_tmp = 0 regardless of dr_val
-        sl_none = build_slayer_inputs(equil, sings, profiles;
+        sl_none = build_slayer_inputs(equil, sings, profiles; P_perp_model=:chi_perp_e,
             bt=2.0, dc_type=:none, dr_val=0.0)
         @test sl_none[1].dc_tmp == 0.0
 
         # dc_type=:fitzpatrick with dr_val = 0 still gives zero
-        sl_rf0 = build_slayer_inputs(equil, sings, profiles;
+        sl_rf0 = build_slayer_inputs(equil, sings, profiles; P_perp_model=:chi_perp_e,
             bt=2.0, dc_type=:fitzpatrick, dr_val=0.0)
         @test sl_rf0[1].dc_tmp == 0.0
 
         # dc_type=:fitzpatrick with dr_val > 0 → nonzero negative offset
-        sl_rf = build_slayer_inputs(equil, sings, profiles;
+        sl_rf = build_slayer_inputs(equil, sings, profiles; P_perp_model=:chi_perp_e,
             bt=2.0, dc_type=:fitzpatrick, dr_val=0.01)
         @test sl_rf[1].dc_tmp < 0
         @test isfinite(sl_rf[1].dc_tmp)
     end
 
     @testset "build_slayer_inputs: empty sings returns empty vector" begin
-        sl = build_slayer_inputs(equil, SingType[], profiles; bt=2.0)
+        sl = build_slayer_inputs(equil, SingType[], profiles; P_perp_model=:chi_perp_e, bt=2.0)
         @test sl isa Vector{SLAYERParameters}
         @test isempty(sl)
     end
@@ -217,8 +236,8 @@
         # where the factor is applied.
         s1 = [_mk_sing(; psi=0.3, q=2.0, q1=1.5, m=2, n=1)]
         s2 = [_mk_sing(; psi=0.3, q=2.0, q1=1.5, m=4, n=2)]
-        sl1 = build_slayer_inputs(equil, s1, profiles; bt=2.0, dr_val=0.0)
-        sl2 = build_slayer_inputs(equil, s2, profiles; bt=2.0, dr_val=0.0)
+        sl1 = build_slayer_inputs(equil, s1, profiles; P_perp_model=:chi_perp_e, bt=2.0, dr_val=0.0)
+        sl2 = build_slayer_inputs(equil, s2, profiles; P_perp_model=:chi_perp_e, bt=2.0, dr_val=0.0)
 
         @test sl2[1].Q_e / sl2[1].tauk ≈ 2 * (sl1[1].Q_e / sl1[1].tauk) rtol = 1e-12
         @test sl2[1].Q_i / sl2[1].tauk ≈ 2 * (sl1[1].Q_i / sl1[1].tauk) rtol = 1e-12
diff --git a/test/runtests_slayer_params.jl b/test/runtests_slayer_params.jl
index 4052dc2c5..f8f5a6081 100644
--- a/test/runtests_slayer_params.jl
+++ b/test/runtests_slayer_params.jl
@@ -12,7 +12,7 @@
             omega_e=-1.0e4, omega_i=5.0e3,
             qval=2.0, sval_r=1.0, bt=2.0,
             rs=0.5, R0=1.7, mu_i=2.0, zeff=1.0,
-            chi_perp=1.0, chi_tor=1.0,
+            chi_perp_e=1.0, chi_tor=1.0,
             m=2, n=1,
             dr_val=dr_val, dgeo_val=0.5, dc_type=dc_type,
             ising=3
@@ -21,7 +21,7 @@
 
     @testset "Test 1: round-trip from dimensional inputs" begin
         @info "Building SLAYERParameters from a reference deuterium case"
-        p = slayer_parameters(; _ref_kwargs()...)
+        p = slayer_parameters(; P_perp_model=:chi_perp_e, _ref_kwargs()...)
 
         # Identity / passthrough
         @test p.ising == 3
@@ -47,7 +47,7 @@
         # Default resistivity closure is neoclassical (Sauter F_33). The
         # trapped-particle correction raises η above plain Spitzer at the
         # same lnΛ: η_neo = η_Sp / F_33 with F_33 ∈ (0,1).
-        p_spitzer = slayer_parameters(; _ref_kwargs()...,
+        p_spitzer = slayer_parameters(; P_perp_model=:chi_perp_e, _ref_kwargs()...,
             resistivity_model=SpitzerModel(), lnLambda_form=:nrl)
         @test p.eta > p_spitzer.eta
         @test 1.1 < p.eta / p_spitzer.eta < 4.0    # banana-regime F_33 ≈ 0.3–0.9
@@ -57,7 +57,7 @@
         # formulas), so this is a genuine cross-check, not a tautology.
         lnLamb_wesson = 24.0 + 3.0 * log(10.0) - 0.5 * log(5.0e19) + log(1000.0)
         eta_wesson = 1.65e-9 * lnLamb_wesson / (1000.0 / 1e3)^1.5
-        p_legacy = slayer_parameters(; _ref_kwargs()...,
+        p_legacy = slayer_parameters(; P_perp_model=:chi_perp_e, _ref_kwargs()...,
             resistivity_model=SpitzerHarmModel(), lnLambda_form=:wesson)
         @test p_legacy.eta ≈ eta_wesson rtol = 3e-2
 
@@ -79,7 +79,7 @@
         # Compressibility is in (0,1) for finite β
         @test 0.0 < p.c_beta < 1.0
 
-        # Prandtl-like ratios are positive and equal here (chi_perp=chi_tor=1)
+        # Prandtl-like ratios are positive and equal here (chi_perp_e=chi_tor=1)
         @test p.P_perp ≈ p.P_tor
         @test p.P_perp > 0
 
@@ -96,12 +96,12 @@
         # All four dc_type branches must produce finite, non-NaN values
         # and respect the signs/structure of the formulas in
         # the SLAYER params dc_tmp formulas.
-        p_none = slayer_parameters(; _ref_kwargs(; dr_val=0.01, dc_type=:none)...)
+        p_none = slayer_parameters(; P_perp_model=:chi_perp_e, _ref_kwargs(; dr_val=0.01, dc_type=:none)...)
         @test p_none.dc_tmp == 0.0   # :none ignores dr_val
 
-        p_lar = slayer_parameters(; _ref_kwargs(; dr_val=0.01, dc_type=:lar)...)
-        p_rf = slayer_parameters(; _ref_kwargs(; dr_val=0.01, dc_type=:fitzpatrick)...)
-        p_tor = slayer_parameters(; _ref_kwargs(; dr_val=0.01, dc_type=:toroidal)...)
+        p_lar = slayer_parameters(; P_perp_model=:chi_perp_e, _ref_kwargs(; dr_val=0.01, dc_type=:lar)...)
+        p_rf = slayer_parameters(; P_perp_model=:chi_perp_e, _ref_kwargs(; dr_val=0.01, dc_type=:fitzpatrick)...)
+        p_tor = slayer_parameters(; P_perp_model=:chi_perp_e, _ref_kwargs(; dr_val=0.01, dc_type=:toroidal)...)
 
         @test isfinite(p_lar.dc_tmp)
         @test isfinite(p_rf.dc_tmp)
@@ -113,12 +113,12 @@
         @test p_tor.dc_tmp < 0
 
         # Sign flips with sign of dr_val
-        p_lar_neg = slayer_parameters(;
+        p_lar_neg = slayer_parameters(; P_perp_model=:chi_perp_e,
             _ref_kwargs(; dr_val=-0.01, dc_type=:lar)...)
         @test sign(p_lar_neg.dc_tmp) == -sign(p_lar.dc_tmp)
 
         # Reject unknown dc_type
-        @test_throws ArgumentError slayer_parameters(;
+        @test_throws ArgumentError slayer_parameters(; P_perp_model=:chi_perp_e,
             _ref_kwargs(; dr_val=0.01, dc_type=:bogus)...)
     end
 
@@ -143,9 +143,9 @@
     @testset "Test 1d: reject diamagnetic frequencies that leave D undefined" begin
         # iota_e = ω_*e/(ω_*e − ω_*i) must be finite and positive for D ∝ √iota_e.
         _drifts(ωe, ωi) = merge(_ref_kwargs(), (omega_e=ωe, omega_i=ωi))
-        @test_throws ArgumentError slayer_parameters(; _drifts(-1.0e4, -1.0e4)...)   # degenerate: iota_e singular
-        @test_throws ArgumentError slayer_parameters(; _drifts(-1.0e4, -2.0e4)...)   # same sign, |ω_*i| > |ω_*e|: iota_e = -1
-        @test_throws ArgumentError slayer_parameters(; _drifts(0.0, 5.0e3)...)       # no electron drift: iota_e = 0
+        @test_throws ArgumentError slayer_parameters(; P_perp_model=:chi_perp_e, _drifts(-1.0e4, -1.0e4)...)   # degenerate: iota_e singular
+        @test_throws ArgumentError slayer_parameters(; P_perp_model=:chi_perp_e, _drifts(-1.0e4, -2.0e4)...)   # same sign, |ω_*i| > |ω_*e|: iota_e = -1
+        @test_throws ArgumentError slayer_parameters(; P_perp_model=:chi_perp_e, _drifts(0.0, 5.0e3)...)       # no electron drift: iota_e = 0
     end
 
     @testset "Test 2: r-based shear conversion" begin
@@ -181,8 +181,8 @@
         # diagnostic keeping the sign. dc_type=:lar with nonzero dr_val exercises the Wd
         # iteration and the critical-Δ square roots as well as tau_h.
         base = _ref_kwargs(; dr_val=-0.1, dc_type=:lar)
-        pos = slayer_parameters(; base...)
-        neg = slayer_parameters(; merge(base, (; sval_r=-1.0))...)
+        pos = slayer_parameters(; P_perp_model=:chi_perp_e, base...)
+        neg = slayer_parameters(; P_perp_model=:chi_perp_e, merge(base, (; sval_r=-1.0))...)
 
         # The sign survives where it is a diagnostic, not a magnitude.
         @test pos.sval_r == 1.0
@@ -199,4 +199,45 @@
         @test neg.lu > 0
         @test neg.tauk > 0
     end
+
+    @testset "Test 4: P_perp_model options" begin
+        base = merge(_ref_kwargs(), (chi_tor=1.7,))
+        ref = slayer_parameters(; P_perp_model=:chi_perp_e, base...)
+        @test ref.P_perp ≈ ref.tau_r * 1.0 / ref.rs^2 rtol = 1e-14
+        @test slayer_parameters(; base..., P_perp_model=:chi_perp_i, chi_perp_i=2.0).P_perp ≈ 2 * ref.P_perp
+        @test slayer_parameters(; base..., P_perp_model=:P_phi).P_perp == ref.P_tor
+        # Eq. 16: with flat temperatures the χ terms vanish and only c_β²η/μ₀ survives, so P_perp = c_β²
+        grads = (dlnn=-1.0, dlnTe=-2.0, dlnTi=-3.0)
+        d0 = slayer_parameters(; base..., P_perp_model=:D_perp, dlnn=-1.0, dlnTe=0.0, dlnTi=0.0, chi_perp_i=2.0)
+        @test d0.P_perp ≈ d0.c_beta^2 rtol = 1e-12
+        # η_e = 2, η_i = 3, T_e = T_i ⇒ τ = 3/4, τ_e = 3/7, τ_i = 4/7
+        d1 = slayer_parameters(; base..., P_perp_model=:D_perp, grads..., chi_perp_i=2.0)
+        D_hand = d1.c_beta^2 * d1.eta / MU_0 + (2 / 3) * (1 - d1.c_beta^2) * ((2 / 3) * (3 / 7) * 1.0 + (3 / 4) * (4 / 7) * 2.0)
+        @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
+        # A flat ion profile is finite (no 0·∞); a vanishing pressure gradient or a missing gradient is a clear error
+        @test_throws "pressure gradient vanishes" slayer_parameters(; base..., P_perp_model=:D_perp, dlnn=0.0, dlnTe=0.0, dlnTi=0.0, chi_perp_i=2.0)
+        @test_throws "needs finite dlnn" slayer_parameters(; base..., P_perp_model=:D_perp, dlnn=-1.0, chi_perp_i=2.0)
+        @test_throws "positive, finite chi_perp_e" slayer_parameters(; base..., P_perp_model=:chi_perp_e, chi_perp_e=-1.0)
+        @test isfinite(slayer_parameters(; base..., P_perp_model=:D_perp, dlnn=0.0, dlnTe=-2.0, dlnTi=0.0, chi_perp_i=2.0).P_perp)
+        # A non-positive or missing χ⊥,i is an error in every model that uses it
+        @test_throws ArgumentError slayer_parameters(; base..., P_perp_model=:chi_perp_i)
+        @test_throws ArgumentError slayer_parameters(; base..., P_perp_model=:chi_perp_i, chi_perp_i=0.0)
+        @test_throws ArgumentError slayer_parameters(; base..., P_perp_model=:D_perp)
+        @test slayer_parameters(; base..., P_perp_model=:c_beta).P_perp == ref.c_beta^2
+        @test_throws ArgumentError slayer_parameters(; base..., P_perp_model=:bogus)
+        # τ̄_E = 0.1 s and a = 0.6 m ⇒ uniform χ_E = (3/16)·0.36/0.1 m²/s sets both Prandtl numbers
+        te = slayer_parameters(; base..., P_perp_model=:tau_E, tau_E=0.1, a_minor=0.6)
+        @test te.P_tor ≈ te.tau_r * ((3 / 16) * 0.6^2 / 0.1) / te.rs^2 rtol = 1e-14
+        @test te.P_perp == te.P_tor
+        @test_throws ArgumentError slayer_parameters(; base..., P_perp_model=:tau_E)
+        @test_throws ArgumentError slayer_parameters(; base..., P_perp_model=:tau_E, tau_E=0.1)
+        # Only P_perp moves; the critical-Δ offset keeps using chi_perp_e.
+        cb = slayer_parameters(; merge(base, (dr_val=-1.0, dc_type=:lar))..., P_perp_model=:c_beta)
+        @test cb.dc_tmp == slayer_parameters(; P_perp_model=:chi_perp_e, merge(base, (dr_val=-1.0, dc_type=:lar))...).dc_tmp
+        @test cb.P_tor == ref.P_tor
+    end
 end
diff --git a/test/runtests_slayer_profile_transport.jl b/test/runtests_slayer_profile_transport.jl
new file mode 100644
index 000000000..0ca137890
--- /dev/null
+++ b/test/runtests_slayer_profile_transport.jl
@@ -0,0 +1,32 @@
+using Test
+using GeneralizedPerturbedEquilibrium.Equilibrium: KineticProfileData, write_kinetic_h5
+using GeneralizedPerturbedEquilibrium.Runner
+
+@testset "SLAYER loads only transport profiles used by its model" begin
+    mktempdir() do dir
+        psi = collect(range(0.0, 1.0; length=5))
+        function load_transport(model; chi_e=ones(5), chi_i=ones(5), chi_phi=ones(5), dc_type=:none)
+            data = KineticProfileData(; psi, n_e=fill(5e19, 5), T_e=1000 .- 500 .* psi,
+                T_i=800 .- 300 .* psi, chi_e, chi_i, chi_phi)
+            write_kinetic_h5(joinpath(dir, "profiles.h5"), data)
+            ctrl = SLAYERControl(; enabled=true, profile_file="profiles.h5", P_perp_model=model, tau_E=0.1, D_c_type=dc_type)
+            return Runner._load_profiles(ctrl, dir)
+        end
+
+        @test load_transport(:tau_E; chi_phi=(-ones(5))).chi_tor === nothing
+        for model in (:chi_perp_e, :P_phi, :c_beta, :tau_E)
+            @test load_transport(model; chi_i=(-ones(5))).chi_perp_i === nothing
+        end
+        for model in (:chi_perp_i, :P_phi, :c_beta, :tau_E)
+            @test load_transport(model; chi_e=(-ones(5))).chi_perp_e === nothing
+            @test_throws "chi_e" load_transport(model; chi_e=(-ones(5)), dc_type=:lar)
+        end
+        @test_throws "chi_i" load_transport(:D_perp; chi_i=(-ones(5)))
+        @test_throws "chi_i" load_transport(:chi_perp_i; chi_i=(-ones(5)))
+        @test_throws "chi_e" load_transport(:D_perp; chi_e=(-ones(5)))
+        @test_throws "chi_e" load_transport(:chi_perp_e; chi_e=(-ones(5)))
+        @test_throws "chi_phi" load_transport(:chi_perp_e; chi_phi=(-ones(5)))
+        @test load_transport(:D_perp; chi_i=zeros(5)).chi_perp_i === nothing
+        @test load_transport(:D_perp).chi_perp_i(0.5) ≈ 1.0
+    end
+end
diff --git a/test/runtests_slayer_riccati.jl b/test/runtests_slayer_riccati.jl
index 48a6263b2..b6ab240de 100644
--- a/test/runtests_slayer_riccati.jl
+++ b/test/runtests_slayer_riccati.jl
@@ -13,12 +13,12 @@
     # ratio is ~2.4 (vs ~0.5 at 1 keV), placing the fixture solidly on the large_D side of the
     # branch boundary. All other inputs unchanged.
     function _ref_params_large_D()
-        return slayer_parameters(;
+        return slayer_parameters(; P_perp_model=:chi_perp_e,
             n_e=5.0e19, t_e=3000.0, t_i=3000.0,
             omega_e=-1.0e4, omega_i=5.0e3,
             qval=2.0, sval_r=1.0, bt=2.0,
             rs=0.5, R0=1.7, mu_i=2.0, zeff=1.0,
-            chi_perp=1.0, chi_tor=1.0,
+            chi_perp_e=1.0, chi_tor=1.0,
             m=2, n=1)
     end
 
diff --git a/test/runtests_slayer_runner.jl b/test/runtests_slayer_runner.jl
index 2c4a42d5a..2d9afd6e9 100644
--- a/test/runtests_slayer_runner.jl
+++ b/test/runtests_slayer_runner.jl
@@ -58,6 +58,7 @@
             "bt" => 1.8,
             "ion_mass" => 2.0,
             "D_R" => 0.01,
+            "P_perp_model" => "D_perp",
             "scan_grid" => Dict{String,Any}(
                 "Q_re_range" => [-5.0, 5.0],
                 "Q_im_range" => [-1.0, 3.0],
@@ -80,6 +81,7 @@
         @test c.msing_max == 2
         @test c.bt === 1.8
         @test c.D_R == 0.01
+        @test c.P_perp_model === :D_perp
         @test c.Q_re_range == (-5.0, 5.0)
         @test c.Q_im_range == (-1.0, 3.0)
         @test c.nre == 50
@@ -93,6 +95,22 @@
         # Unknown keys should raise
         bad = merge(section, Dict{String,Any}("mistyped_key" => 42))
         @test_throws ArgumentError slayer_control_from_toml(bad)
+        # The renamed χ⊥,e key and bad P_perp_model values are rejected
+        @test_throws "chi_perp` is now `chi_perp_e" slayer_control_from_toml(merge(section, Dict{String,Any}("chi_perp" => 1.0)))
+        @test_throws ArgumentError Runner.validate(slayer_control_from_toml(merge(section, Dict{String,Any}("P_perp_model" => "bogus"))))
+        @test_throws ArgumentError Runner.validate(slayer_control_from_toml(merge(section, Dict{String,Any}("P_perp_model" => "tau_E"))))
+        @test_throws ArgumentError Runner.validate(slayer_control_from_toml(merge(section, Dict{String,Any}("P_perp_model" => "tau_E", "tau_E" => Inf))))
+        c_te = slayer_control_from_toml(merge(section, Dict{String,Any}("P_perp_model" => "tau_E", "tau_E" => 0.05)))
+        @test c_te.tau_E == 0.05 && Runner.validate(c_te) === c_te
+        # D_perp is the default, and it refuses to run without an ion diffusivity from the file or the deck.
+        @test SLAYERControl().P_perp_model === :D_perp
+        no_chi = (; chi_perp_e=nothing, chi_perp_i=nothing, chi_tor=nothing)
+        @test_throws "needs the ion diffusivity" Runner._slayer_chi_inputs(no_chi, SLAYERControl(); warn=false)
+        @test Runner._slayer_chi_inputs(no_chi, SLAYERControl(; chi_perp_i=0.7); warn=false).chi_perp_i == 0.7
+        @test Runner._slayer_chi_inputs(no_chi, SLAYERControl(; P_perp_model=:chi_perp_e); warn=false).chi_perp_i === nothing
+        # Scalar diffusivities must be positive.
+        @test_throws "chi_perp_i must be positive" Runner.validate(SLAYERControl(; chi_perp_i=-0.5))
+        @test_throws "chi_perp_e must be positive" Runner.validate(SLAYERControl(; chi_perp_e=0.0))
     end
 
     @testset "TearingProblem solve forwards surfaces and Δ'" begin