diff --git a/examples/DIIID-like_SLAYER_example/gpec.toml b/examples/DIIID-like_SLAYER_example/gpec.toml
index a8a99f4c7..667f83943 100644
--- a/examples/DIIID-like_SLAYER_example/gpec.toml
+++ b/examples/DIIID-like_SLAYER_example/gpec.toml
@@ -73,15 +73,17 @@ coupling_mode           = "uncoupled" # "uncoupled" (per-surface) or "coupled" (
 dc_type                 = "none"     # Critical-Δ offset selector: "none", "lar", "rfitzp", or "toroidal"
 mu_i                    = 2.0        # Ion mass in proton-mass units (2.0 = deuterium)
 zeff                    = 1.0        # Effective charge
-chi_perp                = 1.0       # fallback only; the kinetic file supplies χ⊥(ψ)
+chi_perp_e              = 1.0       # fallback only; the kinetic file supplies χ⊥,e(ψ)
+chi_perp_i              = 1.0       # fallback only; the kinetic file supplies χ⊥,i(ψ)
 chi_tor                 = 1.0       # fallback only; the kinetic file supplies χ_φ(ψ)
+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])
 pole_threshold_adaptive = true       # SLAYER |Δ| spans many orders; adapt to 10·median(|Δ|)
 filter_above_poles      = true       # Discard roots above the highest pole γ
 filter_outside_re       = true       # Condition the above-pole filter on the +γ step exiting the Re(Δ)=0 contour loop
 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 2f084f080..c1a305586 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/src/Equilibrium/KineticProfiles.jl b/src/Equilibrium/KineticProfiles.jl
index 934256554..18037d7e7 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 χ⊥                                                                                 |
+| `chi_e`             | electron perpendicular heat diffusivity χ⊥,e                                                                      |
+| `chi_i`             | ion perpendicular heat diffusivity χ⊥,i                                                                           |
 | `chi_phi`           | toroidal momentum diffusivity χ_φ                                                                                 |
 | `species_densities` | named per-species densities (e.g. `"n_D"`, `"n_T"`) for explicit multi-ion input; `nothing` for ASCII / no extras |
 | `provenance`        | short string recording the source file/format                                                                     |
@@ -80,17 +81,18 @@ 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}}
     species_densities::Union{Nothing,Dict{String,Vector{Float64}}}
     provenance::String
 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,
+    omega_E=nothing, omega_tor=nothing, chi_e=nothing, chi_i=nothing, chi_phi=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), species_densities, String(provenance))
+        _v(omega_E), _v(omega_tor), _v(chi_e), _v(chi_i), _v(chi_phi), species_densities, String(provenance))
 end
 
 const _KINETIC_H5_EXTS = (".h5", ".hdf5", ".he5")
