diff --git a/docs/development/references.md b/docs/development/references.md
index e19b92e32..43f3daf74 100644
--- a/docs/development/references.md
+++ b/docs/development/references.md
@@ -68,6 +68,19 @@ The PerturbedEquilibrium module implements GPEC-style perturbed equilibrium calc
   - Published: Physics of Plasmas **24**, 032505 (2017)
   - Describes: Self-consistent coupling with neoclassical effects
 
+## InnerLayer Module (SLAYER)
+
+The SLAYER slab inner-layer solver and the resistive-layer width diagnostics are based on:
+
+- **Burgess et al. (2026)**: "Tearing Stability Prediction Combining Toroidal Calculations With a Two-Fluid Slab Layer"
+  - 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 (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)
+  - Describes: Diffusive-resistive layer width (Eq. 100) and the resistive-layer-overlap truncation `0 < Ψ < 1 − ε_c` (Sect. 5.9)
+
 ## Resistive MHD Stability Analysis (Future Work)
 
 GPEC will eventually implement resistive MHD stability analysis based on:
diff --git a/docs/resources/2025-Fitzpatrick-Response of a magnetically diverted tokamak plasma to a resonant magnetic perturbation.pdf b/docs/resources/2025-Fitzpatrick-Response of a magnetically diverted tokamak plasma to a resonant magnetic perturbation.pdf
new file mode 100644
index 000000000..b12990f32
Binary files /dev/null and b/docs/resources/2025-Fitzpatrick-Response of a magnetically diverted tokamak plasma to a resonant magnetic perturbation.pdf differ
diff --git a/src/Equilibrium/GridRefinement.jl b/src/Equilibrium/GridRefinement.jl
index e58eb81dc..a600f12d5 100644
--- a/src/Equilibrium/GridRefinement.jl
+++ b/src/Equilibrium/GridRefinement.jl
@@ -57,6 +57,85 @@ const CORE_MODEL_PSI_MAX = 0.03
 const EDGE_MODEL_PSI_MIN = 0.9
 # θ-lines subsample stride for the 2D geometry channels
 const THETA_STRIDE = 8
+# --- shared separatrix edge q-law -------------------------------------------------------------
+# Minimum knots in the edge band before a fit is attempted.
+const EDGE_FIT_MIN_KNOTS = 4
+# Weak absolute floor on the diverging fit. The discrimination is the relative
+# r2_log > r2_linear comparison below (no absolute threshold separates diverted from limited
+# edges); this floor only rejects fits that describe nothing at all.
+const EDGE_FIT_MIN_R2 = 0.90
+
+# Least-squares slope and coefficient of determination for y = a + b*x.
+function _linfit_r2(x::Vector{Float64}, y::Vector{Float64})
+    x_bar = sum(x) / length(x)
+    y_bar = sum(y) / length(y)
+    sxx = sum((x .- x_bar) .^ 2)
+    sxx > 0 || return (NaN, NaN, NaN, NaN)
+    b = sum((x .- x_bar) .* (y .- y_bar)) / sxx
+    ss_res = sum((y .- (y_bar .+ b .* (x .- x_bar))) .^ 2)
+    ss_tot = sum((y .- y_bar) .^ 2)
+    r2 = ss_tot > 0 ? 1 - ss_res / ss_tot : NaN
+    return (b, r2, x_bar, y_bar)
+end
+
+"""
+    edge_q_law(equil; psi_max, psi_min=EDGE_MODEL_PSI_MIN, min_knots=EDGE_FIT_MIN_KNOTS,
+               min_r2=EDGE_FIT_MIN_R2) -> nothing | (; A, q_bar, u_bar, n_knots, r2_log, r2_linear)
+
+Least-squares fit of the separatrix edge law `q = q̄ + A·(ln(1−ψ) − ū)` over the equilibrium's
+outer knots — the single shared statement of that model, used both by the grid-refinement edge
+density floor and by the resistive-layer overlap scan's out-of-grid surface search.
+
+Returns `nothing` when the diverging model does **not** describe this equilibrium's edge, so a
+plasma with finite edge q is never extrapolated as if q blew up:
+
+  - fewer than `min_knots` knots in the band;
+  - `A ≥ 0`, i.e. q not rising toward ψ = 1;
+  - `r2_log < min_r2` — the log law does not actually fit;
+  - `r2_log ≤ r2_linear` — a plain linear-in-ψ fit explains the edge q at least as well, which is
+    what a **limited** plasma looks like. This comparison carries no scale and is what separates
+    the shipped limited decks (Solovev 0.760 vs 0.9996 linear; LAR 0.904 vs 0.998) from the
+    diverted ones (DIII-D 0.996 vs 0.865, 0.999 vs 0.594).
+
+This is a test of the **model**, not a topology classification: it asks whether q diverges
+logarithmically here, not whether an x-point exists. Geometric x-point detection is
+[`classify_topology`](@ref), which is a separate concern.
+"""
+function edge_q_law(equil::PlasmaEquilibrium;
+    psi_max::Real=Float64(equil.profiles.xs[end]),
+    psi_min::Real=EDGE_MODEL_PSI_MIN,
+    min_knots::Int=EDGE_FIT_MIN_KNOTS,
+    min_r2::Real=EDGE_FIT_MIN_R2)
+    xs = collect(Float64, equil.profiles.xs)
+    band = findall(x -> x >= psi_min && x < psi_max, xs)
+    if length(band) < min_knots
+        n_tail = max(min_knots, length(xs) ÷ 10)
+        band = filter(i -> xs[i] < psi_max, collect(max(1, length(xs) - n_tail + 1):length(xs)))
+    end
+    length(band) >= min_knots || return nothing
+
+    q = [Float64(equil.profiles.q_spline(xs[i])) for i in band]
+    u = [log(1.0 - xs[i]) for i in band]
+    all(isfinite, u) && all(isfinite, q) || return nothing
+
+    A, r2_log, u_bar, q_bar = _linfit_r2(u, q)
+    _, r2_linear, _, _ = _linfit_r2([xs[i] for i in band], q)
+    (isfinite(A) && A < 0) || return nothing            # q must rise toward the edge
+    (isfinite(r2_log) && r2_log >= min_r2) || return nothing
+    (isfinite(r2_linear) && r2_log > r2_linear) || return nothing
+
+    return (A=A, q_bar=q_bar, u_bar=u_bar, n_knots=length(band), r2_log=r2_log, r2_linear=r2_linear)
+end
+
+"""
+ψ at which the edge law reaches `q_target`; closed form, no root-finding needed.
+"""
+edge_q_law_psi(fit, q_target::Real) = 1.0 - exp(fit.u_bar + (q_target - fit.q_bar) / fit.A)
+
+"""
+dq/dψ from the edge law: q = q̄ + A·ln(1−ψ) + const ⇒ dq/dψ = −A/(1−ψ).
+"""
+edge_q_law_dqdpsi(fit, psi::Real) = -fit.A / (1.0 - psi)
 # Rational-surface bracketing (Δ′ robustness). The ideal-MHD Δ′ asymptotic matching samples the
 # cubic equilibrium splines' 2nd/3rd derivatives across each rational ψ_s over the matching stencil
 # [ψ_s − dpsi, ψ_s + dpsi], dpsi = singfac_min/|n·q′|. A cubic 3rd derivative is piecewise constant