@@ -98,7 +100,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")
+    "omega_tor" => "rad/s", "chi_e" => "m^2/s", "chi_i" => "m^2/s", "chi_phi" => "m^2/s")
 
 """
     read_kinetic_file(path; group="/") -> KineticProfileData
@@ -137,7 +139,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")
+        standard = ("psi", "n_i", "n_e", "T_i", "T_e", "omega_E", "omega_tor", "chi_e", "chi_i", "chi_phi")
         extras = Dict{String,Vector{Float64}}()
         for k in keys(g)
             (k in standard || !startswith(k, "n_")) && continue
@@ -147,7 +149,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"),
+            chi_e=rd("chi_e"), chi_i=rd("chi_i"), chi_phi=rd("chi_phi"),
             species_densities=(isempty(extras) ? nothing : extras), provenance=prov)
     end
 end
@@ -178,6 +180,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)
         # Named per-species densities (multi-ion input); dataset names must be `n_*` to round-trip.
         if data.species_densities !== nothing
diff --git a/src/InnerLayer/SLAYER/LayerInputs.jl b/src/InnerLayer/SLAYER/LayerInputs.jl
index 252e81ae9..a8530573f 100644
--- a/src/InnerLayer/SLAYER/LayerInputs.jl
+++ b/src/InnerLayer/SLAYER/LayerInputs.jl
@@ -181,8 +181,14 @@ profiles, without an intermediate file round-trip.
 
   - `mu_i`      -- ion mass in proton-mass units (default `2.0` for D).
   - `zeff`      -- effective charge (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` (default `1.0`); used by `P_perp_model = :chi_perp_i` and `:D_perp`.
+  - `tau_E` -- whole-plasma energy confinement time [s], one scalar for every surface;
+    required by `P_perp_model=:tau_E`.
+  - `P_perp_model` -- how `P_perp` is built; see [`slayer_parameters`](@ref)
+    (default `:chi_perp_e`).
   - `chi_tor`   -- toroidal heat diffusivity [m²/s]. Scalar or a callable
     of `psi` (default `1.0`).
   - `dr_val`    -- resistive interchange index `D_R = E + F + H²`
@@ -232,7 +238,10 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles;
     mu_i::Real=2.0,
     zeff::Real=1.0,
     z_i::Real=1.0,
-    chi_perp=1.0,
+    chi_perp_e=1.0,
+    chi_perp_i=1.0,
+    tau_E::Real=NaN,
+    P_perp_model::Symbol=:chi_perp_e,
     chi_tor=1.0,
     dr_val=nothing,
     dgeo_val=nothing,
@@ -262,13 +271,11 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles;
     # 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)
@@ -291,6 +298,7 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles;
         n_res = sing.n[1]
 
         prof = profiles(psi)
+        kin = _kinetic_at(psi)
         # Take ω_*e, ω_*i from the spline derivatives, or from `profiles` when the caller
         # supplies them directly. `run_slayer` supplies zeros, so the latter is a library path.
         ω_e_use, ω_i_use = compute_omega_star ? _omega_star_at(psi, n_res) : (prof.omega_e, prof.omega_i)
@@ -378,8 +386,12 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles;
             omega=prof.omega, 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=zeff,
-            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,
+            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 ae875dea4..73b9e8c22 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, :rfitzp, :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
 
@@ -132,7 +135,7 @@ end
 # Internal: solve the Wd self-consistency loop for the chi_parallel-based
 # critical Δ (Connor-Hastie-Helander 2015). 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 ||
@@ -149,7 +152,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
@@ -165,13 +168,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 === :rfitzp
         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
@@ -179,8 +182,11 @@ end
 """
     slayer_parameters(; n_e, t_e, t_i, omega, 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=1.0, tau_E=NaN,
+                        dlnn=NaN, dlnTe=NaN, dlnTi=NaN,
+                        P_perp_model=:chi_perp_e,
                         dr_val=0.0, dgeo_val=0.0,
                         dc_type=:none, ising=0,
                         resistivity_model=SauterNeoModel(),
@@ -210,7 +216,23 @@ 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_e`, `chi_tor` -- electron perpendicular heat / toroidal momentum diffusivity [m²/s]
+  - `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`)
+  - `P_perp_model` -- how `P_perp` is built:
+      + `:chi_perp_e` (default) -- τ_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)
+      + `:D_perp` -- τ_R/τ_⊥ with τ_⊥ = r_s²/D_⊥ and D_⊥ from Eq. 16 of Fitzpatrick, Phys. Plasmas 29,
+        032507 (2022), with η_⊥ = η_∥:
+        D_⊥ = c_β² η/μ₀ + (2/3)(1 − c_β²)[η_e τ_e/(1+η_e) χ⊥,e + η_i τ_i/(1+η_i) χ⊥,i]
+      + `:c_beta` -- P_perp = C² with C = c_β (J.-K. Park, Phys. Plasmas 29, 072506 (2022), Eq. 10, with the
+        thermal-conduction term K dropped). This is the χ⊥ → 0 limit of `:D_perp` (P_⊥ = μ₀D_⊥/η = c_β²), so it is
+        a floor on P_perp; realistic anomalous transport makes D_⊥ much larger.
+      + `:tau_E`  -- P_perp = P_tor = τ_R/τ_E, overriding χ_φ (Fitzpatrick, Phys. Plasmas 30, 092512 (2023), Eq. 131)
+    `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`, `:rfitzp`, `:toroidal`
@@ -263,8 +285,11 @@ function slayer_parameters(;
     omega::Real, 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::Real=1.0, tau_E::Real=NaN,
+    dlnn::Real=NaN, dlnTe::Real=NaN, dlnTi::Real=NaN,
+    P_perp_model::Symbol=:chi_perp_e,
     dr_val::Real=0.0, dgeo_val::Real=0.0,
     dc_type::Symbol=:none, ising::Integer=0,
     resistivity_model::NeoResistivityModel=SauterNeoModel(),
@@ -366,10 +391,41 @@ function slayer_parameters(;
     c_beta = sqrt(lbeta / (1.0 + lbeta))
 
     # Effective Prandtl-like transport ratios
-    tau_perp = rs^2 / chi_perp
-    P_perp = tau_r / tau_perp
     tau_tor = rs^2 / chi_tor
     P_tor = tau_r / tau_tor
+    if P_perp_model === :tau_E
+        isfinite(tau_E) && tau_E > 0 ||
+            throw(ArgumentError("slayer_parameters: P_perp_model=:tau_E needs a positive tau_E, got $tau_E"))
+        # P_⊥ = P_φ = τ_R/τ_E (Fitzpatrick, Phys. Plasmas 30, 092512 (2023), Eq. 131)
+        P_tor = tau_r / tau_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_⊥ = P_φ (Fitzpatrick, Phys. Plasmas 30, 092512 (2023), Eq. 143)
+        P_tor
+    elseif P_perp_model === :D_perp
+        # 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)
+        isfinite(D_perp) && D_perp > 0 ||
+            throw(ArgumentError("slayer_parameters: P_perp_model=:D_perp gave D_perp=$D_perp; " *
+                                "check dlnn/dlnTe/dlnTi (got $dlnn, $dlnTe, $dlnTi)"))
+        # τ_⊥ = r_s²/D_⊥, P_⊥ = τ_R/τ_⊥
+        tau_r / (rs^2 / D_perp)
+    elseif P_perp_model === :c_beta
+        # C² = P_⊥ with C = c_β (J.-K. Park, Phys. Plasmas 29, 072506 (2022), Eq. 10, K = 0)
+        c_beta^2
+    else
+        throw(ArgumentError("slayer_parameters: unknown P_perp_model=$P_perp_model. " *
+                            "Allowed: $(ALLOWED_P_PERP_MODELS)"))
+    end
 
     # Normalized beta-related width and Δ-normalization
     d_beta = c_beta * d_i
@@ -378,7 +434,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 0b6db9761..431c3c806 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 87781259b..9f630838b 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 033f71eca..9b4d0d9a9 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/Runner/Control.jl b/src/Tearing/Runner/Control.jl
index 1a12ff7e6..98943f1b1 100644
--- a/src/Tearing/Runner/Control.jl
+++ b/src/Tearing/Runner/Control.jl
@@ -34,10 +34,16 @@ constructor.
     callable of `psi` overrides it
   - `mu_i`     -- ion mass in proton-mass units (default 2.0 for D)
   - `zeff`     -- effective charge
-  - `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
+  - `tau_E`    -- whole-plasma energy confinement time [s], required by
+    `P_perp_model = :tau_E` (default `nothing`)
+  - `P_perp_model` -- how the perpendicular Prandtl number is built: `:chi_perp_e`
+    (default), `:chi_perp_i`, `:P_phi` (P_perp = P_tor), `:D_perp` (D⊥ from χ⊥,e, χ⊥,i and the n, T gradients), `:c_beta`
+    (P_perp = c_β²), or `:tau_E` (P_perp = P_tor = τ_R/τ_E); see
+    `InnerLayer.SLAYER.slayer_parameters`
   - `dr_val`, `dgeo_val`  -- critical-Δ formula inputs. `nothing` (default)
     auto-derives them from the equilibrium: `dr_val` from the resistive
     interchange index `D_R = E + F + H²` at each surface, `dgeo_val` from the
@@ -84,9 +90,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 `"/"`)
 
@@ -107,8 +113,11 @@ there is one consistent interface for resistive and kinetic profiles.
     bt::Union{Float64,Nothing} = nothing
     mu_i::Float64 = 2.0
     zeff::Float64 = 1.0
-    chi_perp::Float64 = 1.0
+    chi_perp_e::Float64 = 1.0
     chi_tor::Float64 = 1.0
+    chi_perp_i::Float64 = 1.0
+    tau_E::Union{Float64,Nothing} = nothing
+    P_perp_model::Symbol = :chi_perp_e
     dr_val::Union{Float64,Nothing} = nothing
     dgeo_val::Union{Float64,Nothing} = nothing
     theta_sample::Float64 = 0.0
@@ -161,6 +170,7 @@ const _VALID_INNER_MODELS = (:slayer_fitzpatrick, :ggj_shooting, :ggj_galerkin)
 const _VALID_SCAN_MODES = (:amr, :brute_force)
 const _VALID_COUPLING_MODES = (:uncoupled, :coupled)
 const _VALID_DC_TYPES = (:none, :lar, :rfitzp, :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)
 
@@ -177,6 +187,11 @@ function validate(ctrl::SLAYERControl)
     ctrl.dc_type in _VALID_DC_TYPES ||
         throw(ArgumentError("SLAYERControl: dc_type=$(ctrl.dc_type) " *
                             "not in $(_VALID_DC_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)"))