@@ -248,10 +327,15 @@ function _knot_density(equil::PlasmaEquilibrium; tau::Float64, kin::Union{Nothin
     # nodal data of the smallest flux surfaces is dominated by integration and axis
     # extrapolation error, so measured curvature is not trusted below the core split.
     dlog = (4.0 * tau)^(1 / 3)
+    # The edge floor encodes the DIVERGING edge law q ≈ -A·ln(1-ψ), so it is applied only where
+    # that model actually describes the equilibrium. A limited plasma has finite edge q and must
+    # not be packed as if q blew up. The density itself stays A-independent (that is the point of
+    # the form: uniform relative q′ error regardless of A) -- the fit supplies validity, not slope.
+    edge_diverges = edge_q_law(equil) !== nothing
     @inbounds for i in 1:n
         if xs[i] <= CORE_MODEL_PSI_MAX
             rho_s[i] = 1.0 / (dlog * xs[i])
-        elseif xs[i] >= EDGE_MODEL_PSI_MIN
+        elseif edge_diverges && xs[i] >= EDGE_MODEL_PSI_MIN
             rho_s[i] = max(rho_s[i], 1.0 / (dlog * (1.0 - xs[i])))
         end
         rho_s[i] = max(rho_s[i], 1.0 / H_TARGET_MAX)
@@ -409,9 +493,33 @@ function bracket_mandatory_nodes(grid::Vector{Float64}, centers::Vector{Float64}
     return merged
 end
 
+"""
+    _truncate_density(xs, rho, psihigh) -> (xs_t, rho_t)
+
+Restrict a measured knot density to `[xs[1], psihigh]`, linearly interpolating `rho` at the new
+outer endpoint so the density integral stays continuous in `psihigh`. Errors when `psihigh` lies
+outside the sampled grid, where the density is unmeasured.
+"""
+function _truncate_density(xs, rho::Vector{Float64}, psihigh::Float64)
+    xs_v = collect(Float64, xs)
+    psihigh > xs_v[1] ||
+        error("_truncate_density: psihigh=$psihigh must exceed the inner grid bound $(xs_v[1])")
+    psihigh <= xs_v[end] + 1e-12 ||
+        error(
+            "_truncate_density: psihigh=$psihigh exceeds the pass-1 grid end $(xs_v[end]); " *
+            "the knot density is unmeasured there — form the enlarged domain first, then refine against it"
+        )
+    psihigh >= xs_v[end] - 1e-12 && return (xs_v, rho)
+
+    k = searchsortedlast(xs_v, psihigh)
+    frac = (psihigh - xs_v[k]) / (xs_v[k+1] - xs_v[k])
+    rho_end = rho[k] + frac * (rho[k+1] - rho[k])
+    return (vcat(xs_v[1:k], psihigh), vcat(rho[1:k], rho_end))
+end
+
 """
     refined_psi_grid(equil::PlasmaEquilibrium; tau, kin=nothing, mandatory=Float64[],
-                     singfac_min=1e-4, n_min=1, bracket_coef=BRACKET_COEF,
+                     psihigh=nothing, singfac_min=1e-4, n_min=1, bracket_coef=BRACKET_COEF,
                      min_spacing=MIN_KNOT_SPACING, N_cap=1024) -> Vector{Float64}
 
 Build the refined pass-2 ψ grid from a formed pass-1 equilibrium: measured-curvature knot
@@ -423,11 +531,18 @@ whose pedestal gradients attract knots; `mandatory` lists rational-surface ψ va
 `dpsi = singfac_min/(n_min·|q′|)`, and the bracket half-width is `bracket_coef·dpsi` (floored at
 `min_spacing`). Rational surfaces are bracketed, not pinned: a knot on the surface would make the
 Δ′ extraction's cubic 3rd derivative jump mid-stencil (see `BRACKET_COEF`).
+
+`psihigh` builds the grid for a domain *smaller* than the one `equil` was formed on: the measured
+density is truncated there and the last node lands exactly on it. This lets a pass-1 equilibrium
+supply the density for a re-form on a reduced domain without an extra solve. Passing a `psihigh`
+beyond the pass-1 grid is an error — the density out there is unmeasured, so an enlarged domain
+must be formed first and refined against that.
 """
 function refined_psi_grid(equil::PlasmaEquilibrium;
     tau::Float64,
     kin::Union{Nothing,KineticProfileSplines}=nothing,
     mandatory::Vector{Float64}=Float64[],
+    psihigh::Union{Nothing,Real}=nothing,
     singfac_min::Float64=1e-4,
     n_min::Int=1,
     bracket_coef::Float64=BRACKET_COEF,
@@ -435,6 +550,9 @@ function refined_psi_grid(equil::PlasmaEquilibrium;
     N_cap::Int=REFINED_N_CAP)
     xs = equil.profiles.xs
     rho = _knot_density(equil; tau, kin)
+    if psihigh !== nothing
+        xs, rho = _truncate_density(xs, rho, Float64(psihigh))
+    end
     # Floor the density to a fixed (τ-independent) locally-uniform fine patch around each rational
     # so Δ′ has the resolution to sample 3rd derivatives there at any accuracy target (see
     # RATIONAL_RES_SPACING).
diff --git a/src/ForceFreeStates/CoreTypes.jl b/src/ForceFreeStates/CoreTypes.jl
index 98d4807e1..cbce9f841 100644
--- a/src/ForceFreeStates/CoreTypes.jl
+++ b/src/ForceFreeStates/CoreTypes.jl
@@ -135,6 +135,7 @@ gpec.toml.
   - `numunorms_init::Int` - Initial array size for solution normalization data
   - `singfac_min::Float64` - Fractional distance from rational q at which ideal jump condition is enforced
   - `set_psilim_via_dmlim::Bool` - Truncate the integration domain at `(last_rational_q + dmlim) / n` rather than at `qhigh` / `psihigh`. Fortran STRIDE found that truncating ~20 % above the outermost rational (`dmlim = 0.2`) avoids a numerical kink instability in δW that appears when the integration ends too close to or just below a rational surface. **For diverted equilibria where q → ∞ at the separatrix** (e.g. DIII-D geqdsks, the bulk of production use) this costs negligible physical domain because rationals get arbitrarily dense near the LCFS — `set_psilim_via_dmlim = true` is the safe and recommended default. **For limited circular / analytical equilibria with finite q at the edge** (Solovev, LAR scans), rationals are sparse and 20 % above the last rational chops off too much edge, so set `set_psilim_via_dmlim = false` and let `qhigh` / `psihigh` control the truncation. Multi-`n` runs are not supported by this truncation (the "outermost rational + dmlim / n" depends on which `n`); when `set_psilim_via_dmlim = true` with `nn_low != nn_high`, `sing_lim!` warns and falls back to `qhigh` / `psihigh`. Default `true`.
+  - `psilim_from_layer_overlap::Bool` - Cap the integration domain at the point where adjacent rational surfaces' resistive layers overlap, after Fitzpatrick, Nucl. Fusion 2025 Sect. 5.9. Past that point no surface retains a well-separated inner region, so the matched-asymptotic treatment is not defined there. Applied as an upper bound on `qlim` at the top of `sing_lim!`, so `dmlim` / `qhigh` still select the final surface from inside it, and a bound beyond `psihigh` is inert. The scan needs kinetic profiles and is run whenever they are readable — its result is recorded under `ForceFreeStates/LayerOverlap/` regardless of this flag, so a run always shows whether layer physics would have constrained the domain. Single-`n` only, like `dmlim`. Default `false`.
   - `dmlim::Float64` - Distance beyond last rational surface (normalised ∈ [0,1) in units of 1/n). Only used when `set_psilim_via_dmlim` is true. Fortran STRIDE convention is 0.2 (truncate 20 % of one rational-surface spacing above the last surface), retained here.
   - `sing_order::Int` - Order of singular layer (Frobenius) expansion at rational surfaces. Default 6 (Fortran STRIDE convention for Δ' calculations; lower values trade accuracy for speed).
   - `qhigh::Float64` - Integration terminated at q limit determined by minimum of qhigh and qa from equil
@@ -172,6 +173,7 @@ gpec.toml.
     numunorms_init::Int = 100
     singfac_min::Float64 = 1e-4   # Matches Fortran STRIDE; required nonzero for the Riccati path.
     set_psilim_via_dmlim::Bool = true   # Safe default for diverted equilibria (most production use); set false for limited/analytical (LAR, Solovev). Auto-skipped for multi-n. See docstring.
+    psilim_from_layer_overlap::Bool = false  # Cap psilim where adjacent resistive layers overlap. Opt-in; the scan is recorded either way. See docstring.
     dmlim::Float64 = 0.2
     sing_order::Int = 6
     qhigh::Float64 = 1e3
diff --git a/src/ForceFreeStates/Surfaces/Finding.jl b/src/ForceFreeStates/Surfaces/Finding.jl
index 2c7bdac14..cbe7c805a 100644
--- a/src/ForceFreeStates/Surfaces/Finding.jl
+++ b/src/ForceFreeStates/Surfaces/Finding.jl
@@ -110,7 +110,8 @@ performed to find the corresponding `psilim` to integrate to.
 Note that the Newton iteration will be triggered if either `set_psilim_via_dmlim` is true
 or `ctrl.qhigh < equil.params.qmax`. Otherwise, the equilibrium edge values are used.
 """
-function sing_lim!(intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium)
+function sing_lim!(intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium;
+    psilim_cap::Union{Nothing,Real}=nothing)
 
     profiles = equil.profiles
 
@@ -119,6 +120,20 @@ function sing_lim!(intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl,
     intr.q1lim = profiles.q_deriv(profiles.xs[end]; hint=Ref(profiles.npts_minus_1))
     intr.psilim = equil.params.psihigh_resolved
 
+    # Resistive-layer overlap imposes an UPPER BOUND on the domain: past the first pair of
+    # overlapping layers no surface retains a well-separated inner region, so matched asymptotics
+    # is not defined out there. Applied as a cap on qlim rather than as psilim directly, so the
+    # dmlim / qhigh truncation below still selects the final surface from inside the bound.
+    # Never widens the domain: a cap beyond psihigh is inert by construction.
+    if psilim_cap !== nothing && psilim_cap < intr.psilim
+        q_cap = profiles.q_spline(Float64(psilim_cap))
+        if q_cap < intr.qlim
+            @info "Resistive-layer overlap caps the domain: qlim $(@sprintf("%.3f", intr.qlim)) -> " *
+                  "$(@sprintf("%.3f", q_cap)) (psi $(@sprintf("%.6f", intr.psilim)) -> $(@sprintf("%.6f", Float64(psilim_cap))))"
+            intr.qlim = q_cap
+        end
+    end
+
     # Optionally override qlim based on dmlim (Fortran sas_flag=t equivalent). The cutoff reads
     # the *resolved* toroidal range on `intr`, so callers must assign intr.nlow / intr.nhigh
     # before calling; an unresolved range is an error rather than a silent change of truncation
diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl
index f2ed4ef65..6bca04ea0 100755
--- a/src/GeneralizedPerturbedEquilibrium.jl
+++ b/src/GeneralizedPerturbedEquilibrium.jl
@@ -241,7 +241,9 @@ function main_from_inputs(
     ffs_start = time()
 
     locstab, ballooning_boundary = run_local_stability(ctrl, equil)
-    metric, mats = prepare_force_free_states!(intr, ctrl, equil, kf_ctrl, kinetic_profiles; species=kf_species)
+    metric, mats, layer_overlap = prepare_force_free_states!(intr, ctrl, equil, kf_ctrl, kinetic_profiles;
+        species=kf_species,
+        overlap_profile_file=_overlap_profile_path(inputs, kf_ctrl, intr.dir_path))
     ffs_result = run_force_free_states(ctrl, equil, mats, intr, metric)
 
     if ctrl.write_outputs_to_HDF5
@@ -251,7 +253,8 @@ function main_from_inputs(
             inputs=inputs,
             forcing_modes=forcing_modes_snapshot,
             locstab=locstab,
-            ballooning_boundary=ballooning_boundary
+            ballooning_boundary=ballooning_boundary,
+            layer_overlap=layer_overlap
         )
         @info "Results written to $(ctrl.HDF5_filename)"
     end
@@ -479,8 +482,50 @@ function run_local_stability(ctrl::ForceFreeStatesControl, equil::Equilibrium.Pl
     return locstab, ballooning_boundary
 end
 
+# Kinetic profiles for the overlap scan can come from either the NTV path
+# (`[KineticForces] kinetic_file`) or the tearing path (`[SLAYER] profile_file`); the shipped
+# DIII-D SLAYER deck carries only the latter. Prefer whichever exists on disk.
+function _overlap_profile_path(inputs, kf_ctrl::KineticForces.KineticForcesControl, dir_path::AbstractString)
+    _resolve(f) = isempty(f) ? nothing : (isabspath(f) ? String(f) : joinpath(dir_path, f))
+    for cand in (_resolve(kf_ctrl.kinetic_file),
+        (inputs isa AbstractDict && haskey(inputs, "SLAYER")) ?
+        _resolve(get(inputs["SLAYER"], "profile_file", "")) : nothing)
+        cand !== nothing && isfile(cand) && return cand
+    end
+    return nothing
+end
+
+# Run the resistive-layer overlap scan for the run's toroidal mode number, or return `nothing`
+# when it cannot be run (no kinetic file, multi-n, or the scan itself refuses). Reads the kinetic
+# file directly rather than reusing `kinetic_profiles`, which is a `KineticProfileSplines` built
+# for the NTV path and carries a different field set than the layer builders take.
+function _layer_overlap_scan(path::Union{Nothing,AbstractString},
+    intr::ForceFreeStatesInternal, equil::Equilibrium.PlasmaEquilibrium)
+    (path === nothing || !isfile(path)) && return nothing
+    if intr.nlow != intr.nhigh
+        @info "Layer-overlap scan skipped: the overlap point depends on n, and this is a multi-n run (nn_low=$(intr.nlow), nn_high=$(intr.nhigh))."
+        return nothing
+    end
+    try
+        data = Equilibrium.read_kinetic_file(path)
+        (data.n_e === nothing || data.T_e === nothing || data.T_i === nothing) && return nothing
+        npsi = length(data.psi)
+        omega = data.omega_E === nothing ? zeros(npsi) : collect(Float64, data.omega_E)
+        profiles = Utilities.KineticProfiles(; psi=collect(Float64, data.psi),
+            n_e=collect(Float64, data.n_e), T_e=collect(Float64, data.T_e),
+            T_i=collect(Float64, data.T_i), omega=omega,
+            omega_e=zeros(npsi), omega_i=zeros(npsi))
+        # Eq. (100) is not covariant -- it is anchored to the toroidal-flux label of the paper's
+        # Eq. (30), so the scan is driven in that label regardless of the SLAYER default.
+        return Tearing.resistive_layer_overlap(equil, profiles; n_tor=intr.nlow, rs_method=:flux)
+    catch err
+        @warn "Layer-overlap scan failed; the integration domain is untouched." exception = (err, catch_backtrace())
+        return nothing
+    end
+end
+
 """
-    prepare_force_free_states!(intr, ctrl, equil, kf_ctrl, kinetic_profiles) -> (metric, mats)
+    prepare_force_free_states!(intr, ctrl, equil, kf_ctrl, kinetic_profiles) -> (metric, mats, overlap)
 
 Set up the force-free-states solve on `intr`: integration limits, the surviving singular
 surfaces and their GGJ coefficients, the poloidal mode range, and the metric plus
@@ -492,10 +537,21 @@ function prepare_force_free_states!(
     equil::Equilibrium.PlasmaEquilibrium,
     kf_ctrl::KineticForces.KineticForcesControl,
     kinetic_profiles;
-    species=nothing
+    species=nothing,
+    overlap_profile_file::Union{Nothing,AbstractString}=nothing
 )
+    # Resistive-layer overlap: locate where adjacent rational surfaces' layers run into each
+    # other, and use it as an upper bound on the integration domain. The scan runs whenever
+    # kinetic profiles are readable so `ForceFreeStates/LayerOverlap/` always records the point,
+    # but it only constrains the domain when the user opts in.
+    overlap = _layer_overlap_scan(overlap_profile_file, intr, equil)
+    psilim_cap = (ctrl.psilim_from_layer_overlap && overlap !== nothing) ? overlap.psihigh : nothing
+    if ctrl.psilim_from_layer_overlap && overlap === nothing
+        @warn "psilim_from_layer_overlap = true but no layer-overlap scan was available; the domain is untouched."
+    end
+
     # Determine psilim and qlim (where we will integrate to)
-    sing_lim!(intr, ctrl, equil)
+    sing_lim!(intr, ctrl, equil; psilim_cap=psilim_cap)
 
     # Find all singular surfaces in the equilibrium
     sing_find!(intr, equil)
@@ -594,7 +650,7 @@ function prepare_force_free_states!(
         end
     end
 
-    return metric, mats
+    return metric, mats, overlap
 end
 
 """
@@ -776,7 +832,7 @@ function solve(prob::EulerLagrangeProblem, alg::ForceFreeStates.AbstractIntegrat
     end
 
     locstab, ballooning_boundary = run_local_stability(ctrl, equil)
-    metric, mats = prepare_force_free_states!(intr, ctrl, equil, kf_ctrl, nothing)
+    metric, mats, _ = prepare_force_free_states!(intr, ctrl, equil, kf_ctrl, nothing)
     result = run_force_free_states(ctrl, equil, mats, intr, metric)
 
     if ctrl.write_outputs_to_HDF5
@@ -1051,7 +1107,8 @@ function write_outputs_to_HDF5(
     inputs::Union{Nothing,Dict{String,Any}}=nothing,
     forcing_modes::Union{Nothing,Vector{ForcingTerms.ForcingMode}}=nothing,
     locstab::Union{FastInterpolations.CubicSeriesInterpolant,Nothing}=nothing,
-    ballooning_boundary=(psi=Float64[], alpha=Float64[], alpha_critical=Float64[])
+    ballooning_boundary=(psi=Float64[], alpha=Float64[], alpha_critical=Float64[]),
+    layer_overlap=nothing
 )
 
     ctrl = result.control
@@ -1173,6 +1230,30 @@ function write_outputs_to_HDF5(
         out_h5["LocalStability/alpha"] = ballooning_boundary.alpha
         out_h5["LocalStability/alpha_critical"] = ballooning_boundary.alpha_critical
 
+        # Resistive-layer overlap scan. Written whenever the scan ran, whether or not it
+        # constrained the domain, so a run always shows where layer physics would have cut.
+        if layer_overlap !== nothing
+            lo = layer_overlap
+            out_h5["ForceFreeStates/LayerOverlap/m"] = lo.m
+            out_h5["ForceFreeStates/LayerOverlap/n"] = lo.n
+            out_h5["ForceFreeStates/LayerOverlap/psi"] = lo.psi
+            out_h5["ForceFreeStates/LayerOverlap/r_s"] = lo.rs
+            out_h5["ForceFreeStates/LayerOverlap/delta_s_abs"] = lo.delta_s_m
+            out_h5["ForceFreeStates/LayerOverlap/width_delta_s"] = lo.width_delta_s
+            out_h5["ForceFreeStates/LayerOverlap/width_visco"] = lo.width_visco
+            out_h5["ForceFreeStates/LayerOverlap/width_dr"] = lo.width_dr
+            out_h5["ForceFreeStates/LayerOverlap/extrapolated"] = Int.(lo.extrapolated)
+            out_h5["ForceFreeStates/LayerOverlap/psilim_overlap"] =
+                lo.psihigh === nothing ? NaN : lo.psihigh
+            out_h5["ForceFreeStates/LayerOverlap/first_overlap_index"] =
+                lo.first_overlap === nothing ? -1 : lo.first_overlap
+            # "applied" records that the bound was active going into sing_lim!, not that it set
+            # the final psilim -- dmlim/qhigh may truncate deeper inside it.
+            out_h5["ForceFreeStates/LayerOverlap/applied"] =
+                Int(ctrl.psilim_from_layer_overlap && lo.psihigh !== nothing &&
+                    lo.psihigh < equil.params.psihigh_resolved)
+        end
+
         # Write integration data: the ψ trace and integrator diagnostics from the raw ODE state,
         # the ξ profiles from the solution. Either may be absent (Galerkin has no ODE state;
         # Riccati has no ξ solution), in which case the datasets are written empty.
diff --git a/src/HDF5Schema.jl b/src/HDF5Schema.jl
index 274f54f14..b9d9a62b8 100644
--- a/src/HDF5Schema.jl
+++ b/src/HDF5Schema.jl
@@ -134,6 +134,33 @@ const MAIN_H5_ANNOTATIONS = [
         (; long_name="Glasser-Greene-Johnson resistive interchange criterion D_R", dims=("psi",), attach=(1 => "LocalStability/psi",)),
     "LocalStability/ballooning_Delta_prime" =>
         (; long_name="high-n ballooning Δ' (distinct from the tearing Δ')", dims=("psi",), attach=(1 => "LocalStability/psi",)),
+    # Resistive-layer overlap scan (Fitzpatrick, Nucl. Fusion 2025, Sect. 5.9). Its own surface
+    # axis: the list includes surfaces extrapolated beyond the integration domain, so it does not
+    # align with the SingularSurfaces/ axis.
+    "ForceFreeStates/LayerOverlap/psi" =>
+        (; long_name="rational-surface location scanned for resistive-layer overlap", scale="surface_overlap"),
+    "ForceFreeStates/LayerOverlap/m" =>
+        (; long_name="poloidal mode number of each scanned surface", dims=("surface_overlap",), attach=(1 => "ForceFreeStates/LayerOverlap/psi",)),
+    "ForceFreeStates/LayerOverlap/n" =>
+        (; long_name="toroidal mode number of each scanned surface", dims=("surface_overlap",), attach=(1 => "ForceFreeStates/LayerOverlap/psi",)),
+    "ForceFreeStates/LayerOverlap/r_s" =>
+        (; long_name="minor radius of each scanned surface in the Fitzpatrick flux label", units="m", dims=("surface_overlap",), attach=(1 => "ForceFreeStates/LayerOverlap/psi",)),
+    "ForceFreeStates/LayerOverlap/delta_s_abs" =>
+        (; long_name="Riccati resistive layer thickness |δ_s|", units="m", dims=("surface_overlap",), attach=(1 => "ForceFreeStates/LayerOverlap/psi",)),
+    "ForceFreeStates/LayerOverlap/width_delta_s" =>
+        (; long_name="|δ_s| expressed as a normalized-flux width", dims=("surface_overlap",), attach=(1 => "ForceFreeStates/LayerOverlap/psi",)),
+    "ForceFreeStates/LayerOverlap/width_visco" =>
+        (; long_name="visco-resistive comparison scale as a normalized-flux width", dims=("surface_overlap",), attach=(1 => "ForceFreeStates/LayerOverlap/psi",)),
+    "ForceFreeStates/LayerOverlap/width_dr" =>
+        (; long_name="diffusive-resistive layer width (Fitzpatrick 2025 Eq. 100) as a normalized-flux width; the channel the criterion uses", dims=("surface_overlap",), attach=(1 => "ForceFreeStates/LayerOverlap/psi",)),
+    "ForceFreeStates/LayerOverlap/extrapolated" =>
+        (; long_name="1 where the surface was located beyond the equilibrium grid via the separatrix edge q-law", dims=("surface_overlap",), attach=(1 => "ForceFreeStates/LayerOverlap/psi",)),
+    "ForceFreeStates/LayerOverlap/psilim_overlap" =>
+        (; long_name="domain limit implied by the first layer overlap, at the last surface with a well-separated inner region (NaN when no overlap was found in range)"),
+    "ForceFreeStates/LayerOverlap/first_overlap_index" =>
+        (; long_name="index into this group's surface axis of the first surface overlapping its inner neighbour (-1 when none do)"),
+    "ForceFreeStates/LayerOverlap/applied" =>
+        (; long_name="1 when the overlap bound was handed to sing_lim! as an active cap on qlim (dmlim/qhigh may still truncate deeper); 0 when recorded only"),
     "LocalStability/ballooning_psi" =>
         (; long_name="normalized poloidal flux ψ_N of the ballooning α boundary scan", scale="psi_ballooning"),
     "LocalStability/alpha" =>
diff --git a/src/InnerLayer/SLAYER/LayerThickness.jl b/src/InnerLayer/SLAYER/LayerThickness.jl
index 0b6db9761..9a70a2d38 100644
--- a/src/InnerLayer/SLAYER/LayerThickness.jl
+++ b/src/InnerLayer/SLAYER/LayerThickness.jl
@@ -137,6 +137,26 @@ is built from, retained as a drift-scale reference.
   - `delta_s`   -- complex layer thickness `δ_s = dels_db · d_β` [m]
   - `delta_s_m` -- `|δ_s|`, the resistive layer thickness in meters (primary)
   - `d_beta`    -- β-weighted ion scale `c_β·d_i` in meters (drift reference)
+  - `delta_norm`  -- layer normalization length `r_s · S^(-1/3)` in meters, the
+    length by which Δ' is made dimensionless (`Δ̂' = Δ'·delta_norm`). Shared by all
+    four regimes of Burgess et al. (2026) Eqs. (10)-(13), not specific to any one
+    of them, and **not** the classic non-rotating Furth-Killeen-Rosenbluth (1963)
+    constant-ψ tearing width, which scales as `S^(-2/5)` and is not computed here.
+    Exactly the reciprocal of the `delta_n` Δ-normalization already carried on
+    `SLAYERParameters`, restated as a length for comparison with `delta_s_m`
+  - `delta_dr` -- diffusive-resistive layer thickness in meters, Fitzpatrick (2025)
+    Eq. (100). Strong magnetic shear near the separatrix forces every resonant layer in
+    that region into the diffusive-resistive (`nu = 1/4`) regime, so this is the width the
+    edge overlap criterion compares against surface spacing. In the paper's normalized
+    radius, `delta_hat = tau_A^(1/2) / (tau_R^(1/4) tau_E^(1/4) d_beta_hat^(1/2) (n|s|)^(1/2))`;
+    with `lu = tau_R/tau_A` and `P_perp = tau_R/tau_E` that is
+    `lu^(-1/2) P_perp^(1/4) / (d_beta_hat^(1/2) (n|s|)^(1/2))`, scaled by `rs` for meters.
+    Distinct from `delta_visco`, which is the 2/3 power of the same timescale grouping and
+    carries neither the `d_beta` nor the shear factor
+  - `delta_visco` -- viscous-resistive scale `delta_norm · P_perp^(1/6)` in meters,
+    from the `P^(1/6)` broadening of the VR-regime growth rate, Burgess et al.
+    (2026) Eq. (11). The paper gives the growth rate rather than an explicit width;
+    this is the corresponding length
 
 `delta_s_m` should sit within a few orders of magnitude of `d_beta` for a
 well-posed surface (`dels_db` is O(1)); a large gap flags a normalisation
@@ -150,6 +170,9 @@ struct LayerWidths
     delta_s::ComplexF64
     delta_s_m::Float64
     d_beta::Float64
+    delta_norm::Float64
+    delta_visco::Float64
+    delta_dr::Float64
 end
 
 """
@@ -161,11 +184,27 @@ surface.
 Runs [`riccati_del_s`](@ref) for the dimensionless `δ_s / d_β` and scales
 by `p.d_beta` to obtain `δ_s` in meters. Keyword arguments are forwarded
 to `riccati_del_s`.
+
+The two algebraic comparison scales `delta_norm` and `delta_visco` come from
+`p` directly and do not depend on the Riccati solve.
 """
 function slayer_layer_thickness(p::SLAYERParameters; kwargs...)
     dels_db = riccati_del_s(p; kwargs...)
     delta_s = dels_db * p.d_beta
+    # Layer normalization length: exactly 1/delta_n (delta_n = S^(1/3)/r_s), computed
+    # from rs and lu so it reads as a length. Burgess et al. (2026), the Δ̂' = Δ'/S^(1/3)
+    # normalization preceding Eq. (10).
+    delta_norm = p.rs * p.lu^(-1.0 / 3.0)
+    # Viscous-resistive broadening, P^(1/6) coefficient of Burgess et al. (2026) Eq. (11).
+    delta_visco = delta_norm * p.P_perp^(1.0 / 6.0)
+    # Diffusive-resistive width, Fitzpatrick (2025) Eq. (100), in meters. d_beta_hat is the
+    # normalized ion scale d_beta/rs; the shear is the r-based |s| SLAYER already carries.
+    # The paper's tau_A (its Eq. 74) carries no shear, whereas SLAYER's tau_h divides by n*s,
+    # so lu = tau_R/tau_h = (n|s|) * (tau_R/tau_A). Substituting that into Eq. (100) cancels its
+    # explicit (n|s|)^(-1/2) exactly, leaving no shear dependence in terms of lu.
+    d_beta_hat = p.d_beta / p.rs
+    delta_dr = d_beta_hat > 0 ? p.rs * p.lu^(-0.5) * p.P_perp^(0.25) / sqrt(d_beta_hat) : NaN
     return LayerWidths(p.ising, p.m, p.n,
         dels_db, delta_s, abs(delta_s),
-        p.d_beta)
+        p.d_beta, delta_norm, delta_visco, delta_dr)
 end
diff --git a/src/Tearing/LayerOverlap.jl b/src/Tearing/LayerOverlap.jl
new file mode 100644
index 000000000..1aac9e572
--- /dev/null
+++ b/src/Tearing/LayerOverlap.jl
@@ -0,0 +1,380 @@
+# LayerOverlap.jl
+#
+# Resistive-layer overlap scan: how far out in ψ the equilibrium domain needs
+# to extend before adjacent rational surfaces' resistive layers run into each
+# other. Once two neighbouring layers overlap, neither surface has a
+# well-separated inner region, so the matched-asymptotic treatment stops being
+# meaningful there and extending `psihigh` past that point buys nothing.
+# Criterion: Fitzpatrick, Nucl. Fusion 2025 (doi 10.1088/1741-4326/ae4fdd) Sect. 5.9 —
+# retain surfaces with psi < 1 - eps_c, where overlap is judged with the
+# diffusive-resistive layer width of its Eq. (100), the regime strong shear forces
+# near the separatrix.
+#
+# Lives in `Tearing` rather than `InnerLayer/SLAYER` for the same reason
+# `build_ggj_inputs` does: it needs `ForceFreeStates._find_rational_surfaces`,
+# and `InnerLayer` loads before `ForceFreeStates`.
+#
+# Surfaces beyond the equilibrium's own ψ grid are located on the separatrix edge
+# law q ≈ -A·ln(1-ψ), via the shared `Equilibrium.edge_q_law` helper that also gates
+# the grid-refinement edge density model -- one statement of the model, one fit, one
+# limited-plasma guard. The log form is the only one that survives the extrapolation:
+# a polynomial in ψ (on q or on ι = 1/q) saturates instead of diverging, undershooting
+# the outermost surfaces enough to report no overlap where there is one.
+#
+# The extrapolated surfaces inform the *choice* of domain only; the final
+# equilibrium is always re-formed and solved on the accepted domain.
+
+using ..Utilities: KineticProfiles
+using ..ForceFreeStates: _find_rational_surfaces, SingType
+using ..InnerLayer.SLAYER: SLAYERParameters, build_slayer_inputs, slayer_layer_thickness,
+    surface_minor_radius, surface_da_dpsi, radial_label
+using ..Utilities.NeoclassicalResistivity: NeoResistivityModel, SauterNeoModel
+using ..Equilibrium: edge_q_law, edge_q_law_psi, edge_q_law_dqdpsi, InverseIngest
+using FastInterpolations: cubic_interp, ExtendExtrap, DerivOp
+using Roots: find_zero, Brent
+using Printf: @sprintf
+
+"""
+    LayerOverlapScan
+
+Per-surface resistive layer widths and the `psihigh` they imply, from
+[`resistive_layer_overlap`](@ref).
+
+All widths are in **normalized flux**, converted from the metre-valued
+`LayerWidths` scales by `Δψ = w / |da/dψ|` so they can be compared against
+surface *separations* in ψ. Comparing the metre values directly against ψ is a
+dimensional error.
+
+# Fields
+
+  - `m`, `n` -- resonant mode numbers per surface, ordered by increasing ψ
+  - `psi` -- surface location in normalized flux
+  - `rs` -- minor radius at the surface in meters
+  - `delta_s_m` -- `|δ_s|`, the Riccati resistive layer width, in meters
+  - `width_delta_s` -- the same width in normalized flux, `|δ_s|/|da/dψ|`
+  - `width_visco` -- the visco-resistive comparison scale as a Δψ width
+  - `width_dr` -- the diffusive-resistive width of Fitzpatrick (2025) Eq. (100) as a Δψ
+    width. This is the criterion the recommendation uses: the paper's Sect. 5.6 shows that
+    strong shear near the separatrix forces every layer there into the DR regime
+  - `extrapolated` -- true when the surface was located outside the
+    equilibrium's ψ grid via the edge q-law
+  - `psihigh_delta_s`, `psihigh_visco`, `psihigh_dr` -- domain implied by each width
+    channel, `nothing` when that channel never overlaps within the scanned range
+  - `psihigh` -- the recommendation, equal to `psihigh_dr`: the DR width is the criterion
+    (Sect. 5.6 puts every near-separatrix layer in that regime). The other two channels are
+    reported for comparison and do not set the domain
+  - `first_overlap` -- index of the first surface that overlaps its inner
+    neighbour, `nothing` when none do
+  - `notes` -- surfaces that were located but could not be scored, and why
+"""
+struct LayerOverlapScan
+    m::Vector{Int}
+    n::Vector{Int}
+    psi::Vector{Float64}
+    rs::Vector{Float64}
+    delta_s_m::Vector{Float64}
+    width_delta_s::Vector{Float64}
+    width_visco::Vector{Float64}
+    width_dr::Vector{Float64}
+    extrapolated::Vector{Bool}
+    psihigh_delta_s::Union{Float64,Nothing}
+    psihigh_visco::Union{Float64,Nothing}
+    psihigh_dr::Union{Float64,Nothing}
+    psihigh::Union{Float64,Nothing}
+    first_overlap::Union{Int,Nothing}
+    notes::Vector{String}
+end
+
+# Inverse equilibria (CHEASE) are solved on a PRESCRIBED boundary: `InverseIngest.sq_xs` spans
+# [0, 1] and `sq_fs[:, 3]` is the code's own q there, finite at ψ = 1 (6.90 on the shipped
+# fixture). So beyond `psihigh` there is nothing to model -- the real q is already known, and the
+# scan reads it instead of extrapolating. That also means the search runs all the way to ψ = 1
+# rather than stopping short of a separatrix that a fixed-boundary equilibrium does not have.
+#
+# The direct/EFIT path cannot do this: a g-file also carries q on [0, 1], but it is the
+# reconstruction's q, which goes to a finite qa where a diverted plasma's q diverges. GPEC
+# recomputes q by field-line tracing out to psihigh, and past that nothing has been traced.
+_inverse_q_spline(ingest::InverseIngest) =
+    cubic_interp(collect(Float64, ingest.sq_xs), collect(Float64, ingest.sq_fs[:, 3]);
+        extrap=ExtendExtrap())
+
+# Locate q = q_target on the real q profile, between `psi_lo` and `psi_hi`. Returns nothing when
+# the target is outside the range this equilibrium actually reaches.
+function _real_q_surface(qspl, q_target::Real, psi_lo::Float64, psi_hi::Float64)
+    q_lo, q_hi = Float64(qspl(psi_lo)), Float64(qspl(psi_hi))
+    (q_lo - q_target) * (q_hi - q_target) <= 0 || return nothing
+    return find_zero(p -> Float64(qspl(p)) - q_target, (psi_lo, psi_hi), Brent())
+end
+
+"""
+    resistive_layer_overlap(equil, profiles::KineticProfiles; n_tor, kwargs...)
+        -> LayerOverlapScan
+
+Locate the q = m/`n_tor` rational surfaces, compute each one's resistive layer
+width, and report the outermost `psihigh` at which adjacent layers are still
+separated.
+
+Surfaces inside the equilibrium grid come from
+`ForceFreeStates._find_rational_surfaces` (Brent root-finding segmented between
+q-extrema, so reverse shear is handled). Surfaces beyond the grid come from the
+separatrix edge law `q ≈ -A·ln(1-ψ)` fitted over the outer knots when
+`extrapolate=true`, which is what lets the scan recommend a domain *larger* than
+the one it was handed. See the file header for why a cubic extrapolation is not
+usable here.
+
+Layer widths come from [`slayer_layer_thickness`](@ref) via
+[`build_slayer_inputs`](@ref), so the scan uses the same plasma inputs and the
+same resistivity closure as the SLAYER analysis itself.
+
+# Arguments
+
+  - `equil` -- a formed `PlasmaEquilibrium`
+  - `profiles` -- `KineticProfiles` spanning the scanned ψ range
+
+# Keyword arguments
+
+  - `n_tor` -- toroidal mode number; surfaces are q = m/`n_tor`
+
+  - `m_max` -- runaway backstop on the poloidal mode number (default `2000`), **not** the
+    physics limit: the outward search should terminate at `psi_cap` (indistinguishable from
+    the separatrix), and if `m_max` is what stops it instead, the scan says so in `notes`
+    because "no overlap found" would then rest on a non-physical cutoff
+  - `max_layer_solves` -- bound on how many surfaces get a Riccati layer solve (default `64`,
+    innermost kept). The first overlap decides the answer, so surfaces far beyond it only cost
+    time; when the bound bites, `notes` records that a "no overlap" verdict is inconclusive
+  - `psihigh_safe` -- treat this as the outer edge of trusted equilibrium data
+    (default `equil.profiles.xs[end]`)
+  - `psi_cap` -- never place a surface beyond this ψ on the **direct** path (default
+    `1 - 1e-9`, purely numerical: it stops the outward search from crowding indefinitely onto
+    ψ = 1 where surfaces become indistinguishable). Inverse equilibria ignore it: they have a
+    prescribed boundary and the search runs to ψ = 1 on the real q. Geometry needs no clamp of
+    its own -- `radial_label` derivatives are analytic and valid under extrapolation
+  - `extrapolate` -- search past `psihigh_safe` on the edge q-law (default `true`)
+  - `theta` -- poloidal angle for the minor-radius chord (default `0.0`)
+  - `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
+  - `mu_i`, `zeff`, `chi_perp`, `chi_tor`, `resistivity_model`, `lnLambda_form`
+    -- passed through to `build_slayer_inputs`
+
+The overlap criterion is that surface `k` overlaps its inner neighbour when
+`ψ_k − w_k/2 < ψ_{k-1} + w_{k-1}/2`. When that happens **both** surfaces are
+contaminated, so the recommendation cuts at the inner edge of surface `k-1`,
+leaving `k-2` as the outermost trustworthy surface.
+"""
+function resistive_layer_overlap(equil, profiles::KineticProfiles;
+    n_tor::Integer,
+    m_max::Integer=2000,
+    max_layer_solves::Integer=64,
+    psihigh_safe::Real=Float64(equil.profiles.xs[end]),
+    psi_cap::Real=1.0 - 1e-9,
+    extrapolate::Bool=true,
+    theta::Real=0.0,
+    rs_method::Symbol=:midplane,
+    mu_i::Real=2.0,
+    zeff::Real=1.0,
+    chi_perp=1.0,
+    chi_tor=1.0,
+    resistivity_model::NeoResistivityModel=SauterNeoModel(),
+    lnLambda_form::Symbol=:nrl)
+
+    n_tor > 0 || throw(ArgumentError("resistive_layer_overlap: n_tor must be positive, got $n_tor"))
+    # Widths come back in metres of whichever radial label `rs_method` selects, so the
+    # conversion into normalised flux must use that same label's Jacobian.
+    _, _dr_dpsi = radial_label(equil; rs_method=rs_method, theta=theta)
+    psihigh_safe = Float64(psihigh_safe)
+    psi_cap = Float64(psi_cap)
+    notes = String[]
+
+    # In-grid surfaces for this n, ordered outward.
+    found = [(m=s.m, psi=s.psifac, extrap=false)
+             for s in _find_rational_surfaces(equil, Int(n_tor), Int(n_tor))
+             if s.psifac <= psihigh_safe && s.m <= m_max]
+    sort!(found; by=s -> s.psi)
+
+    # Surfaces past the grid. An inverse equilibrium already knows its real q out to ψ = 1, so
+    # it is read rather than modelled; only the direct path needs the edge law.
+    edge_fit = nothing
+    inv_q = (getfield(equil, :ingest) isa InverseIngest) ? _inverse_q_spline(equil.ingest) : nothing
+    if extrapolate && inv_q !== nothing
+        psi_top = 1.0                      # a prescribed boundary, not a separatrix
+        m_start = isempty(found) ? 1 : maximum(s.m for s in found) + 1
+        reached_end = false
+        for m in m_start:Int(m_max)
+            psi_m = _real_q_surface(inv_q, m / n_tor, psihigh_safe, psi_top)
+            if psi_m === nothing
+                reached_end = true         # q never reaches m/n inside this equilibrium
+                break
+            end
+            psi_m > psihigh_safe || continue
+            push!(found, (m=m, psi=psi_m, extrap=true))
+        end
+        !reached_end && m_start <= m_max &&
+            push!(notes,
+                "outward search stopped at the m_max = $m_max backstop rather than at the boundary " *
+                "q; a \"no overlap\" result here is not conclusive -- raise m_max")
+        push!(notes, "inverse equilibrium: surfaces beyond ψ = $(@sprintf("%.6f", psihigh_safe)) " *
+                     "read from the real q profile out to ψ = 1 (no edge-law extrapolation)")
+    elseif extrapolate && psihigh_safe < psi_cap
+        edge_fit = edge_q_law(equil; psi_max=psihigh_safe)
+        if edge_fit === nothing
+            push!(
+                notes,
+                "edge q-law rejected: the diverging model q ~ -A*ln(1-ψ) does not describe " *
+                "this edge (limited plasma with finite edge q, q not rising, or too few outer " *
+                "knots), so no surfaces are placed beyond ψ = $(@sprintf("%.6f", psihigh_safe))"
+            )
+        else
+            m_start = isempty(found) ? 1 : maximum(s.m for s in found) + 1
+            # Terminates on psi_cap -- surfaces indistinguishable from the separatrix -- not on
+            # m_max. Reaching m_max means the scan was cut short for a non-physical reason and
+            # any "no overlap" verdict from it is unsafe, so record that.
+            reached_cap = false
+            for m in m_start:Int(m_max)
+                psi_m = edge_q_law_psi(edge_fit, m / n_tor)
+                psi_m > psihigh_safe || continue   # already inside the grid
+                if psi_m >= psi_cap
+                    reached_cap = true             # separatrix reached: the physical end of the scan
+                    break
+                end
+                push!(found, (m=m, psi=psi_m, extrap=true))
+            end
+            if !reached_cap && m_start <= m_max
+                push!(notes,
+                    "outward search stopped at the m_max = $m_max backstop rather than at the " *
+                    "separatrix; a \"no overlap\" result here is not conclusive -- raise m_max")
+            end
+        end
+    end
+
+    isempty(found) && return LayerOverlapScan(Int[], Int[], Float64[], Float64[], Float64[], Float64[],
+        Float64[], Float64[], Bool[], nothing, nothing, nothing, nothing, nothing,
+        push!(notes, "no q = m/$n_tor surfaces found"))
+
+    ms, ns, psis, rss = Int[], Int[], Float64[], Float64[]
+    dels_m, w_dels, w_visc, w_dr, extraps = Float64[], Float64[], Float64[], Float64[], Bool[]
+
+    # Bound the per-surface Riccati work, and say so when the bound bites: a "no overlap"
+    # verdict from a truncated list is not conclusive.
+    if length(found) > max_layer_solves
+        push!(notes,
+            "surface list truncated from $(length(found)) to the innermost $max_layer_solves for " *
+            "the layer solve (max_layer_solves); a \"no overlap\" verdict here is not conclusive")
+        found = found[1:max_layer_solves]
+    end
+
+    for s in found
+        q_val = s.m / n_tor
+        # dq/dψ from the equilibrium spline in-grid, from the edge law outside it.
+        q1_val = if !s.extrap
+            Float64(equil.profiles.q_deriv(s.psi))
+        elseif inv_q !== nothing
+            Float64(inv_q(s.psi; deriv=DerivOp(1)))   # real profile, not a model
+        else
+            edge_q_law_dqdpsi(edge_fit, s.psi)
+        end
+        da_dpsi = _dr_dpsi(s.psi)
+        if !isfinite(da_dpsi) || da_dpsi == 0.0
+            push!(notes, "m=$(s.m): da/dψ = $da_dpsi at ψ=$(@sprintf("%.6f", s.psi)); cannot convert width to flux units")
+            continue
+        end
+
+        sing = SingType(; m=[s.m], n=[Int(n_tor)], psifac=s.psi, rho=sqrt(s.psi), q=q_val, q1=q1_val)
+        # Built one surface at a time so a single degenerate surface (e.g. ω_*e == ω_*i,
+        # 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,
+                mu_i=mu_i, zeff=zeff, chi_perp=chi_perp, chi_tor=chi_tor,
+                dr_val=0.0, dc_type=:none, theta=theta,
+                resistivity_model=resistivity_model, lnLambda_form=lnLambda_form)
+        catch err
+            # Any single-surface failure is recorded and skipped, never fatal. A DomainError is
+            # the common one: profiles that stop short of the extrapolated surfaces, or a
+            # negative-shear surface reaching a sqrt/fractional power in the layer parameters.
+            err isa Union{ArgumentError,DomainError} || rethrow()
+            msg = err isa ArgumentError ? err.msg : sprint(showerror, err)
+            push!(notes, "m=$(s.m) at ψ=$(@sprintf("%.6f", s.psi)): skipped -- $(first(split(msg, '\n')))")
+            continue
+        end
+
+        lw = slayer_layer_thickness(params[1])
+        if !isfinite(lw.delta_s_m)
+            push!(notes, "m=$(s.m) at ψ=$(@sprintf("%.6f", s.psi)): del_s Riccati did not converge")
+            continue
+        end
+        # Every width channel must be finite: a NaN width compares false against every
+        # neighbour in _first_overlap_limit, silently turning that channel into "no overlap".
+        if !isfinite(lw.delta_dr)
+            push!(notes, "m=$(s.m) at ψ=$(@sprintf("%.6f", s.psi)): Eq. (100) width is not finite; surface excluded")
+            continue
+        end
+        if !isfinite(lw.delta_visco)
+            push!(notes, "m=$(s.m) at ψ=$(@sprintf("%.6f", s.psi)): viscous width is not finite; surface excluded")
+            continue
+        end
+
+        # Metres → normalized flux. This conversion is the whole reason the scan
+        # can compare a layer width against a surface separation.
+        push!(ms, s.m)
+        push!(ns, Int(n_tor))
+        push!(psis, s.psi)
+        push!(rss, params[1].rs)
+        push!(dels_m, lw.delta_s_m)
+        push!(w_dels, lw.delta_s_m / abs(da_dpsi))
+        push!(w_visc, lw.delta_visco / abs(da_dpsi))
+        push!(w_dr, lw.delta_dr / abs(da_dpsi))
+        push!(extraps, s.extrap)
+    end
+
+    ph_dels, k_dels = _first_overlap_limit(psis, w_dels)
+    ph_visc, k_visc = _first_overlap_limit(psis, w_visc)
+    ph_dr, k_dr = _first_overlap_limit(psis, w_dr)
+    # The DR width is the criterion: near the separatrix the shear is strong enough that
+    # Fitzpatrick (2025) Sect. 5.6 puts every layer in that regime. The delta_s ODE and the
+    # viscous scale stay reported so the three can be compared, but they do not set the domain.
+    recommended = ph_dr
+
+    return LayerOverlapScan(ms, ns, psis, rss, dels_m, w_dels, w_visc, w_dr, extraps,
+        ph_dels, ph_visc, ph_dr, recommended, k_dr, notes)
+end
+
+# Walk outward and return (recommended psilim, index of first overlapping surface).
+# Overlap at k contaminates both k and k-1, so the limit is surface k-1's own position
+# (Fitzpatrick 2025 Sect. 5.9 retains 0 < psi < 1 - eps_c with psi[k-1] as 1 - eps_c);
+# offsetting by ±w[k-1]/2 either drops the last clean surface or can exceed psi = 1.
+function _first_overlap_limit(psi::Vector{Float64}, w::Vector{Float64})
+    for k in 2:length(psi)
+        if psi[k] - w[k] / 2 < psi[k-1] + w[k-1] / 2
+            return (psi[k-1], k)
+        end
+    end
+    return (nothing, nothing)
+end
+
+function Base.show(io::IO, ::MIME"text/plain", s::LayerOverlapScan)
+    println(io, "LayerOverlapScan: $(length(s.psi)) surface(s), n = ",
+        isempty(s.n) ? "-" : string(s.n[1]))
+    if !isempty(s.psi)
+        println(io, rpad("q", 8), rpad("psi", 13), rpad("r_s [m]", 11),
+            rpad("dpsi(del_s)", 14), rpad("dpsi(visco)", 14), rpad("dpsi(DR)", 14), "source")
+        for i in eachindex(s.psi)
+            println(io,
+                rpad("$(s.m[i])/$(s.n[i])", 8),
+                rpad(@sprintf("%.7f", s.psi[i]), 13),
+                rpad(@sprintf("%.5f", s.rs[i]), 11),
+                rpad(@sprintf("%.4e", s.width_delta_s[i]), 14),
+                rpad(@sprintf("%.4e", s.width_visco[i]), 14),
+                rpad(@sprintf("%.4e", s.width_dr[i]), 14),
+                s.extrapolated[i] ? "extrapolated" : "in-grid")
+        end
+    end
+    _fmt(x) = x === nothing ? "none (no overlap in range)" : @sprintf("%.7f", x)
+    println(io, "  psihigh from |del_s| : ", _fmt(s.psihigh_delta_s))
+    println(io, "  psihigh from visco   : ", _fmt(s.psihigh_visco))
+    println(io, "  psihigh from Eq.(100): ", _fmt(s.psihigh_dr))
+    println(io, "  recommended psihigh  : ", _fmt(s.psihigh))
+    for note in s.notes
+        println(io, "  note: ", note)
+    end
+    return nothing
+end
diff --git a/src/Tearing/Tearing.jl b/src/Tearing/Tearing.jl
index 745a30857..80bce7b23 100644
--- a/src/Tearing/Tearing.jl
+++ b/src/Tearing/Tearing.jl
@@ -22,6 +22,7 @@ using ..Utilities
 import ..InnerLayer as InnerLayer
 
 include("LayerInputs.jl")
+include("LayerOverlap.jl")
 include("Dispersion/Dispersion.jl")
 include("Runner/Runner.jl")
 
@@ -30,5 +31,6 @@ import .Runner as Runner
 
 export InnerLayer, Dispersion, Runner
 export build_ggj_inputs
+export resistive_layer_overlap, LayerOverlapScan
 
 end # module Tearing
diff --git a/test/runtests.jl b/test/runtests.jl
index 7845eea60..540d99a28 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -40,6 +40,7 @@ else
     include("./runtests_slayer_params.jl")
     include("./runtests_slayer_riccati.jl")
     include("./runtests_slayer_inputs.jl")
+    include("./runtests_layer_overlap_fitzpatrick.jl")
     include("./runtests_dispersion_residual.jl")
     include("./runtests_dispersion_coupled.jl")
     include("./runtests_dispersion_coupled_full.jl")
diff --git a/test/runtests_layer_overlap_fitzpatrick.jl b/test/runtests_layer_overlap_fitzpatrick.jl
new file mode 100644
index 000000000..aaf0ea198
--- /dev/null
+++ b/test/runtests_layer_overlap_fitzpatrick.jl
@@ -0,0 +1,193 @@
+# External-reference check of the diffusive-resistive layer width against a published
+# result: Fitzpatrick, Nucl. Fusion 2025 (doi 10.1088/1741-4326/ae4fdd), Sect. 5.8-5.9.
+#
+# The paper reports that for its model JET equilibrium the resistive layers of adjacent
+# rational surfaces first overlap at Psi = 0.9985 for n = 1 (its Fig. 9) and Psi = 0.9952
+# for n = 4 (its Fig. 10). This rebuilds that equilibrium from the paper's own definitions
+# and drives GPEC's `delta_dr` through `slayer_parameters`, so the width formula and the
+# tau_R / tau_A / d_beta chain feeding it are checked against a number in print rather than
+# against GPEC's own history.
+@testset "Layer overlap vs Fitzpatrick (2025) JET model" begin
+    using GeneralizedPerturbedEquilibrium.InnerLayer.SLAYER: slayer_parameters,
+        slayer_layer_thickness, SpitzerHarmModel
+    using QuadGK: quadgk
+    using Roots: find_zero, Bisection
+
+    # Sect. 5.8 parameters.
+    B0, R0, A_MIN = 3.45, 2.96, 1.25
+    ZEFF, MNUM, CHI = 10.0, 2.0, 1.0
+    Q0, Q95, Q105 = 1.01, 3.5, 4.0
+
+    # Model safety factor, Eqs. (32)-(35). alpha- depends on rhat95, which depends on the
+    # profile through Eq. (36), so the pair is solved self-consistently.
+    q_in(r, am) = Q0 - am * log(1 - r^2)
+    q_out(r, ap) = -ap * log(r^2 - 1)
+    q_any(r, am, ap) = r < 1 ? q_in(r, am) : q_out(r, ap)
+    psi_raw(r, am, ap) = quadgk(x -> x / q_any(x, am, ap), 0.0, r; rtol=1e-11)[1]
+
+    am, ap = 1.0, 1.0
+    fp_converged = false
+    for _ in 1:200
+        psi_sep = psi_raw(1.0 - 1e-12, am, ap)
+        r95 = find_zero(r -> psi_raw(r, am, ap) / psi_sep - 0.95, (0.5, 1 - 1e-9), Bisection())
+        f105(r) = (psi_sep + quadgk(x -> x / q_out(x, ap), 1 + 1e-12, r; rtol=1e-10)[1]) / psi_sep - 1.05
+        hi = 1.0 + 1e-9
+        while hi < 1.40 && f105(hi) < 0     # q_out turns negative past r = sqrt(2)
+            hi += 0.005
+        end
+        hi < 1.40 || error("Fitzpatrick JET model: failed to bracket the Psi = 1.05 surface before q_out degenerates")
+        r105 = find_zero(f105, (1 + 1e-9, hi), Bisection())
+        am_new = -(Q95 - Q0) / log(1 - r95^2)
+        ap_new = -Q105 / log(r105^2 - 1)
+        fp_converged = abs(am_new - am) < 1e-12 && abs(ap_new - ap) < 1e-12
+        am, ap = am_new, ap_new
+        fp_converged && break
+    end
+    fp_converged || error("Fitzpatrick JET model: the alpha-/alpha+ fixed point did not converge in 200 iterations")
+    psi_sep = psi_raw(1.0 - 1e-12, am, ap)
+    q_of(r) = q_any(r, am, ap)
+    Psi_of(r) = r < 1 ? psi_raw(r, am, ap) / psi_sep :
+                (psi_sep + quadgk(x -> x / q_out(x, ap), 1 + 1e-12, r; rtol=1e-10)[1]) / psi_sep
+    dq_dr(r) = r < 1 ? am * 2r / (1 - r^2) : -ap * 2r / (r^2 - 1)
+    s_of(r) = r * dq_dr(r) / q_of(r)
+
+    # mtanh edge profiles read off the paper's Fig. 8, anchored at Psi = 0.96 and at the
+    # separatrix (Sect. 5.9 states Te ~ 100 eV there).
+    function make_tanh(f096, f100, f_ped, f_sol)
+        a = (f_ped - f_sol) / 2
+        x1 = atanh(clamp(1 - (f096 - f_sol) / a, -0.999, 0.999))
+        x2 = atanh(clamp(1 - (f100 - f_sol) / a, -0.999, 0.999))
+        d = 0.04 / (x2 - x1)
+        p0 = 0.96 - x1 * d
+        return P -> f_sol + a * (1 - tanh((P - p0) / d))
+    end
+    te_of = make_tanh(400.0, 100.0, 600.0, 20.0)
+    ne_of = make_tanh(3.5e19, 1.2e19, 5.0e19, 3.0e18)
+
+    # Width in units of rhat, through the shipped code path.
+    function width_hat(r, n)
+        P = Psi_of(r)
+        te = te_of(P)
+        ne = ne_of(P)
+        p = slayer_parameters(; n_e=ne, t_e=te, t_i=te,
+            omega=0.0, 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,
+            resistivity_model=SpitzerHarmModel(), lnLambda_form=:nrl)
+        return slayer_layer_thickness(p).delta_dr / A_MIN
+    end
+
+    function surfaces(n; rmin=0.90, rmax=1.35)
+        out = Tuple{Int,Float64}[]
+        for m in 1:400
+            qt = m / n
+            qt < q_of(rmin) && continue
+            r = if qt < q_of(1 - 1e-13)
+                find_zero(x -> q_of(x) - qt, (rmin, 1 - 1e-13), Bisection())
+            elseif qt > q_of(1 + 1e-13) && qt < q_of(rmax)
+                find_zero(x -> q_of(x) - qt, (1 + 1e-13, rmax), Bisection())
+            else
+                nothing
+            end
+            r === nothing || push!(out, (m, r))
+        end
+        sort!(out; by=t -> t[2])
+        return out
+    end
+
+    # Inner boundary of the overlap region: the first surface whose layer runs into its
+    # inner neighbour, reported at that neighbour's location.
+    function overlap_psi(n)
+        S = surfaces(n)
+        for k in 2:length(S)
+            gap = S[k][2] - S[k-1][2]
+            (width_hat(S[k-1][2], n) + width_hat(S[k][2], n)) / 2 >= gap && return Psi_of(S[k-1][2])
+        end
+        return nothing
+    end
+
+    # The self-consistent profile must reproduce the paper's own inputs before the widths
+    # mean anything.
+    @test isapprox(q_of(0.0), Q0; atol=1e-10)
+    @test am > 0 && ap > 0
+
+    psi1 = overlap_psi(1)
+    @test psi1 !== nothing
+    psi1 === nothing || @info "Fitzpatrick 2025 overlap, n = 1" psi=psi1 paper=0.9985 deviation=abs(psi1 - 0.9985) atol=1e-4
+    # Paper: Psi = 0.9985 (n = 1); measured deviation 1.1e-5. The bound covers the pedestal read
+    # off Fig. 8 and the resistivity closure (Spitzer-Harm here against the paper's Eqs. 70-72),
+    # and still leaves ~9x margin. Pure quadrature and root-finding, no BLAS, so platform spread
+    # sits far below this.
+    @test isapprox(psi1, 0.9985; atol=1e-4)
+
+    psi4 = overlap_psi(4)
+    @test psi4 !== nothing
+    psi4 === nothing || @info "Fitzpatrick 2025 overlap, n = 4" psi=psi4 paper=0.9952 deviation=abs(psi4 - 0.9952) atol=1e-3
+    # Paper: Psi = 0.9952 (n = 4); measured deviation 2.3e-4. Looser than n = 1 because these
+    # surfaces sit further in, where the digitised pedestal shape matters more.
+    @test isapprox(psi4, 0.9952; atol=1e-3)
+    # The overlap region must move inward with n, which is the paper's reported scaling.
+    @test psi4 < psi1
+end
+
+# Wiring of the overlap cap into the integration domain. The scan itself is checked against the
+# paper above; this checks that the cap reaches `sing_lim!` and behaves as an upper bound.
+@testset "Layer-overlap cap wiring" begin
+    using GeneralizedPerturbedEquilibrium.ForceFreeStates: sing_lim!, ForceFreeStatesInternal,
+        ForceFreeStatesControl
+    using GeneralizedPerturbedEquilibrium.Equilibrium
+    using TOML
+
+    dir_path = joinpath(dirname(@__DIR__), "examples", "Solovev_ideal_example")
+    inputs = TOML.parsefile(joinpath(dir_path, "gpec.toml"))
+    equil = Equilibrium.setup_equilibrium(
+        Equilibrium.EquilibriumConfig(inputs["Equilibrium"], dir_path),
+        Equilibrium.SolovevConfig(inputs["SOL_INPUT"]))
+
+    _fresh() = (i = ForceFreeStatesInternal(); i.nlow = 1; i.nhigh = 1; i)
+    ctrl = ForceFreeStatesControl(; set_psilim_via_dmlim=false, qhigh=1e3)
+
+    # Baseline: no cap.
+    base = _fresh()
+    sing_lim!(base, ctrl, equil)
+
+    # A cap inside the domain must pull qlim down and move psilim inward.
+    inside = _fresh()
+    cap = 0.5 * (equil.profiles.xs[1] + base.psilim)
+    sing_lim!(inside, ctrl, equil; psilim_cap=cap)
+    @test inside.qlim < base.qlim
+    @test inside.psilim < base.psilim
+    @test isapprox(inside.psilim, cap; rtol=1e-6)
+
+    # A cap beyond psihigh is inert by construction -- the bound never widens the domain.
+    beyond = _fresh()
+    sing_lim!(beyond, ctrl, equil; psilim_cap=base.psilim + 1e-3)
+    @test beyond.qlim == base.qlim
+    @test beyond.psilim == base.psilim
+
+    # nothing is the same as not passing it at all.
+    none = _fresh()
+    sing_lim!(none, ctrl, equil; psilim_cap=nothing)
+    @test none.qlim == base.qlim
+    @test none.psilim == base.psilim
+
+    # The control flag exists and is opt-in.
+    @test ForceFreeStatesControl().psilim_from_layer_overlap == false
+
+    # The no-surfaces path must return an empty scan with a note, not throw (m_max=0 with
+    # extrapolation off forces it deterministically).
+    @testset "resistive_layer_overlap: empty surface list" begin
+        using GeneralizedPerturbedEquilibrium: Tearing, Utilities
+        npsi = 8
+        psi_grid = collect(range(0.0, 1.0; length=npsi))
+        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=zeros(npsi), omega_e=zeros(npsi), omega_i=zeros(npsi))
+        scan = Tearing.resistive_layer_overlap(equil, profs; 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
+        @test any(contains("no q ="), scan.notes)
+    end
+end