+    ctrl.P_perp_model !== :tau_E || (ctrl.tau_E !== nothing && ctrl.tau_E > 0) ||
+        throw(ArgumentError("SLAYERControl: P_perp_model=tau_E requires a positive tau_E [s]"))
     ctrl.resistivity_model in _VALID_RESISTIVITY_MODELS ||
         throw(ArgumentError("SLAYERControl: resistivity_model=$(ctrl.resistivity_model) " *
                             "not in $(_VALID_RESISTIVITY_MODELS)"))
@@ -242,11 +257,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, :dc_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, :dr_val, :dgeo_val)
+        elseif sym in (:bt, :dr_val, :dgeo_val, :tau_E)
             # Allow explicit nothing (auto-derive) or a number (override)
             kwargs[sym] = v === nothing ? nothing : Float64(v)
         elseif sym === :boxes
diff --git a/src/Tearing/Runner/run_slayer.jl b/src/Tearing/Runner/run_slayer.jl
index 2459087cf..37cc6d1e0 100644
--- a/src/Tearing/Runner/run_slayer.jl
+++ b/src/Tearing/Runner/run_slayer.jl
@@ -18,9 +18,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`).
+# `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`).
 function _load_profiles(control::SLAYERControl, dir_path::AbstractString)
     isempty(control.profile_file) &&
         error("run_slayer: [SLAYER] profile_file is empty — point it at a " *
@@ -46,14 +46,15 @@ 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_perp_e = _chi_spline(data.chi_e)
+    chi_perp_i = _chi_spline(data.chi_i)
     chi_tor = _chi_spline(data.chi_phi)
-    return (profiles=profiles, chi_perp=chi_perp, chi_tor=chi_tor)
+    return (profiles=profiles, chi_perp_e=chi_perp_e, chi_perp_i=chi_perp_i, chi_tor=chi_tor)
 end
 
 # ---------------------------------------------------------------------
@@ -407,7 +408,7 @@ Orchestrate the full SLAYER analysis against a `ForceFreeStates.ForceFreeStatesR
 reading its equilibrium, singular surfaces and Δ' matrix. Kinetic profiles are
 read from `control.profile_file` (relative to `dir_path`) through the shared
 `Equilibrium.read_kinetic_file` reader; when the file carries `chi_e`/`chi_phi`
-profiles they set χ⊥(ψ)/χ_φ(ψ), otherwise the scalar `control.chi_perp`/
+profiles they set χ⊥(ψ)/χ_φ(ψ), otherwise the scalar `control.chi_perp_e`/
 `chi_tor` fallbacks are used.
 
 The toroidal field comes from `control.bt`; leaving it unset (the default) makes
@@ -454,18 +455,25 @@ function run_slayer(equil, surfaces::AbstractVector, delta_prime_matrix::Abstrac
         # 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_perp_e = loaded.chi_perp_e === nothing ? control.chi_perp_e : loaded.chi_perp_e
         chi_tor = loaded.chi_tor === nothing ? control.chi_tor : loaded.chi_tor
-        (loaded.chi_perp === nothing || loaded.chi_tor === nothing) && @warn(
+        (loaded.chi_perp_e === nothing || loaded.chi_tor === nothing) && @warn(
             "SLAYER: kinetic file has no usable chi_e/chi_phi profile(s) " *
             "(dataset absent or all-zero); using the scalar " *
-            "control.chi_perp/chi_tor fallback for the missing one(s).")
+            "control.chi_perp_e/chi_tor fallback for the missing one(s).")
+        chi_perp_i = loaded.chi_perp_i === nothing ? control.chi_perp_i : loaded.chi_perp_i
+        (control.P_perp_model in (:chi_perp_i, :D_perp) && loaded.chi_perp_i === nothing) && @warn(
+            "SLAYER: kinetic file has no usable chi_i profile; using the scalar " *
+            "control.chi_perp_i fallback.")
         params = build_slayer_inputs(equil, surfaces, profiles;
             bt=bt,
             mu_i=control.mu_i,
             zeff=control.zeff,
-            chi_perp=chi_perp,
+            chi_perp_e=chi_perp_e,
             chi_tor=chi_tor,
+            chi_perp_i=chi_perp_i,
+            tau_E=something(control.tau_E, NaN),
+            P_perp_model=control.P_perp_model,
             dr_val=control.dr_val,
             dgeo_val=control.dgeo_val,
             dc_type=control.dc_type,
diff --git a/test/runtests_dispersion_coupled.jl b/test/runtests_dispersion_coupled.jl
index 1b950016e..30efddbd4 100644
--- a/test/runtests_dispersion_coupled.jl
+++ b/test/runtests_dispersion_coupled.jl
@@ -24,7 +24,7 @@
             omega=0.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
diff --git a/test/runtests_dispersion_residual.jl b/test/runtests_dispersion_residual.jl
index e8af06da4..2d68531a2 100644
--- a/test/runtests_dispersion_residual.jl
+++ b/test/runtests_dispersion_residual.jl
@@ -24,7 +24,7 @@
             omega=0.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 scale defaults" begin
diff --git a/test/runtests_kinetic.jl b/test/runtests_kinetic.jl
index d502385f7..a8a8a8399 100644
--- a/test/runtests_kinetic.jl
+++ b/test/runtests_kinetic.jl
@@ -605,8 +605,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
@@ -634,6 +635,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_slayer_inputs.jl b/test/runtests_slayer_inputs.jl
index fce001036..26b51932c 100644
--- a/test/runtests_slayer_inputs.jl
+++ b/test/runtests_slayer_inputs.jl
@@ -142,24 +142,24 @@
         @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)
+            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)
+            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)
+            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
diff --git a/test/runtests_slayer_params.jl b/test/runtests_slayer_params.jl
index 1f665fb1d..f3f8cd5df 100644
--- a/test/runtests_slayer_params.jl
+++ b/test/runtests_slayer_params.jl
@@ -12,7 +12,7 @@
             omega=0.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,
             dr_val=dr_val, dgeo_val=0.5, dc_type=dc_type,
             ising=3
@@ -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
 
@@ -199,4 +199,31 @@
         @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(; 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 χ⊥,e = χ⊥,i = 0 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, grads..., chi_perp_e=0.0, chi_perp_i=0.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
+        @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)
+        te = slayer_parameters(; base..., P_perp_model=:tau_E, tau_E=0.1)
+        @test te.P_tor == te.tau_r / 0.1
+        @test te.P_perp == te.P_tor
+        @test_throws ArgumentError slayer_parameters(; base..., P_perp_model=:tau_E)
+        # 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(; 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_riccati.jl b/test/runtests_slayer_riccati.jl
index 9d4a95462..0755d51fd 100644
--- a/test/runtests_slayer_riccati.jl
+++ b/test/runtests_slayer_riccati.jl
@@ -18,7 +18,7 @@
             omega=0.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 ee114e37f..f94e9dc8f 100644
--- a/test/runtests_slayer_runner.jl
+++ b/test/runtests_slayer_runner.jl
@@ -57,6 +57,7 @@
             "bt" => 1.8,
             "mu_i" => 2.0,
             "dr_val" => 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],
@@ -79,6 +80,7 @@
         @test c.msing_max == 2
         @test c.bt === 1.8
         @test c.dr_val == 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
@@ -92,6 +94,12 @@
         # 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 ArgumentError 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"))))
+        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
     end
 
     @testset "run_slayer: result-facing form forwards surfaces and Δ'" begin