diff --git a/.gitignore b/.gitignore
index 7e1ff4f6..0729e730 100644
--- a/.gitignore
+++ b/.gitignore
@@ -16,9 +16,13 @@
 !docs/resources/*.pdf
 *.png
 !docs/src/assets/**/*.png
+# committed documentation figures (generated by docs/src/figures/**/make_*.jl)
+!docs/src/figures/**/*.png
 *.jld2
 .gitattributes
 Manifest.toml
+# the docs figure-provenance manifest is not a Julia Manifest.toml (case-insensitive FS collides)
+!docs/src/figures/manifest.toml
 *.log
 docs/build/
 anaconda_projects/
diff --git a/Project.toml b/Project.toml
index fd2c512d..478b05c1 100644
--- a/Project.toml
+++ b/Project.toml
@@ -25,6 +25,7 @@ PlotlyJS = "f0f68f2c-4968-5e81-91da-67840de0976a"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
@@ -57,14 +58,11 @@ QuadGK = "2.11.3"
 Roots = "2.2.13"
 SparseArrays = "1"
 SpecialFunctions = "2.5.1"
+Random = "1"
 StaticArrays = "1.9.15"
 Statistics = "1"
 TOML = "1"
 Test = "1"
 julia = "1.11"
 
-[extras]
-Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
-[targets]
-test = ["Random"]
diff --git a/docs/DOC_STANDARD.md b/docs/DOC_STANDARD.md
new file mode 100644
index 00000000..a1af2043
--- /dev/null
+++ b/docs/DOC_STANDARD.md
@@ -0,0 +1,126 @@
+# GPEC documentation standard
+
+This is the authoring standard for the GPEC user documentation under
+`docs/src/`. The goal is that each physics module page reads like a compact
+*Journal of Computational Physics* methods paper: the governing equations, the
+numerical method, and the validation that the method works — with figures whose
+provenance is explicit. Those module pages **are** built by Documenter and
+published on the public documentation site.
+
+This standard file (`docs/DOC_STANDARD.md`) is itself the one exception: it lives
+outside `docs/src/` and is not listed in `make.jl`, so Documenter does not
+render it to the website. It is a contributor-facing reference, read here in the
+repository — which is why `developer_notes.md` links to it by GitHub URL rather
+than as a site-internal page.
+
+The `Ballooning and Mercier Local Stability` and `Inner Layer Module` pages are
+the reference exemplars; new or reworked module pages should follow their shape.
+
+## Module-page template
+
+A module reference page (`docs/src/<module>.md`) has these sections, in order.
+Not every section is large for every module, but the skeleton is the same.
+
+1. **Overview** — what the module computes and where it sits in the
+   equilibrium → stability → perturbed-equilibrium pipeline.
+2. **Governing equations** — the physics the module discretizes, written as
+   numbered display equations. Each equation or block cites the paper it comes
+   from (the PDFs in `docs/resources/`), e.g. *(GWP2016 Eq. 11)*. Citing
+   equation numbers is required, not optional — it is what makes the
+   implementation auditable against theory.
+3. **Numerical method** — the discretization and algorithm: what is expanded in
+   what basis, how the linear/eigenvalue/BVP problem is assembled and solved,
+   and the *design choices* that matter (why this contour, this ordering, this
+   precision). This is the part that distinguishes a methods page from an API
+   dump.
+4. **Validation & benchmarks** — evidence the method is correct: a table of
+   pinned values cross-checked against an independent code or analytic result,
+   and figures (see below) showing convergence, invariance, or agreement.
+5. **Practical usage** — the configuration keys and the public entry points a
+   user calls, with a short worked example.
+6. **API Reference** — an `@autodocs` block for the module so every exported
+   docstring is rendered and `checkdocs=:exports` stays satisfied.
+
+## Figures
+
+### Where figures live
+
+- **Content figures** (anything shown on a page) live in
+  `docs/src/figures/<module>/`. Each committed `<name>.png` sits next to the
+  script `make_<name>.jl` that produced it — script and figure are always
+  committed together, so it is always clear how a figure was made.
+- `docs/src/assets/` is reserved for **Documenter chrome only** (the site logo,
+  themes, custom CSS). Content figures never go there.
+- The shared helper `docs/figure_tools.jl` (outside `src/`, so it is not
+  published) provides `save_doc_figure`, provenance stamping, the manifest, and
+  the `step_series` spectrum helper. Every `make_<name>.jl` `include`s it.
+
+### Where a figure's data comes from
+
+Documentation figures get their **own** generator scripts under
+`docs/src/figures/`. They are not forced to route through `benchmarks/`. When a
+figure genuinely *is* a benchmark or Fortran cross-check result, its
+`make_<name>.jl` may call the relevant `benchmarks/` script, and its `depends`
+list records that coupling — but the committed PNG and a runnable script still
+live together under `docs/src/figures/<module>/`.
+
+`benchmarks/` remains the home of performance comparisons and Fortran
+cross-checks (its outputs are git-ignored and may be heavy). `docs/src/figures/`
+is the home of the small, purpose-built, committed figures that illustrate a
+page.
+
+### Figures do not build with the docs
+
+Figure scripts run **manually**, in the root project:
+
+```bash
+julia --project=. docs/src/figures/<module>/make_<name>.jl
+```
+
+They are **never** run at Documenter build time. The build only embeds the
+committed PNGs, so `docs/Project.toml` stays minimal (no `Plots`) and the build
+stays fast. This is also why we do not regenerate every figure on every docs
+change — see the regeneration policy below.
+
+### Provenance
+
+`save_doc_figure` stamps a small `GPEC <shorthash> · <date>` mark in the corner
+of every figure (so a figure's age is visible when browsing the rendered site)
+and writes a machine-readable entry to `docs/src/figures/manifest.toml`:
+
+```toml
+[inner_layer.rotated_ray_contour]
+script  = "make_rotated_ray_contour.jl"
+commit  = "3a0837e3"
+date    = "2026-07-08"
+depends = ["src/InnerLayer/GGJ/Ray.jl"]
+```
+
+`depends` lists the repo-relative source files whose *numbers* the figure
+visualizes. Generate figures as the last step before committing a change; the
+stamp then reflects the state they were made against (a `-dirty` suffix marks a
+figure generated with an uncommitted working tree, which is honest, not an
+error).
+
+### When to regenerate a figure
+
+Regenerate a figure (and its manifest entry) only when:
+
+1. a file in its `depends` list changed the numbers it shows — the regression
+   harness is the trigger: if a tracked quantity for that module moved, its
+   figures are suspect;
+2. the physics or method the figure illustrates changed; or
+3. the figure script itself changed.
+
+Do **not** regenerate for prose edits, formatting, or changes to unrelated
+modules. Because each figure carries its generating commit and a `depends`
+list, a reviewer (or a future `docs/check_figures.jl`) can tell when a figure
+predates the last change to a file it depends on and is therefore stale.
+
+### Legacy figures
+
+A figure migrated from before this system, whose original generator was not
+preserved and which cannot be faithfully reproduced from current code, is
+registered with `register_legacy_figure` (`commit = "legacy"`, a `note`
+explaining why). This is honest bookkeeping for genuinely irreproducible
+figures — it is never a substitute for writing a generator for a new figure.
diff --git a/docs/figure_tools.jl b/docs/figure_tools.jl
new file mode 100644
index 00000000..ca8001be
--- /dev/null
+++ b/docs/figure_tools.jl
@@ -0,0 +1,166 @@
+# figure_tools.jl
+#
+# Shared helpers for GPEC documentation figures. Every figure generator
+# (`docs/src/figures/<module>/make_<name>.jl`) `include`s this file and calls
+# `save_doc_figure` to stamp, save, and register its output. See
+# `docs/DOC_STANDARD.md` for the figure-organization and provenance policy.
+#
+# This file lives OUTSIDE `docs/src/` so Documenter does not publish it. The
+# `make_<name>.jl` scripts do live under `docs/src/figures/` (co-located with
+# their PNGs) and are published alongside the figures for reproducibility.
+#
+# Figure scripts run MANUALLY in the root project, never at docs-build time:
+#     julia --project=. docs/src/figures/<module>/make_<name>.jl
+# so the Documenter build stays fast and Plots-free.
+
+using Plots
+using Dates
+using TOML
+
+const _DOCS_DIR = @__DIR__                                   # .../docs
+const _REPO_ROOT = normpath(joinpath(_DOCS_DIR, ".."))       # repo root
+const FIGURES_ROOT = joinpath(_DOCS_DIR, "src", "figures")   # docs/src/figures
+const MANIFEST_PATH = joinpath(FIGURES_ROOT, "manifest.toml")
+
+"""
+    figure_provenance() -> (; commit, date, dirty)
+
+Short git commit and ISO date stamped onto every figure. `commit` is the
+`git rev-parse --short HEAD` hash (reusing the idiom from
+`benchmarks/benchmark_git_branches.jl`), with a `-dirty` suffix when the
+working tree has uncommitted changes so a figure never claims a cleaner
+provenance than it has. Falls back to `"unknown"` outside a git checkout.
+"""
+function figure_provenance()
+    date = string(Dates.today())
+    try
+        h = strip(read(`git -C $(_REPO_ROOT) rev-parse --short HEAD`, String))
+        dirty = !isempty(strip(read(`git -C $(_REPO_ROOT) status --porcelain`, String)))
+        return (; commit=(dirty ? "$h-dirty" : h), date=date, dirty=dirty)
+    catch
+        return (; commit="unknown", date=date, dirty=false)
+    end
+end
+
+"""
+    stamp!(p; commit, date, subplot=length(p.subplots))
+
+Annotate a small `GPEC <commit> · <date>` provenance mark in the bottom-right
+corner of subplot `subplot` (default: the last panel), so the figure's age is
+visible when browsing the rendered docs. Works on linear or log axes
+(positions in data coordinates). Never throws: a plotting hiccup warns and
+leaves the figure unstamped rather than losing the plot.
+"""
+function stamp!(p::Plots.Plot; commit::AbstractString, date::AbstractString,
+    subplot::Int=length(p.subplots))
+    try
+        sp = p[subplot]
+        xl = Plots.xlims(sp)
+        yl = Plots.ylims(sp)
+        xpos = xl[2]                                          # right edge (data coords)
+        ypos = yl[1]                                          # bottom edge
+        txt = Plots.text("GPEC $(commit) · $(date)", 6, RGBA(0.5, 0.5, 0.5, 0.9),
+            :right, :bottom)
+        annotate!(sp, xpos, ypos, txt)
+    catch err
+        @warn "stamp! failed; saving figure without provenance mark" exception = err
+    end
+    return p
+end
+
+"""
+    save_doc_figure(p, mod, name; script, depends, npx...) -> String
+
+Stamp `p` with the current git provenance, save it to
+`docs/src/figures/<mod>/<name>.png`, and upsert its `manifest.toml` entry
+(`script`, `commit`, `date`, `depends`). `depends` lists the repo-relative
+source files whose numbers the figure visualizes — the regeneration policy in
+`docs/DOC_STANDARD.md` keys off them. Returns the absolute PNG path (printed
+so it can be opened directly). Extra keyword args are forwarded to `savefig`.
+"""
+function save_doc_figure(p::Plots.Plot, mod::AbstractString, name::AbstractString;
+    script::AbstractString, depends::AbstractVector{<:AbstractString}=String[])
+    prov = figure_provenance()
+    stamp!(p; commit=prov.commit, date=prov.date)
+
+    outdir = joinpath(FIGURES_ROOT, mod)
+    mkpath(outdir)
+    outpath = joinpath(outdir, name * ".png")
+    savefig(p, outpath)
+
+    _update_manifest!(mod, name; script=String(script), commit=prov.commit,
+        date=prov.date, depends=String.(depends))
+
+    println("Saved doc figure: $(abspath(outpath))")
+    return abspath(outpath)
+end
+
+# Load / merge-write the figures manifest, keyed [<mod>.<name>]. Kept sorted and
+# human-diffable; one entry per committed figure.
+function _load_manifest()
+    return isfile(MANIFEST_PATH) ? TOML.parsefile(MANIFEST_PATH) : Dict{String,Any}()
+end
+
+function _write_manifest(manifest)
+    open(MANIFEST_PATH, "w") do io
+        println(io, "# GPEC documentation-figure provenance manifest.")
+        println(io, "# One [<module>.<figure>] entry per committed PNG under docs/src/figures/.")
+        println(io, "# Regenerate a figure (and this entry) only when a file in its `depends`")
+        println(io, "# list changes the numbers it shows — see docs/DOC_STANDARD.md.")
+        println(io)
+        TOML.print(io, manifest; sorted=true)
+    end
+    return manifest
+end
+
+function _update_manifest!(mod::AbstractString, name::AbstractString;
+    script::String, commit::String, date::String, depends::Vector{String})
+    manifest = _load_manifest()
+    modtbl = get!(manifest, mod, Dict{String,Any}())
+    modtbl[name] = Dict{String,Any}(
+        "script" => script,
+        "commit" => commit,
+        "date" => date,
+        "depends" => depends
+    )
+    return _write_manifest(manifest)
+end
+
+"""
+    register_legacy_figure(mod, name; note, depends=String[])
+
+Record a manifest entry for a figure that predates this provenance system and
+whose original generator script was not preserved (`commit = "legacy"`). Use
+this only for migrated figures that cannot be faithfully reproduced from the
+current code — never to skip writing a generator for a new figure. `note`
+states why the figure is legacy (e.g. "compared the since-removed Mercier.jl").
+"""
+function register_legacy_figure(mod::AbstractString, name::AbstractString;
+    note::AbstractString, depends::AbstractVector{<:AbstractString}=String[])
+    manifest = _load_manifest()
+    modtbl = get!(manifest, mod, Dict{String,Any}())
+    modtbl[name] = Dict{String,Any}(
+        "script" => "(not preserved)",
+        "commit" => "legacy",
+        "date" => "unknown",
+        "depends" => String.(depends),
+        "note" => String(note)
+    )
+    _write_manifest(manifest)
+    println("Registered legacy figure: $(mod).$(name)")
+    return manifest
+end
+
+"""
+    step_series(m_vals, amps) -> (m_ext, amp_ext)
+
+Canonical spectrum-plot helper (CLAUDE.md plotting convention): pad a zero on
+each end of a discrete-mode series so a `seriestype=:steppre` plot draws clean
+boxes that fall to zero at the edges. Centralized here so doc figures and
+benchmarks share one definition.
+"""
+function step_series(m_vals, amps)
+    m_ext = [m_vals[1] - 1; m_vals; m_vals[end] + 1]
+    amp_ext = [0.0; amps; 0.0]
+    return m_ext, amp_ext
+end
diff --git a/docs/src/ballooning.md b/docs/src/ballooning.md
index 6207b225..3313e535 100644
--- a/docs/src/ballooning.md
+++ b/docs/src/ballooning.md
@@ -854,7 +854,7 @@ numerically fragile. The comparison below shows that the `Baloo.f` style, which
 sets the displacement to vanish at ``\pm\theta_{\max}`` and computes ``\Delta'``
 from the resulting solution, gives the most stable continuous profile.
 
-![Comparison of numerical stability among methods for extracting ``\Delta'`` and ``c_{a1}``.](assets/ballooning/delta_prime_matching.png)
+![Comparison of numerical stability among methods for extracting ``\Delta'`` and ``c_{a1}``.](figures/ballooning/delta_prime_matching.png)
 
 `Bal.jl` therefore adopts the `Baloo.f`-style boundary replacement for the
 ballooning ``\Delta'`` calculation.
@@ -905,7 +905,7 @@ The previous standalone `Mercier.jl` path can therefore be removed from the loca
 stability calculation. The comparison below shows agreement between the previous
 Mercier calculation and the new `Bal.jl` calculation.
 
-![Comparison of Mercier criterion calculations from `Mercier.jl` and `Bal.jl`.](assets/ballooning/mercier_comparison.png)
+![Comparison of Mercier criterion calculations from `Mercier.jl` and `Bal.jl`.](figures/ballooning/mercier_comparison.png)
 
 The local-stability output now stores ballooning ``\Delta'`` in the fourth
 `locstab_fs` entry. In the HDF5 output this is written as
diff --git a/docs/src/citations.md b/docs/src/citations.md
index dfb65fb6..143b9a49 100644
--- a/docs/src/citations.md
+++ b/docs/src/citations.md
@@ -104,7 +104,7 @@ Provides an efficient algorithm for computing the full Δ' matrix (coupling betw
 > *Physics of Plasmas* **27**, 012506 (2020).
 > DOI: [10.1063/1.5134999](https://doi.org/10.1063/1.5134999)
 
-Showcases a different basis that significantly aids convergence of the Galerkin solver used in the GGJ InnerLayer module.
+Constructs the Wasow large-``x`` asymptotic basis (the `inps` kernel) that supplies the far-field boundary condition for every GGJ InnerLayer backend — the Galerkin solver and the rotated-contour collocation (`:ray`) solver alike.
 
 ---
 
diff --git a/docs/src/developer_notes.md b/docs/src/developer_notes.md
index 9a77b751..95f586da 100644
--- a/docs/src/developer_notes.md
+++ b/docs/src/developer_notes.md
@@ -16,6 +16,18 @@ CODE - TAG - Detailed message
 
 Where CODE is the module name (EQUIL, ForceFreeStates, VAC, PERTURBED EQUILIBRIUM, etc.) and TAG describes the type of change (WIP, MINOR, IMPROVEMENT, BUG FIX, NEW FEATURE, REFACTOR, CLEANUP, etc.). This format is used for compiling release notes — tags should be human-readable but are not enforced to a fixed set.
 
+## Documentation standard
+
+Module reference pages follow a shared *Journal of Computational Physics*-style
+structure (governing equations → numerical method → validation figures → API),
+and every documentation figure is committed together with the script that made
+it and a provenance stamp. The full policy — page template, figure
+organization under `docs/src/figures/<module>/`, provenance, and the
+regenerate-only-when-`depends`-change rule — is in
+[`docs/DOC_STANDARD.md`](https://github.com/OpenFUSIONToolkit/GPEC/blob/develop/docs/DOC_STANDARD.md).
+The `Ballooning and Mercier Local Stability` and `Inner Layer Module` pages are
+the reference exemplars.
+
 ## Regression Testing
 
 The regression harness **must be run on every pull request before merging into `develop`**. It is the project's primary safeguard for tracking how numerical results evolve across changes, so it is only useful if every PR exercises it. When you open a PR, paste the regression report into the PR thread so reviewers can see what moved (and what did not). If your change touches a quantity that is not yet tracked, add a new regression case — or extend an existing one — in the same PR.
diff --git a/docs/src/assets/ballooning/delta_prime_matching.png b/docs/src/figures/ballooning/delta_prime_matching.png
similarity index 100%
rename from docs/src/assets/ballooning/delta_prime_matching.png
rename to docs/src/figures/ballooning/delta_prime_matching.png
diff --git a/docs/src/assets/ballooning/mercier_comparison.png b/docs/src/figures/ballooning/mercier_comparison.png
similarity index 100%
rename from docs/src/assets/ballooning/mercier_comparison.png
rename to docs/src/figures/ballooning/mercier_comparison.png
diff --git a/docs/src/figures/inner_layer/backend_regime_map.png b/docs/src/figures/inner_layer/backend_regime_map.png
new file mode 100644
index 00000000..d86221bc
Binary files /dev/null and b/docs/src/figures/inner_layer/backend_regime_map.png differ
diff --git a/docs/src/figures/inner_layer/convergence_Sinvariance.png b/docs/src/figures/inner_layer/convergence_Sinvariance.png
new file mode 100644
index 00000000..d8b55082
Binary files /dev/null and b/docs/src/figures/inner_layer/convergence_Sinvariance.png differ
diff --git a/docs/src/figures/inner_layer/make_backend_regime_map.jl b/docs/src/figures/inner_layer/make_backend_regime_map.jl
new file mode 100644
index 00000000..163b5049
--- /dev/null
+++ b/docs/src/figures/inner_layer/make_backend_regime_map.jl
@@ -0,0 +1,65 @@
+# make_backend_regime_map.jl
+#
+# Figure: accuracy of the three GGJ backends along the imaginary-Q axis. The
+# :ray backend is the reference (its own error stays below ~1e-5 out to
+# |Q| = 500 — see the convergence figure). The relative error of :shooting and
+# :galerkin against it is plotted versus |Q|: :shooting holds to |Q| ~ 1,
+# :galerkin to |Q| ~ 4, and both then lose all accuracy, while :ray continues.
+# The crossovers of each curve with the 1% line are the practical reach of each
+# method — this is what motivates the :ray backend. A few tens of solves; ~1 min.
+#
+# Run manually:  julia --project=. docs/src/figures/inner_layer/make_backend_regime_map.jl
+
+include(joinpath(@__DIR__, "..", "..", "..", "figure_tools.jl"))
+
+using GeneralizedPerturbedEquilibrium
+using LinearAlgebra
+using Printf
+
+const IL = GeneralizedPerturbedEquilibrium.InnerLayer
+const GGJ = IL.GGJ
+
+p = IL.q4_surface_benchmark()
+γ(aq) = im * aq * GGJ.q0(p)                       # imaginary axis: Q = i|Q|
+
+# Worst-case relative error over the two parity components vs the :ray reference.
+relerr(Δ, Δref) = maximum(abs.(Δ .- Δref) ./ abs.(Δref))
+
+absQ = 10 .^ range(-1, log10(500); length=22)
+ref = [IL.solve_inner(IL.GGJModel(; solver=:ray), p, γ(aq)) for aq in absQ]
+
+gal_x, gal_e = Float64[], Float64[]
+shoot_x, shoot_e = Float64[], Float64[]
+for (k, aq) in enumerate(absQ)
+    try
+        Δ = IL.solve_inner(IL.GGJModel(; solver=:galerkin), p, γ(aq))
+        all(isfinite, Δ) && (push!(gal_x, aq); push!(gal_e, max(relerr(Δ, ref[k]), 1e-16)))
+    catch
+    end
+    aq > 10 && continue                            # :shooting is meaningless past its regime
+    try
+        Δ = IL.solve_inner(IL.GGJModel(; solver=:shooting), p, γ(aq))
+        all(isfinite, Δ) && (push!(shoot_x, aq); push!(shoot_e, max(relerr(Δ, ref[k]), 1e-16)))
+    catch
+    end
+end
+
+@printf("computed %d :ray refs, %d :galerkin, %d :shooting points\n",
+    length(absQ), length(gal_x), length(shoot_x))
+
+plt = plot(; size=(820, 560), legend=:bottomright, xscale=:log10, yscale=:log10,
+    xlabel="|Q|  (imaginary axis, Q = i|Q|)",
+    ylabel="relative error vs :ray reference",
+    title="Backend accuracy along the imaginary-Q axis  (q=4 surface)",
+    ylims=(1e-7, 3e0), left_margin=9Plots.mm, bottom_margin=5Plots.mm)
+
+plot!(plt, shoot_x, shoot_e; lw=2, marker=:diamond, ms=5, color=:seagreen, label=":shooting")
+plot!(plt, gal_x, gal_e; lw=2, marker=:circle, ms=5, color=:orange, label=":galerkin")
+hline!(plt, [1e-2]; color=:red, ls=:dash, lw=1.5, label="1% accuracy")
+annotate!(plt, 8.0, 3e-6,
+    Plots.text(":ray reference\n(error < 1e-5 to |Q| = 500)", 8, :center, RGB(0.2, 0.3, 0.6)))
+
+save_doc_figure(plt, "inner_layer", "backend_regime_map";
+    script=basename(@__FILE__),
+    depends=["src/InnerLayer/GGJ/Ray.jl", "src/InnerLayer/GGJ/Galerkin.jl",
+        "src/InnerLayer/GGJ/Shooting.jl"])
diff --git a/docs/src/figures/inner_layer/make_convergence_Sinvariance.jl b/docs/src/figures/inner_layer/make_convergence_Sinvariance.jl
new file mode 100644
index 00000000..3e2a71b2
--- /dev/null
+++ b/docs/src/figures/inner_layer/make_convergence_Sinvariance.jl
@@ -0,0 +1,45 @@
+# make_convergence_Sinvariance.jl
+#
+# Figure: honest numerical error bar on the matching data Δ at Q = 500i. The
+# `delta_convergence` battery re-solves with each numerical knob perturbed on
+# an independent axis (contour angle θ, spectral order p, series order/radius,
+# refinement depth, march tolerance, handoff radius, purification e-folds). Δ
+# is an analytic invariant of the contour, so the worst-case spread across
+# these orthogonal perturbations bounds the true error. Eight solves; ~30 s.
+#
+# Run manually:  julia --project=. docs/src/figures/inner_layer/make_convergence_Sinvariance.jl
+
+include(joinpath(@__DIR__, "..", "..", "..", "figure_tools.jl"))
+
+using GeneralizedPerturbedEquilibrium
+using Printf
+
+const IL = GeneralizedPerturbedEquilibrium.InnerLayer
+
+p = IL.q4_surface_benchmark()
+Q = 500.0im
+
+conv = IL.delta_convergence(p, Q; verbose=false)
+names = [r.name for r in conv.table]
+d1 = [max(r.d1, 1e-16) for r in conv.table]      # relative change of Δ_odd
+d2 = [max(r.d2, 1e-16) for r in conv.table]      # relative change of Δ_even
+n = length(names)
+
+plt = plot(; size=(820, 540), legend=:topright, yscale=:log10,
+    xticks=(1:n, names), xrotation=30, ylims=(1e-12, 1e-2),
+    ylabel="relative change of Δ per knob", xlabel="perturbed numerical knob",
+    title=@sprintf("Convergence & contour-invariance of Δ,  Q = %gi  (q=4)", imag(Q)),
+    left_margin=10Plots.mm, bottom_margin=14Plots.mm)
+
+scatter!(plt, 1:n, d1; marker=:circle, ms=7, color=1, label="δΔ_odd")
+scatter!(plt, 1:n, d2; marker=:diamond, ms=7, color=2, label="δΔ_even")
+
+# Worst-case spread — the reported error bar on each parity.
+hline!(plt, [conv.spread[1]]; color=1, ls=:dash, lw=1.5,
+    label=@sprintf("worst-case Δ_odd  %.1e", conv.spread[1]))
+hline!(plt, [conv.spread[2]]; color=2, ls=:dash, lw=1.5,
+    label=@sprintf("worst-case Δ_even  %.1e", conv.spread[2]))
+
+save_doc_figure(plt, "inner_layer", "convergence_Sinvariance";
+    script=basename(@__FILE__),
+    depends=["src/InnerLayer/GGJ/Ray.jl"])
diff --git a/docs/src/figures/inner_layer/make_rotated_ray_contour.jl b/docs/src/figures/inner_layer/make_rotated_ray_contour.jl
new file mode 100644
index 00000000..15abe167
--- /dev/null
+++ b/docs/src/figures/inner_layer/make_rotated_ray_contour.jl
@@ -0,0 +1,65 @@
+# make_rotated_ray_contour.jl
+#
+# Figure: the rotated integration ray in the complex layer-coordinate (x) plane.
+# Shows why the :ray backend rotates the contour by θ = arg(Q)/4: the real-axis
+# contour (θ = 0) runs straight through the pseudo-resonance at
+# x² = −Q²(G + K F) (real and large on the imaginary-Q axis), while the rotated
+# ray clears it. Geometry only — no BVP solve — so it is cheap to regenerate.
+#
+# Run manually:  julia --project=. docs/src/figures/inner_layer/make_rotated_ray_contour.jl
+
+include(joinpath(@__DIR__, "..", "..", "..", "figure_tools.jl"))
+
+using GeneralizedPerturbedEquilibrium
+using Printf
+
+const IL = GeneralizedPerturbedEquilibrium.InnerLayer
+const GGJ = IL.GGJ
+
+# q = 4 rational-surface benchmark on the imaginary-Q axis (the regime that
+# defeats :galerkin). Q is the scaled growth rate; θ makes the parabolic
+# exponent real.
+p = IL.q4_surface_benchmark()
+Q = 500.0im
+θ = angle(Q) / 4                                   # = π/8 = 22.5°
+
+# Pseudo-resonance location: x² = −Q²(G + K F). On the imaginary-Q axis this is
+# real and large, so it sits ON the real-x contour.
+x_pr = sqrt(-Q^2 * (p.G + p.K * p.F))
+xpr = real(x_pr)                                   # imaginary part ≈ 0 here
+
+# Matching radius from the series-residual criterion along the ray (annotation).
+S, _, ok = IL.pick_smax(p, Q; θ=θ)
+
+# Draw both contours out to ~1.8× the pseudo-resonance radius so it is in frame.
+smax = 1.8 * xpr
+s = range(0, smax; length=400)
+ray = cis(θ) .* s                                  # x = e^{iθ} s
+
+plt = plot(; size=(720, 560), legend=:topleft, framestyle=:zerolines,
+    xlabel="Re x", ylabel="Im x", aspect_ratio=:equal,
+    title="Rotated integration ray, Q = 500i  (q=4 surface)",
+    left_margin=6Plots.mm, bottom_margin=4Plots.mm)
+
+# Real-axis contour (θ = 0) and the rotated ray.
+plot!(plt, [0, smax], [0, 0]; lw=2, ls=:dash, color=:gray,
+    label="real-axis contour (θ = 0)")
+plot!(plt, real.(ray), imag.(ray); lw=3, color=1,
+    label=@sprintf("rotated ray  x = e^{iθ}s,  θ = %.1f°", rad2deg(θ)))
+
+# Pseudo-resonance on the real axis — the point the rotation steps around.
+scatter!(plt, [xpr], [0.0]; marker=:xcross, ms=9, color=:red, msw=3,
+    label="pseudo-resonance  x² = −Q²(G+KF)")
+
+# Clearance note above the ray near the pseudo-resonance radius, plus the true
+# matching radius. Positions are tied to xpr so they stay inside the frame.
+annotate!(plt, 1.02 * xpr, 0.56 * xpr,
+    Plots.text("ray stays clear of the\nreal-axis pseudo-resonance", 8, :left,
+        RGB(0.15, 0.15, 0.15)))
+annotate!(plt, 0.45 * xpr, -0.14 * xpr,
+    Plots.text(@sprintf("matching radius S ≈ %.2g%s", S, ok ? "" : " (series tol not met)"),
+        8, :gray, :left))
+
+save_doc_figure(plt, "inner_layer", "rotated_ray_contour";
+    script=basename(@__FILE__),
+    depends=["src/InnerLayer/GGJ/Ray.jl"])
diff --git a/docs/src/figures/inner_layer/make_solution_profiles.jl b/docs/src/figures/inner_layer/make_solution_profiles.jl
new file mode 100644
index 00000000..2165f358
--- /dev/null
+++ b/docs/src/figures/inner_layer/make_solution_profiles.jl
@@ -0,0 +1,60 @@
+# make_solution_profiles.jl
+#
+# Figure: the inner-layer fields (Ψ, Ξ, Υ) reconstructed on the rotated ray for
+# the q=4 benchmark surface. The collocation solution on [0, s_m] (solid) and
+# the analytic u_small + Δ·u_big asymptotic representation for s ≥ S (dashed)
+# share the same power-law tail — the numeric↔asymptotic overlap the outer-
+# region matching relies on. One BVP solve; a few minutes at large imaginary Q.
+#
+# Run manually:  julia --project=. docs/src/figures/inner_layer/make_solution_profiles.jl
+
+include(joinpath(@__DIR__, "..", "..", "..", "figure_tools.jl"))
+
+using GeneralizedPerturbedEquilibrium
+using Printf
+
+const IL = GeneralizedPerturbedEquilibrium.InnerLayer
+const GGJ = IL.GGJ
+
+p = IL.q4_surface_benchmark()
+Q = 2.0im                  # imaginary axis, beyond the |Q|≪1 shooting regime,
+#                            small enough that the collocation domain reaches
+#                            the series radius directly (no march) — so the
+#                            numeric and asymptotic segments join seamlessly.
+isol = 1                   # "odd" parity solution (Ψ'(0)=Ξ(0)=Υ(0)=0)
+
+res = IL.solve_ray(p, Q)
+s_m = res.breaks[end]
+@printf("solve_ray: Q=%s  θ=%.3f  S=%.4g  s_m=%.4g  resid=%.1e\n",
+    string(Q), res.θ, res.S, s_m, res.resid)
+
+# Collocation solution over the BVP domain [0, s_m].
+prof = IL.solution_profile(res; npc=8)
+
+# Analytic tail on [S, few·S] where the inps series is trusted.
+srange = 10 .^ range(log10(res.S), log10(4 * res.S); length=80)
+asy = IL.asymptotic_profile(p, res, srange)
+
+fields = ((:Ψ, prof.Ψ, asy.Ψ, 1), (:Ξ, prof.Ξ, asy.Ξ, 2), (:Υ, prof.Υ, asy.Υ, 3))
+
+plt = plot(; size=(760, 560), legend=:bottomleft, xscale=:log10, yscale=:log10,
+    ylims=(1e-3, 1e2),
+    xlabel="ray parameter s   (x = e^{iθ}s)", ylabel="|field|  (odd parity)",
+    title=@sprintf("Inner-layer fields on the rotated ray, Q = %gi  (q=4)", imag(Q)),
+    left_margin=8Plots.mm, bottom_margin=5Plots.mm)
+
+for (nm, col, acol, ci) in fields
+    plot!(plt, prof.s[2:end], abs.(col[2:end, isol]); lw=2.5, color=ci,
+        label="|$(nm)|  collocation")
+    plot!(plt, asy.s, abs.(acol[:, isol]); lw=2, ls=:dash, color=ci, label="|$(nm)|  asymptotic")
+end
+
+# Matching radius: with no march the BVP edge s_m coincides with the series
+# radius S, so numeric and asymptotic meet at a single point.
+vline!(plt, [res.S]; color=:black, ls=:dot, lw=1,
+    label=@sprintf("S = s_m ≈ %.1f  (match point)", res.S))
+
+save_doc_figure(plt, "inner_layer", "solution_profiles";
+    script=basename(@__FILE__),
+    depends=["src/InnerLayer/GGJ/Ray.jl", "src/InnerLayer/GGJ/RayAsymptotics.jl",
+        "src/InnerLayer/GGJ/InnerAsymptotics.jl"])
diff --git a/docs/src/figures/inner_layer/rotated_ray_contour.png b/docs/src/figures/inner_layer/rotated_ray_contour.png
new file mode 100644
index 00000000..27884cf1
Binary files /dev/null and b/docs/src/figures/inner_layer/rotated_ray_contour.png differ
diff --git a/docs/src/figures/inner_layer/solution_profiles.png b/docs/src/figures/inner_layer/solution_profiles.png
new file mode 100644
index 00000000..3a1c2ef5
Binary files /dev/null and b/docs/src/figures/inner_layer/solution_profiles.png differ
diff --git a/docs/src/figures/manifest.toml b/docs/src/figures/manifest.toml
new file mode 100644
index 00000000..2ad9ef5d
--- /dev/null
+++ b/docs/src/figures/manifest.toml
@@ -0,0 +1,43 @@
+# GPEC documentation-figure provenance manifest.
+# One [<module>.<figure>] entry per committed PNG under docs/src/figures/.
+# Regenerate a figure (and this entry) only when a file in its `depends`
+# list changes the numbers it shows — see docs/DOC_STANDARD.md.
+
+
+[ballooning.delta_prime_matching]
+commit = "legacy"
+date = "unknown"
+depends = ["src/ForceFreeStates/Ballooning.jl"]
+note = "Predates the provenance system; compared methods for extracting the ballooning Delta_prime and c_a1. Original generator script not preserved."
+script = "(not preserved)"
+
+[ballooning.mercier_comparison]
+commit = "legacy"
+date = "unknown"
+depends = ["src/ForceFreeStates/Ballooning.jl"]
+note = "Predates the provenance system; compared the since-removed standalone Mercier.jl against Bal.jl. Not reproducible from current code (Mercier.jl was removed once Bal.jl subsumed it)."
+script = "(not preserved)"
+
+[inner_layer.backend_regime_map]
+commit = "3a0837e3-dirty"
+date = "2026-07-08"
+depends = ["src/InnerLayer/GGJ/Ray.jl", "src/InnerLayer/GGJ/Galerkin.jl", "src/InnerLayer/GGJ/Shooting.jl"]
+script = "make_backend_regime_map.jl"
+
+[inner_layer.convergence_Sinvariance]
+commit = "3a0837e3-dirty"
+date = "2026-07-08"
+depends = ["src/InnerLayer/GGJ/Ray.jl"]
+script = "make_convergence_Sinvariance.jl"
+
+[inner_layer.rotated_ray_contour]
+commit = "3a0837e3-dirty"
+date = "2026-07-08"
+depends = ["src/InnerLayer/GGJ/Ray.jl"]
+script = "make_rotated_ray_contour.jl"
+
+[inner_layer.solution_profiles]
+commit = "3a0837e3-dirty"
+date = "2026-07-08"
+depends = ["src/InnerLayer/GGJ/Ray.jl", "src/InnerLayer/GGJ/RayAsymptotics.jl", "src/InnerLayer/GGJ/InnerAsymptotics.jl"]
+script = "make_solution_profiles.jl"
diff --git a/docs/src/inner_layer.md b/docs/src/inner_layer.md
index cc80f52e..2d9a2246 100644
--- a/docs/src/inner_layer.md
+++ b/docs/src/inner_layer.md
@@ -1,17 +1,272 @@
 # Inner Layer Module
 
-The InnerLayer module provides abstract scaffolding for resistive inner-layer
-models used in matched asymptotic expansions for resistive MHD stability
-analysis. It currently includes the GGJ (Glasser–Greene–Johnson) shooting
-method for computing the inner-layer response.
+The `InnerLayer` module solves the **resistive inner-region problem** of
+matched-asymptotic resistive-MHD stability: the thin layer around a rational
+surface where resistivity, inertia and the ideal singularity balance. Its
+output is the pair of parity-projected matching data ``(\Delta_\mathrm{odd},
+\Delta_\mathrm{even})`` that the outer ideal region (DCON/GPEC) matches to in
+order to obtain resistive growth rates, tearing ``\Delta'``, and resistive
+interchange stability.
 
-## InnerLayer
+The module provides an abstract [`InnerLayerModel`](@ref InnerLayer.InnerLayerModel) interface and one
+concrete model so far, the **GGJ (Glasser–Greene–Johnson)** layer, with three
+interchangeable solver backends. Future inner-layer models (SLAYER, kinetic
+layers) plug in through the same interface.
+
+Two source papers define the equations and the asymptotic construction, and are
+cited by their equation numbers throughout the code and below:
+
+- **GWP2016** — A. H. Glasser, Z. R. Wang and J.-K. Park, *Computation of
+  resistive instabilities by matched asymptotic expansions*, Phys. Plasmas
+  **23**, 112506 (2016).
+- **GW2020** — A. H. Glasser and Z. R. Wang, *Asymptotic solutions and
+  convergence studies of the resistive inner region equations*, Phys. Plasmas
+  **27**, 012506 (2020).
+
+## Governing equations
+
+At a rational surface the layer equations couple three fields — the perturbed
+flux ``\Psi`` and two auxiliary layer variables ``\Xi`` and ``\Upsilon`` — as
+functions of the scaled distance ``x`` from the surface. In the GWP2016 form
+(Eq. 11 ≡ GW2020 Eq. 1) they are
+
+```math
+\begin{aligned}
+\Psi''    &= H\,\Upsilon' + Q\,(\Psi - x\,\Xi),\\[2pt]
+\Xi''     &= \frac{1}{Q^{2}}\Big[\,Q x^{2}\,\Xi - Q x\,\Psi - (E+F)\,\Upsilon - H\,\Psi'\,\Big],\\[2pt]
+\Upsilon''&= \frac{1}{Q}\Big[\,x^{2}\,\Upsilon - x\,\Psi - Q^{2}\big(G(\Xi-\Upsilon) - K(E\Xi + F\Upsilon + H\Psi')\big)\Big],
+\end{aligned}
+```
+
+where ``E, F, G, H, K`` are the flux-surface-averaged equilibrium coefficients
+of GWP2016 Eq. (A8) and ``Q`` is the dimensionless growth rate defined below.
+These coefficients are collected in [`GGJParameters`](@ref InnerLayer.GGJParameters).
+
+### Dimensionless scales
+
+The layer is scaled by the local Lundquist number ``S = \tau_R/\tau_A`` (ratio
+of resistive to Alfvén time). The displacement and growth-rate scales are
+
+```math
+X_0 = S^{-1/3}, \qquad Q_0 = X_0/\tau_A \qquad \text{(GWP2016 Eqs. A14–A15)},
+```
+
+so a physical growth rate ``\gamma`` maps to the scaled eigenvalue
+``Q = \gamma/Q_0`` ([`inner_Q`](@ref InnerLayer.inner_Q)) that appears in the equations above, and
+the scaled matching data are returned to physical units by the
+``X_0^{-2\sqrt{-D_I}} = S^{\,2\sqrt{-D_I}/3}`` rescaling (together with the
+``v_1^{\,2\sqrt{-D_I}}`` linear-scale factor; [`rescale_delta`](@ref InnerLayer.rescale_delta)).
+
+### Mercier index and matching exponents
+
+The premise of an inner-layer model is local Mercier stability. The Mercier
+interchange index and its resistive counterpart are
+
+```math
+D_I = E + F + H - \tfrac14, \qquad
+D_R = E + F + H^2 = D_I + \big(H-\tfrac12\big)^2 \qquad \text{(GWP2016 A9–A10)},
+```
+
+([`mercier_di`](@ref InnerLayer.mercier_di), [`mercier_dr`](@ref InnerLayer.mercier_dr)). For ``D_I < 0`` the large-``x``
+solutions are the two power laws ``x^{r_\pm}`` with the Frobenius exponents
+
+```math
+r_\pm = \tfrac32 \pm \sqrt{-D_I}, \qquad p_1 \equiv \sqrt{-D_I}
+\qquad \text{(GW2020 Eq. 49)}.
+```
+
+The matching datum ``\Delta`` is the amplitude of the large solution
+``x^{r_+}`` when the small solution ``x^{r_-}`` is normalized to unity (i.e. the
+ratio of large to small coefficient), in the physical (outer-matching)
+normalization. Imposing the two parities at ``x=0`` (odd:
+``\Psi'(0)=\Xi(0)=\Upsilon(0)=0``; even: ``\Psi(0)=\Xi'(0)=\Upsilon'(0)=0``)
+gives the pair ``(\Delta_\mathrm{odd}, \Delta_\mathrm{even})`` of GWP2016
+Eqs. (34)–(35).
+
+## Numerical method
+
+!!! note "Benchmark surface used in the figures"
+    Every figure on this page is computed on the ``q = 4`` rational-surface
+    benchmark ([`q4_surface_benchmark`](@ref InnerLayer.q4_surface_benchmark)) —
+    a fixed set of GGJ layer coefficients (``S = \tau_R/\tau_A \approx 4.6\times10^{6}``,
+    ``D_I \approx -0.31``) representative of a physical ``q = 4`` surface. It is a
+    stand-alone inner-layer coefficient set, not extracted from a particular
+    equilibrium run such as the DIII-D-like example. The well-conditioned
+    real-``Q`` cross-check in the Validation section instead uses the
+    Glasser & Wang (2020) Eq. (55) surface.
+
+### Solver backends
+
+The `GGJ` model exposes three solvers through the `solver` type parameter of
+[`GGJModel`](@ref InnerLayer.GGJModel); all consume the same large-``x`` asymptotic basis and return
+``\Delta`` in the same convention:
+
+| backend | method | robust regime |
+|:--------|:-------|:--------------|
+| `:shooting` | backward stable shoot from ``X_\mathrm{max}\to 0`` | ``\lvert Q\rvert \ll 1`` |
+| `:galerkin` | Hermite-cubic finite-element (real axis) | ``\lvert Q\rvert \lesssim 1`` |
+| `:ray` (default) | rotated-contour spectral-element collocation | ``\lvert Q\rvert \sim 500`` on/near the imaginary axis |
+
+The difficulty is that ``\Delta(Q)`` has poles on the imaginary-``Q`` axis, and
+both older backends lose accuracy off the real axis: `:shooting` through the
+exponential dichotomy of the backward shoot, `:galerkin` through real-axis
+oscillation and an on-axis pseudo-resonance. Measured against the `:ray`
+reference along the imaginary axis, `:shooting` holds to ``|Q|\sim 1`` and
+`:galerkin` to ``|Q|\sim 4`` before both lose all accuracy:
+
+![Relative error of the :shooting and :galerkin backends against the :ray reference along the imaginary-Q axis, on the q=4 benchmark surface. Each curve's crossing of the 1% line marks that method's practical reach.](figures/inner_layer/backend_regime_map.png)
+
+The `:ray` backend was written to reach the large-``|Q|`` imaginary-axis regime
+(rotation and resistivity scans) where these fail. The remainder of this section
+describes it.
+
+### Entire-solution formulation
+
+`:ray` works with the **plain** state ``v = (\Psi, \Xi, \Upsilon, \Psi',
+\Xi', \Upsilon')`` and writes the layer equations as the first-order system
+
+```math
+\frac{dv}{dx} = M(x)\,v .
+```
+
+The coefficient matrix (the `GGJ` internal `ode_matrix`) is **polynomial** in ``x``, so
+``x = 0`` is an ordinary point — the ``x^{-2}/x^{-4}`` singularities of the
+GW2020 Eq. (2) scaled form ``(x\Psi,\ \Psi'/x,\ \dots)`` are artifacts of that
+scaling, not of the equations. Because the coefficients are entire, the system
+continues analytically to complex ``x``, which is what makes the contour
+rotation below legitimate.
+
+### The rotated ray
+
+The equations are continued onto the ray
+
+```math
+x = e^{i\theta}\, s, \qquad s \in [0, S], \qquad \theta = \tfrac14\arg Q,
+```
+
+The angle ``\theta = \tfrac14\arg Q`` makes the parabolic-cylinder exponent of
+the outer solutions exactly real and clears the pseudo-resonance at
+``x^2 = -Q^2(G + K F)``, which on the imaginary-``Q`` axis is real and large and
+therefore sits directly on the un-rotated (real-``x``) contour. Rotating the
+contour lifts it off that point; ``\theta = 0`` recovers a real-axis solve.
+
+![The rotated integration ray in the complex layer-coordinate plane at Q = 500i. The real-axis contour runs through the pseudo-resonance x² = −Q²(G+KF); the ray at θ = arg(Q)/4 = 22.5° clears it.](figures/inner_layer/rotated_ray_contour.png)
+
+### Spectral-element collocation BVP
+
+On ``[0, s_m]`` the system is discretized by a **global Chebyshev
+spectral-element collocation** boundary-value problem: right-biased (Radau-like)
+collocation at the Chebyshev–Lobatto nodes of each cell, three parity conditions
+at the ordinary-point origin ``s = 0``, and six matching conditions at the outer
+edge,
+
+```math
+v(S) - \Delta\, U_b - c_1 E_1 - c_2 E_2 = U_s ,
+```
+
+with the matching datum ``\Delta`` and the two decaying-mode amplitudes
+``c_1, c_2`` carried as **bordered unknowns**. Here ``U_s, U_b`` are the small
+and large power-like solutions and ``E_{1,2}`` the forward-decaying exponential
+pair. No quantity is ever propagated across the layer, so the exponential
+dichotomy that limits the shooting backend never enters; the boundary condition
+splits it exactly.
+
+The two parities differ only in their three ``s=0`` rows, i.e. by a rank-3
+update, so both are obtained from a **single sparse LU factorization** plus a
+Woodbury correction (the `GGJ` internal `_solve_parities`) rather than two
+factorizations. A residual-driven bisection refinement adds cells until the
+collocation residual meets tolerance, with a roundoff-plateau guard.
+
+### Far-field boundary and the inward march
+
+The large-``x`` boundary data use the **same** `inps` Wasow asymptotic basis as
+the other backends (GW2020 Eqs. 3–52), evaluated at complex ``x`` by
+`RayAsymptotics.jl` and applied at the series radius ``S`` where the series is
+trusted (residual below tolerance; [`pick_smax`](@ref InnerLayer.pick_smax)). For large ``|Q|`` the
+trusted radius ``S`` can be far outside the collocation domain ``s_m``, so the
+power-pair data are transported inward from ``S`` to ``s_m`` by an **L-stable
+2-stage Radau IIA march** in the quotient modulo the decaying exponential pair —
+the subspace in which ``\Delta`` is defined (the `GGJ` internal
+`march_boundary`). An L-stable implicit method is essential: the decaying pair
+grows under backward integration, so any explicit marcher is stability-limited
+to ``O(\rho S^2)`` steps, while the Radau march *damps* the unresolvable
+backward-growing directions instead of amplifying them.
+
+The damped-zone march runs in `Complex{Double64}` extended precision. At large
+``S`` the near-parallel power-pair geometry amplifies the structured backward
+error of the ill-conditioned implicit solves into ``\Delta``-mixing of order
+``10^{-4}`` at ``|Q| = 500`` in `Float64`; extended precision removes this
+floor, while the well-conditioned resolved band stays in `Float64`.
+
+The result is a seamless numeric↔asymptotic solution: the collocation solution
+on ``[0, s_m]`` and the analytic ``u_\mathrm{small} + \Delta\,u_\mathrm{big}``
+representation for ``s \ge S`` share the same power-law tail — the overlap the
+outer-region matching relies on.
+
+![Inner-layer fields Ψ, Ξ, Υ on the rotated ray for the q=4 surface at Q = 2i. The collocation solution (solid) joins the asymptotic representation (dashed) seamlessly at the match point S = s_m.](figures/inner_layer/solution_profiles.png)
+
+## Validation and benchmarks
+
+**Cross-check against the Fortran `rmatch` solver.** At the Glasser & Wang
+(2020) Eq. (55) operating point ``Q = 0.1234`` (real, well-conditioned) the
+Julia GGJ solvers and the Fortran `rmatch deltac` solver agree to ``\sim
+10^{-8}``:
+
+| quantity | value (`:galerkin` and `:ray`) |
+|:---------|:-------------------------------|
+| ``\Delta_\mathrm{odd}`` | ``3.698368\times 10^{4}`` |
+| ``\Delta_\mathrm{even}`` | ``14.759721`` |
+
+**Physical benchmark on the imaginary axis.** On the ``q=4`` rational-surface
+benchmark ([`q4_surface_benchmark`](@ref InnerLayer.q4_surface_benchmark), ``S \approx 4.58\times10^6``,
+``D_I \approx -0.312``) the `:ray` backend is pinned at ``Q = 500i``, a regime
+entirely beyond `:galerkin`:
+
+| quantity | value at ``Q = 500i`` |
+|:---------|:----------------------|
+| ``\Delta_\mathrm{odd}`` | ``2.47206 + 13.35412\,i`` |
+| ``\Delta_\mathrm{even}`` | ``0.137497 + 0.742755\,i`` |
+
+**Convergence and contour invariance.** Because ``\Delta`` is an analytic
+invariant of the contour angle, re-solving with each numerical knob perturbed on
+an independent axis ([`delta_convergence`](@ref InnerLayer.delta_convergence)) gives an honest error bar: at
+``Q = 500i`` the worst-case spread across all knobs is ``\sim 5\times10^{-6}``
+for ``\Delta_\mathrm{odd}`` and ``\sim 6\times10^{-7}`` for
+``\Delta_\mathrm{even}``.
+
+![Relative change of Δ at Q = 500i under independent perturbations of each numerical knob (contour angle, spectral order, series order/radius, refinement depth, march tolerance, handoff radius, purification depth). The worst-case spread is the reported error bar.](figures/inner_layer/convergence_Sinvariance.png)
+
+## Choosing a backend
+
+`:ray` is the **default** — `GGJModel()` constructs `GGJModel{:ray}()`. It is
+the correct choice for ``|Q| \gtrsim 1`` and for any ``Q`` near the imaginary
+axis. `:galerkin` remains available and may be faster for very small real
+``|Q|``, but note the backends take different numerical-knob keywords: pass
+`GGJModel(solver=:galerkin)` explicitly when supplying the Galerkin `nx`/`xfac`
+knobs, since those are silently ignored by `:ray`.
+
+```julia
+using GeneralizedPerturbedEquilibrium.InnerLayer
+
+p = q4_surface_benchmark()          # GGJParameters for the q=4 benchmark surface
+γ = 500im * GGJ.q0(p)               # physical rate placing Q on the imaginary axis at 500i
+
+Δ = solve_inner(GGJModel(), p, γ)   # (Δ_odd, Δ_even) with the default :ray backend
+
+# Full result (raw Δ, contour, mesh, nodal solution, residuals) via solve_ray:
+res = solve_ray(p, GGJ.inner_Q(p, γ))
+res.Δ, res.resid, res.bc_cond
+```
+
+## API Reference
+
+### InnerLayer
 
 ```@autodocs
 Modules = [GeneralizedPerturbedEquilibrium.InnerLayer]
 ```
 
-## GGJ
+### GGJ
 
 ```@autodocs
 Modules = [GeneralizedPerturbedEquilibrium.InnerLayer.GGJ]
diff --git a/regression-harness/cases/ggj_ray_q500i.toml b/regression-harness/cases/ggj_ray_q500i.toml
new file mode 100644
index 00000000..967fe379
--- /dev/null
+++ b/regression-harness/cases/ggj_ray_q500i.toml
@@ -0,0 +1,53 @@
+# Regression case: GGJ inner-layer rotated-ray collocation backend at Q = 500i.
+# A "computed" case (kind = "computed", no example_dir) that calls solve_inner with the
+# :ray backend on the physical q=4 benchmark surface at Q = 500i — a regime beyond the
+# :galerkin backend — and tracks the Δ_odd / Δ_even matching coefficients. Guards the
+# rotated-ray solver (and its extended-precision damped-zone march) against numerical
+# drift across changes. Each [quantities.*] block names an HDF5 path in the computed
+# output, how to extract it, and the noise floor for diffs.
+[case]
+name = "ggj_ray_q500i"
+description = "GGJ inner-layer :ray backend, q=4 surface, Q = 500i (beyond :galerkin)"
+kind = "computed"
+
+# Δ_odd / Δ_even from solve_inner(GGJModel(:ray), q4_surface_benchmark(), 500im·Q₀).
+# Complex-valued; both real and imaginary parts are tracked separately.
+[quantities.delta_odd_real]
+h5path = "ggj/delta_odd_real"
+type = "real_scalar"
+extract = "value"
+label = "Re(Δ_odd)"
+noise_threshold = 1e-8
+order = 10
+
+[quantities.delta_odd_imag]
+h5path = "ggj/delta_odd_imag"
+type = "real_scalar"
+extract = "value"
+label = "Im(Δ_odd)"
+noise_threshold = 1e-8
+order = 11
+
+[quantities.delta_even_real]
+h5path = "ggj/delta_even_real"
+type = "real_scalar"
+extract = "value"
+label = "Re(Δ_even)"
+noise_threshold = 1e-8
+order = 12
+
+[quantities.delta_even_imag]
+h5path = "ggj/delta_even_imag"
+type = "real_scalar"
+extract = "value"
+label = "Im(Δ_even)"
+noise_threshold = 1e-8
+order = 13
+
+[quantities.runtime]
+h5path = ""
+type = "runtime"
+extract = "value"
+label = "Runtime (s)"
+noise_threshold = 0.0
+order = 90
diff --git a/regression-harness/src/runner.jl b/regression-harness/src/runner.jl
index 8705cb4d..84e80838 100644
--- a/regression-harness/src/runner.jl
+++ b/regression-harness/src/runner.jl
@@ -79,6 +79,31 @@ open(ARGS[2], "w") do f
 end
 """
 
+# GGJ rotated-ray backend at Q = 500i on the q=4 benchmark surface — a regime beyond the
+# :galerkin backend. Builds the physical rate γ = 500i·Q₀ so inner_Q lands exactly on the
+# imaginary axis at 500i, then writes the parity matching data. Runtime to ARGS[2] as usual.
+const COMPUTED_GGJ_RAY_SCRIPT_TEMPLATE = """
+using Pkg
+%INSTANTIATE%
+using GeneralizedPerturbedEquilibrium
+using GeneralizedPerturbedEquilibrium.InnerLayer
+using HDF5
+p = q4_surface_benchmark()
+γ = 500.0im * InnerLayer.GGJ.q0(p)
+t_start = time()
+Δ = solve_inner(GGJModel(solver=:ray), p, γ)
+elapsed = time() - t_start
+h5open(ARGS[1], "w") do fid
+    fid["ggj/delta_odd_real"]  = real(Δ[1])
+    fid["ggj/delta_odd_imag"]  = imag(Δ[1])
+    fid["ggj/delta_even_real"] = real(Δ[2])
+    fid["ggj/delta_even_imag"] = imag(Δ[2])
+end
+open(ARGS[2], "w") do f
+    println(f, elapsed)
+end
+"""
+
 # Self-contained separatrix-finder regression (PR #296). Loads a fixed-boundary EFIT whose
 # computational box hugs the LCFS (eps=0.05 TokaMaker aspect-scan g-file): outside the prescribed
 # LCFS the coil-vacuum flux turns back above the boundary value before the grid edge, so the old
@@ -139,6 +164,8 @@ Pick the subprocess script template for a `kind="computed"` case.
 function _computed_script_template(case_spec::CaseSpec)
     if case_spec.name == "ggj_reference"
         return COMPUTED_GGJ_SCRIPT_TEMPLATE
+    elseif case_spec.name == "ggj_ray_q500i"
+        return COMPUTED_GGJ_RAY_SCRIPT_TEMPLATE
     elseif case_spec.name == "efit_fixedbdy_separatrix"
         return COMPUTED_SEPARATRIX_SCRIPT_TEMPLATE
     end
diff --git a/src/InnerLayer/GGJ/GGJ.jl b/src/InnerLayer/GGJ/GGJ.jl
index e7f0488a..289d3c38 100644
--- a/src/InnerLayer/GGJ/GGJ.jl
+++ b/src/InnerLayer/GGJ/GGJ.jl
@@ -1,15 +1,20 @@
 # GGJ.jl
 #
-# Glasser–Greene–Johnson resistive inner-layer model. Provides two
+# Glasser–Greene–Johnson resistive inner-layer model. Provides three
 # interchangeable solvers selected via the `solver` type-parameter of
 # `GGJModel`:
 #
-#   - `:shooting`  – stable backward shoot from X_max → 0 (Phase 3)
-#   - `:galerkin`  – Hermite-cubic finite element method (Phase 4)
+#   - `:shooting`  – stable backward shoot from X_max → 0;
+#                    numerically stable only for |Q| ≪ 1, should not be used
+#   - `:galerkin`  – Hermite-cubic finite element method;
+#                    real-axis method, degrades for |Q| ≳ 1 off the real axis
+#   - `:ray`       – rotated-contour spectral-element collocation (Ray.jl);
+#                    robust to |Q| ~ 500 on/near the imaginary axis
 #
-# Both solvers share the same `inps` Wasow asymptotic-basis kernel
-# (`InnerAsymptotics.jl`) for the large-x boundary condition. They return
-# the parity-projected matching data `(Δ_odd, Δ_even)` of GWP2016 Eqs. (34)–(35).
+# All solvers share the same `inps` Wasow asymptotic-basis kernel
+# (`InnerAsymptotics.jl`, complex-x extension in `RayAsymptotics.jl`) for the
+# large-x boundary condition. They return the parity-projected matching data
+# `(Δ_odd, Δ_even)` of GWP2016 Eqs. (34)–(35) in the same (deltac) convention.
 #
 # Equation references throughout this module use two source papers:
 #
@@ -31,6 +36,10 @@ module GGJ
 
 using LinearAlgebra
 using StaticArrays
+using SparseArrays
+using Random
+using Printf
+using DoubleFloats: Double64
 
 import ..InnerLayerModel, ..solve_inner
 
@@ -38,25 +47,35 @@ import ..InnerLayerModel, ..solve_inner
     GGJModel{S} <: InnerLayerModel
 
 Glasser–Greene–Johnson resistive inner-layer model. The type parameter `S`
-selects the solver: `:galerkin` (default) for the Hermite-cubic finite element
-solver and `:shooting` for the backward stable-shoot solver. Both
-implementations consume the same `inps` asymptotic-basis kernel and return
-the parity-projected matching data.
+selects the solver: `:ray` (default) for the rotated-contour collocation
+solver (robust at large |Q| on/near the imaginary axis; agrees with
+`:galerkin` within its error bar at moderate Q), `:galerkin` for the
+Hermite-cubic finite element solver (real-axis method; degrades for
+|Q| ≳ 1), and `:shooting` for the backward stable-shoot solver (|Q| ≪ 1
+only). All implementations consume the same `inps` asymptotic-basis kernel
+and return the parity-projected matching data in the same convention.
+Note the backends take different numerical-knob keywords: pass
+`GGJModel(solver=:galerkin)` explicitly when using `nx`/`xfac`-style
+Galerkin knobs.
 """
 struct GGJModel{S} <: InnerLayerModel end
 
-GGJModel(; solver::Symbol=:galerkin) = GGJModel{solver}()
+GGJModel(; solver::Symbol=:ray) = GGJModel{solver}()
 
 include("GGJParameters.jl")
 include("InnerAsymptotics.jl")
+include("RayAsymptotics.jl")
 include("Reference.jl")
 include("Shooting.jl")
 include("Galerkin.jl")
+include("Ray.jl")
 
 export GGJModel, GGJParameters
 export mercier_di, mercier_dr, inner_Q, rescale_delta
 export build_asymptotics, evaluate_asymptotics, pick_xmax
 export InnerAsymptoticsCache
 export glasser_wang_2020_eq55
+export solve_ray, RaySolveResult, pick_smax, physical_ua_dua
+export delta_convergence, solution_profile, asymptotic_profile, q4_surface_benchmark
 
 end # module GGJ
diff --git a/src/InnerLayer/GGJ/GGJParameters.jl b/src/InnerLayer/GGJ/GGJParameters.jl
index 732ab781..8ab243e8 100644
--- a/src/InnerLayer/GGJ/GGJParameters.jl
+++ b/src/InnerLayer/GGJ/GGJParameters.jl
@@ -112,8 +112,10 @@ inner_Q(p::GGJParameters, γ::Number) = ComplexF64(γ) / q0(p)
     rescale_delta(Δ, p::GGJParameters) -> SVector{2,ComplexF64}
 
 Map the scaled inner-layer matching data back to physical Δ at the rational
-surface by the `X₀^(2√(−D_I))` rescaling implied by the power-like matching
-(GWP2016 Sec. IV; the inner solution ~ X^{r±} converts to physical x = X₀ X).
+surface by the `X₀^(−2√(−D_I)) = S^(2√(−D_I)/3)` rescaling, together with the
+`v₁^(2√(−D_I))` linear-scale factor, implied by the power-like matching
+(GWP2016 Sec. IV; the inner solution ~ X^{r±} converts to physical x = X₀ X, so
+the large/small amplitude ratio Δ scales by X₀^{−(r₊−r₋)} = X₀^(−2√(−D_I))).
 Operates element-wise on a 2-vector of `(Δ_odd, Δ_even)`.
 """
 function rescale_delta(Δ::AbstractVector, p::GGJParameters)
diff --git a/src/InnerLayer/GGJ/Ray.jl b/src/InnerLayer/GGJ/Ray.jl
new file mode 100644
index 00000000..174aebf8
--- /dev/null
+++ b/src/InnerLayer/GGJ/Ray.jl
@@ -0,0 +1,1183 @@
+# Ray.jl
+#
+# Rotated-contour spectral-element collocation backend (`GGJModel{:ray}`)
+# for the GGJ inner-layer matching data — robust at large |Q| on and near
+# the imaginary axis, where both `:shooting` (exponential dichotomy) and
+# `:galerkin` (real-axis oscillation + on-axis pseudo-resonance) degrade.
+#
+# Ported 2026-07-08 from the standalone GGJRay package (~/Projects/GGJ_test),
+# validated there against the Fortran rmatch pins, the `:galerkin` backend at
+# moderate Q, and physical benchmarks to Q = 500i. Complete methods paper:
+# GGJ_test/docs/METHOD.md. In one paragraph:
+#
+#   The equations are continued analytically onto the ray x = e^{iθ}s with
+#   θ = arg(Q)/4 (parabolic-cylinder exponent exactly real; pseudo-resonance
+#   cleared). On [0, s_m] a global Chebyshev spectral-element collocation BVP
+#   in the plain variables v = (Ψ, Ξ, Υ, Ψ', Ξ', Υ') is solved with parity
+#   conditions at the (ordinary-point) origin and Δ, c₁, c₂ as bordered
+#   unknowns of one sparse LU per parity (rank-3 Woodbury for the second).
+#   The far-field condition uses the SAME inps Wasow basis as the other
+#   backends (RayAsymptotics.jl evaluates it at complex x; the rdcon
+#   normalization of Δ is retained), applied at the series radius S and
+#   transported S → s_m by an L-stable 2-stage Radau IIA march in the
+#   quotient modulo the decaying exponential pair. The damped-zone march
+#   arithmetic runs in Complex{Double64}: at large S the near-parallel
+#   power-pair geometry amplifies the structured backward error of the
+#   ill-conditioned implicit solves (Σ eps·z ≈ eps·ρS²/2) into Δ-mixing
+#   ~1e-4 at |Q| = 500 in Float64 (METHOD.md §8).
+#
+# Sections below mirror the GGJRay source layout:
+#   1. system   — plain-variable first-order ODE matrix, parity sets
+#   2. mesh     — Chebyshev cells, barycentric interpolation, WKB grading
+#   3. boundary — decaying pair, Radau IIA quotient march
+#   4. solve    — assembly, bordered solve, refinement, driver, diagnostics
+
+"""
+    q4_surface_benchmark() -> GGJParameters
+
+Physical q = 4 rational-surface benchmark point (S = τ_R/τ_A ≈ 4.58×10⁶,
+D_I ≈ −0.31166, α = √(−D_I) ≈ 0.5583). Primary validation point for the
+rotated-ray backend on the imaginary-Q axis (pinned at Q = 100i, 500i in
+the test suite).
+"""
+function q4_surface_benchmark()
+    return GGJParameters(;
+        E=-0.13733, F=0.022202, G=7.60633, H=0.053468, K=14.66987,
+        M=30.26883, taua=2.11226e-7, taur=0.968219, v1=1.55009
+    )
+end
+
+# system.jl
+#
+# The GGJ inner-region equations (GWP2016 Eq. 11 ≡ GW2020 Eq. 1) as a plain
+# first-order system in v = (Ψ, Ξ, Υ, Ψ', Ξ', Υ'), dv/dx = M(x)·v.
+# All coefficients are POLYNOMIAL in x, so x = 0 is an ordinary point (the
+# x⁻²/x⁻⁴ singularities of the GW2020 Eq. 2 form are artifacts of the
+# (xΨ, Ψ'/x) scaling) and the system continues analytically to complex x.
+#
+#   Ψ'' = H Υ' + Q(Ψ − xΞ)
+#   Ξ'' = [Q x² Ξ − Q x Ψ − (E+F) Υ − H Ψ'] / Q²
+#   Υ'' = [x² Υ − x Ψ − Q²(G(Ξ−Υ) − K(EΞ + FΥ + HΨ'))] / Q
+#
+# Row-by-row this matches the deltac `inpso_get_uv` matrices (I, V, U) after
+# dividing each equation by its diagonal second-derivative coefficient
+# (1, Q², Q).
+
+"""
+    ode_matrix([CT,] p::GGJParameters, Q, x) -> SMatrix{6,6,CT}
+
+Coefficient matrix `M(x)` of `dv/dx = M v`, `v = (Ψ, Ξ, Υ, Ψ', Ξ', Υ')`,
+valid for real or complex `x`. The optional leading type argument selects
+the element type (e.g. `Complex{Double64}` for the extended-precision
+damped-zone march); default `ComplexF64`.
+"""
+ode_matrix(p::GGJParameters, Q::ComplexF64, x::Number) =
+    ode_matrix(ComplexF64, p, Q, x)
+
+function ode_matrix(::Type{CT}, p::GGJParameters, Q::Number, x::Number) where {CT<:Complex}
+    e = CT(p.E)
+    f = CT(p.F)
+    g = CT(p.G)
+    h = CT(p.H)
+    k = CT(p.K)
+    q = CT(Q)
+    q2 = q * q
+    xc = CT(x)
+    x2 = xc * xc
+
+    m = zeros(MMatrix{6,6,CT})
+    # rows 1–3: (Ψ, Ξ, Υ)' = (Ψ', Ξ', Υ')
+    m[1, 4] = 1
+    m[2, 5] = 1
+    m[3, 6] = 1
+    # row 4: Ψ'' = Q Ψ − Q x Ξ + H Υ'
+    m[4, 1] = q
+    m[4, 2] = -q * xc
+    m[4, 6] = h
+    # row 5: Ξ'' = −(x/Q) Ψ + (x²/Q) Ξ − ((E+F)/Q²) Υ − (H/Q²) Ψ'
+    m[5, 1] = -xc / q
+    m[5, 2] = x2 / q
+    m[5, 3] = -(e + f) / q2
+    m[5, 4] = -h / q2
+    # row 6: Υ'' = −(x/Q) Ψ − Q(G−KE) Ξ + (x²/Q + Q(G+KF)) Υ + HKQ Ψ'
+    m[6, 1] = -xc / q
+    m[6, 2] = -q * (g - k * e)
+    m[6, 3] = x2 / q + q * (g + k * f)
+    m[6, 4] = h * k * q
+
+    return SMatrix{6,6,CT}(m)
+end
+
+"""
+    parity_rows(isol) -> SVector{3,Int}
+
+Component indices of `v = (Ψ, Ξ, Υ, Ψ', Ξ', Υ')` that vanish at s = 0 for
+each parity solution, mirroring deltac_set_boundary:
+
+  - `isol = 1` ("odd"):  Ψ'(0) = 0, Ξ(0) = 0, Υ(0) = 0  → components (4, 2, 3)
+  - `isol = 2` ("even"): Ψ(0) = 0, Ξ'(0) = 0, Υ'(0) = 0 → components (1, 5, 6)
+"""
+parity_rows(isol::Int) = isol == 1 ? SVector(4, 2, 3) : SVector(1, 5, 6)
+# mesh.jl
+#
+# Graded spectral-element mesh on the ray parameter s ∈ [0, S], plus the
+# Chebyshev–Lobatto reference-cell machinery (nodes, differentiation matrix,
+# barycentric interpolation).
+#
+# The initial grading tracks the local stiffness scales of the rotated
+# system:
+#   - the fast Ohm's-law pair e^{±√Q x}: active near the origin, decaying
+#     at rate Re(√Q e^{iθ}) — resolved on a fine zone [0, s₁];
+#   - the parabolic-cylinder pair e^{±ρ s²/2}, ρ = Re(e^{2iθ}/√Q): decaying
+#     content is above relative machine floor only for s ≲ s_d = √(2D/ρ) —
+#     resolved on a rate-adapted zone [s₁, s_d];
+#   - beyond s_d the physical solution is a smooth power law: geometric
+#     (log-spaced) cells suffice out to S.
+# Residual-driven bisection refinement (solve.jl) then repairs whatever the
+# initial guess missed.
+
+"""
+    cheblobatto(p) -> (t, D)
+
+Chebyshev–Lobatto nodes on [-1, 1] in ASCENDING order and the spectral
+differentiation matrix `D` for those nodes (Trefethen's `cheb`, reflected).
+"""
+function cheblobatto(p::Int)
+    p >= 1 || throw(ArgumentError("polynomial order p must be ≥ 1"))
+    # Standard descending Trefethen nodes x_j = cos(jπ/p), j = 0..p.
+    x = [cos(pi * j / p) for j in 0:p]
+    c = [j == 0 || j == p ? 2.0 : 1.0 for j in 0:p] .* [(-1.0)^j for j in 0:p]
+    D = zeros(p + 1, p + 1)
+    for i in 1:(p+1), j in 1:(p+1)
+        if i != j
+            D[i, j] = (c[i] / c[j]) / (x[i] - x[j])
+        end
+    end
+    for i in 1:(p+1)
+        D[i, i] = -sum(D[i, j] for j in 1:(p+1) if j != i)
+    end
+    # Reflect to ascending order: y = -x is ascending with the SAME index
+    # order (x is descending), and for g(y) := f(-y) sampled on y the chain
+    # rule gives D_y = -D. No index permutation.
+    t = -x
+    Dt = -D
+    return t, Dt
+end
+
+"""
+    barycentric_weights(t) -> w
+
+Barycentric interpolation weights for nodes `t` (any distribution).
+"""
+function barycentric_weights(t::AbstractVector{Float64})
+    n = length(t)
+    w = ones(Float64, n)
+    for j in 1:n
+        for k in 1:n
+            k == j && continue
+            w[j] /= (t[j] - t[k])
+        end
+    end
+    return w
+end
+
+"""
+    barycentric_eval(t, w, vals, τ) -> value
+
+Evaluate the interpolant of `vals` (matrix: nodes × components) at `τ`.
+Returns a vector over components.
+"""
+function barycentric_eval(t::AbstractVector{Float64}, w::AbstractVector{Float64},
+    vals::AbstractMatrix{ComplexF64}, τ::Float64)
+    n = length(t)
+    for j in 1:n
+        if τ == t[j]
+            return vals[j, :]
+        end
+    end
+    num = zeros(ComplexF64, size(vals, 2))
+    den = 0.0
+    for j in 1:n
+        fac = w[j] / (τ - t[j])
+        den += fac
+        @views num .+= fac .* vals[j, :]
+    end
+    return num ./ den
+end
+
+"""
+    initial_breaks(params, Q, θ, S, p; cfl=0.4, drop=40.0, ratio=1.3,
+                   hmax_frac=0.08) -> Vector{Float64}
+
+Build the initial cell-boundary vector `0 = s₀ < s₁ < … = S` using the
+three-zone grading described in the header. `cfl` ≈ resolved-exponential
+rate per node; `drop` = e-folds after which decaying exponential content is
+considered machine-dead; `ratio` = geometric growth in the outer zone.
+"""
+function initial_breaks(params::GGJParameters, Q::ComplexF64, θ::Float64,
+    S::Float64, p::Int; cfl::Float64=0.4, drop::Float64=40.0,
+    ratio::Float64=1.3, hmax_frac::Float64=0.08)
+    sq = sqrt(Q)
+    ρ = real(cis(2θ) / sq)                 # parabolic pair: exponent ρ s²/2
+    ρ > 0 || throw(ArgumentError("contour angle θ=$θ gives non-decaying parabolic pair (ρ=$ρ ≤ 0); use θ ≈ arg(Q)/4"))
+    # Fast local WKB scales: the Ohm pair |√Q| and the Υ-family
+    # √|Q(G+KF)| (dominant on extreme-coefficient surfaces, K ~ 10⁶ —
+    # dormant almost everywhere but must be resolved; see solve.jl header).
+    μ = abs(sq) + sqrt(abs(Q) * abs(params.G + params.K * params.F))
+    μdec = max(real(sq * cis(θ)), 0.1 * μ) # its decay rate along the ray
+
+    s1 = min(drop / μdec, 0.7 * S)         # fast pair dead beyond s1
+    sd = min(sqrt(2 * drop / ρ), 0.8 * S)  # parabolic pair dead beyond sd
+    sd = max(sd, min(1.05 * s1, S))
+    hmax = hmax_frac * S
+
+    # Single grading rule: every cell resolves the full local WKB rate
+    # μ + ρs (fast pair + parabolic pair), whether the corresponding mode
+    # amplitudes are alive or dormant — unresolved dormant modes alias into
+    # spurious bounded discrete modes and poison the global solve. The
+    # outer power-law region is not discretized at all (handled by the
+    # boundary march), so this never gets expensive: the BVP domain ends at
+    # s_m ~ max(s_d, s₁) where μ dominates the rate budget.
+    _ = (s1, sd, ratio)  # zone markers retained in signature for tuning
+    breaks = [0.0]
+    while breaks[end] < S
+        s = breaks[end]
+        h = cfl * p / (μ + ρ * s)
+        h = min(h, hmax, S / 8)
+        snew = s + h
+        if snew >= S - 0.25 * h
+            snew = S
+        end
+        push!(breaks, min(snew, S))
+    end
+    # Guard against a degenerate final sliver.
+    if length(breaks) >= 3 && (breaks[end] - breaks[end-1]) < 0.2 * (breaks[end-1] - breaks[end-2])
+        deleteat!(breaks, length(breaks) - 1)
+    end
+    return breaks
+end
+# boundary.jl
+#
+# Large-s boundary data at the matching point on the ray x = e^{iθ}s.
+#
+# The admissible (bounded-as-s→∞) solution space is spanned by
+#   { u_small, u_big }  — the two power-like Mercier solutions, taken from
+#                         the inps series in the standard normalization, and
+#   { E₁, E₂ }          — the two forward-DECAYING exponential solutions.
+#
+# The inps series is only trusted at radius S (residual ≤ smax_tol), which
+# for extreme-coefficient surfaces (K ~ 10⁶) can be S ~ 10⁴–10⁵, while the
+# collocation BVP ends at s_m ~ 30–150. The gap is bridged by MARCHING the
+# power-pair data W = [u_small, u_big] backward S → s_m.
+#
+# Marching design (the hard-won part):
+#   - Backward, the decaying pair GROWS at rate ρs (ρ = Re(e^{2iθ}/√Q)), up
+#     to ~10³/unit at large s. Any explicit integrator — including
+#     projected/continuous-QR variants, whose unprojected stage points
+#     reintroduce the transverse stiffness — is stability-limited to
+#     h ≲ 3/(ρs), i.e. O(ρS²) steps: hopeless at S ~ 10⁵.
+#   - The system is LINEAR, so we use a stiffly-accurate L-STABLE implicit
+#     method (2-stage Radau IIA, order 3; one 12×12 solve per step). For
+#     resolved slow content it is 3rd-order accurate; for the unresolvable
+#     backward-growing exponential directions |R(z)| < 1 at large |z| —
+#     the method DAMPS what it cannot resolve instead of amplifying it.
+#     That is normally an accuracy vice; here it is exactly the required
+#     behavior, because Δ is defined modulo the decaying pair.
+#   - In the partially-resolved band (|z| = ρs·h ≲ 10) genuine backward
+#     growth of rounding-injected decaying-pair content still occurs, so W
+#     is PURGED against a locally-extracted decaying frame every
+#     Φ = ∫ρs ds ≈ 8 e-folds. Purging only ever removes span{E} content —
+#     legitimate by the definition of Δ — so the transported coefficients
+#     of u_small/u_big are exact up to integrator tolerance.
+
+"""
+    decaying_pair(params, Q, θ, S; growth=30.0, seed=20260707) -> E::Matrix{ComplexF64} (6×2)
+
+Orthonormal basis of the forward-decaying pair at s = S on the ray, via a
+SHORT backward RK4 integration over [S, S+ΔS], with ΔS chosen so the pair is
+amplified by ≈ e^growth relative to everything else (backward-attracting
+purification of random seeds). The integration is fully resolved
+(h·rate ≈ 0.05), so plain RK4 is stable over this short stretch.
+"""
+function decaying_pair(params::GGJParameters, Q::ComplexF64, θ::Float64, S::Float64;
+    growth::Float64=30.0, seed::Int=20260707, dp_res::Float64=20.0)
+    sq = sqrt(Q)
+    ρ = real(cis(2θ) / sq)
+    ρ > 0 || throw(ArgumentError("θ=$θ gives ρ=$ρ ≤ 0; decaying pair undefined"))
+    ph = cis(θ)
+    # Purification window: BOTH backward-growing directions must separate
+    # from the power pair by ≥ `growth` e-folds over [S, S+ΔS], so the
+    # budget is set by the SLOWEST grower's rate margin over the algebraic
+    # power rates (from the frozen spectrum) — not by ρS. At moderate |Q|
+    # the fast-Ohm grower can be ~0.6/unit while ρS ~ 10²: budgeting on ρS
+    # leaves the second frame column O(1)-contaminated with power content,
+    # which then poisons every march purge that uses the frame.
+    g2 = Inf
+    gmax_loc = 0.0
+    for λ in eigvals(Matrix(ph * ode_matrix(params, Q, ph * S)))
+        g = -real(λ)
+        g > 20.0 / S && (g2 = min(g2, g))
+        gmax_loc = max(gmax_loc, abs(real(λ)))
+    end
+    gsep = isfinite(g2) ? max(g2 - 20.0 / S, 0.05 * g2) : ρ * S
+    ΔS = min(growth / gsep, 3.0 * S)
+
+    # Step resolution must use the TRUE local spectrum, not the analytic
+    # ρs + |√Q| estimate: on extreme-coefficient surfaces the Υ-family rate
+    # √|Q(G+KF)| dominates at small s and the analytic estimate under-counts
+    # it 10×, leaving h·λ ~ 0.2 — stable but phase-inaccurate, which corrupts
+    # the extracted frame directions (observed as O(1) θ-drift downstream).
+    for λ in eigvals(Matrix(ph * ode_matrix(params, Q, ph * (S + ΔS))))
+        gmax_loc = max(gmax_loc, abs(real(λ)))
+    end
+    rate = max(ρ * (S + ΔS) + abs(sq), gmax_loc)
+    nsteps = clamp(ceil(Int, dp_res * rate * ΔS), 200, 2_000_000)
+    h = -ΔS / nsteps
+
+    rng = MersenneTwister(seed)
+    v = randn(rng, ComplexF64, 6, 2)
+    v = Matrix(qr(v).Q)
+
+    f(s, y) = (ph * ode_matrix(params, Q, ph * s)) * y
+
+    s = S + ΔS
+    for istep in 1:nsteps
+        k1 = f(s, v)
+        k2 = f(s + h / 2, v + (h / 2) * k1)
+        k3 = f(s + h / 2, v + (h / 2) * k2)
+        k4 = f(s + h, v + h * k3)
+        v = v + (h / 6) * (k1 + 2k2 + 2k3 + k4)
+        s += h
+        if istep % 50 == 0
+            v = Matrix(qr(v).Q)
+        end
+    end
+    return Matrix(qr(v).Q)
+end
+
+"""
+    march_boundary(params, Q, θ, S, s_m, Us, Ub; rtol=1e-9, growth=30.0, seed=20260707)
+        -> (Us_m, Ub_m, F)
+
+Transport the power-pair boundary data from the series radius `S` inward to
+`s_m` along the ray, modulo the decaying exponential pair (the quotient in
+which the matching data Δ is defined). Adaptive 2-stage Radau IIA (L-stable,
+stiffly accurate) on the raw linear system, with Φ-budgeted purges against
+locally-extracted decaying frames; see the file header for why explicit
+(including projected-QR) marchers cannot do this job at large S.
+"""
+function march_boundary(params::GGJParameters, Q::ComplexF64, θ::Float64,
+    S::Float64, s_m::Float64, Us::Vector{ComplexF64}, Ub::Vector{ComplexF64};
+    rtol::Float64=1e-9, growth::Float64=30.0, seed::Int=20260707,
+    ztrk::Float64=0.3, purge_budget::Float64=2.5, dp_res::Float64=20.0)
+    ph = cis(θ)
+    sq = sqrt(Q)
+    ρ = real(cis(2θ) / sq)
+    𝔐(s) = ph * ode_matrix(params, Q, ph * s)
+
+    W = hcat(copy(Us), copy(Ub))
+    function purge!(W, s)
+        Floc = decaying_pair(params, Q, θ, s; growth=growth, seed=seed, dp_res=dp_res)
+        W .-= Floc * (Floc' * W)
+        return Floc
+    end
+
+    I6 = Matrix{ComplexF64}(I, 6, 6)
+    # One Radau IIA step: A = [5/12 -1/12; 3/4 1/4], c = (1/3, 1), b = A[2,:].
+    function radau_step(s, W, h)
+        M1 = 𝔐(s + h / 3)
+        M2 = 𝔐(s + h)
+        Ablk = [I6-(5h/12)*M1 (h/12)*M1; -(3h/4)*M2 I6-(h/4)*M2]
+        rhs = [M1 * W; M2 * W]
+        K = Ablk \ rhs
+        return W + h * ((3 / 4) * @view(K[1:6, :]) + (1 / 4) * @view(K[7:12, :]))
+    end
+
+    # Extended-precision (Double64) Radau step for the damped-zone W march.
+    # At z = ρ·c·s² ≫ 1 the step matrices have cond ~ z and the LU backward
+    # error has a SYSTEMATIC direction set by the elimination structure — it
+    # does not dither away, accumulating a relative error Σ eps·z ≈ eps·ρS²/2
+    # on the power columns that is INDEPENDENT of the step law (c cancels).
+    # The near-parallel power-pair geometry at large S turns that into
+    # u_small/u_big coefficient mixing β ~ eps·ρS²/2 · (‖u_big‖/‖u_small‖)(S)
+    # ~ 1e-4 at |Q| = 500 — saturating Δ of the near-pole parity at ~1/β and
+    # flooring the other at |Δ|·β (the ε_mix defect; docs/BUGHUNT_epsmix.md).
+    # Double64 arithmetic (eps ~ 5e-33) removes it outright; the resolved
+    # band stays Float64 (z ≤ ztrk there, solves well-conditioned).
+    CTd = Complex{Double64}
+    phd = CTd(cos(Double64(θ)), sin(Double64(θ)))
+    I6d = Matrix{CTd}(I, 6, 6)
+    𝔐d(s) = phd * Matrix(ode_matrix(CTd, params, Q, phd * Double64(s)))
+    function radau_step_d(s, Wd, h)
+        hd = Double64(h)
+        M1 = 𝔐d(s + h / 3)
+        M2 = 𝔐d(s + h)
+        Ablk = [I6d-(5hd/12)*M1 (hd/12)*M1; -(3hd/4)*M2 I6d-(hd/4)*M2]
+        rhs = [M1 * Wd; M2 * Wd]
+        K = Ablk \ rhs
+        return Wd + hd * ((CTd(3) / 4) * K[1:6, :] + (CTd(1) / 4) * K[7:12, :])
+    end
+
+    # --- Schedule: geometric where damped, growth-capped where resolved. ---
+    # No adaptive controller (three generations died on estimator artifacts);
+    # the step law is derived from the local mode structure instead. At each
+    # step the eigenvalues of 𝔐(s) give the two backward-growing exponential
+    # rates g (Re λ < 0, excluding the algebraic power-pair rates ~ r/s):
+    #
+    #   damped zone (g_min·c·s > 8 for all growers): h = -c·s. All growing
+    #     content — including the O(smax_tol) fast-direction contamination of
+    #     the series data — sits at z ≫ 1 where L-stable Radau DAMPS it
+    #     (|R| ~ 6/z²). Slow-mode global error ≈ 10²c³ → c = (rtol/100)^{1/3};
+    #     step count ~ ln(S/s_m)/c independent of S.
+    #
+    #   resolved band (some grower at z ≲ 8): genuine backward growth must be
+    #     tracked and removed. h = -min(c·s, 0.3/g_max) keeps z ≤ 0.3 so the
+    #     DISCRETE growing eigendirection matches the continuous one to
+    #     O(z²) ~ 1% (an orthogonal purge against a frame misaligned by δ
+    #     leaks δ·content — with e-fold budget Φ ≤ 2.5 the cycle gain
+    #     e^Φ·δ < 0.2 and contamination stays at the rounding floor). The
+    #     frame F is carried CONTINUOUSLY with the same Radau steps + QR
+    #     (backward-attracting ⇒ self-correcting), so no re-extraction cost.
+    c = clamp(cbrt(rtol / 100.0), 2e-5, 2e-3)
+
+    # --- Phase 1: damped zone, W in Double64 (see radau_step_d comment). ---
+    # First purge (at S) and all damped-zone steps run in extended precision;
+    # on resolved-band entry W drops back to Float64 for phase 2.
+    Wd = CTd.(W)
+    let Floc = decaying_pair(params, Q, θ, S; growth=growth, seed=seed, dp_res=dp_res)
+        Fd = CTd.(Floc)
+        Wd .-= Fd * (Fd' * Wd)
+    end
+    s = S
+    nstep = 0
+    while s > s_m + 1e-12 * S
+        λs = eigvals(Matrix(𝔐(s)))
+        gmin = Inf
+        for λ in λs
+            g = -real(λ)
+            g > 10.0 / s && (gmin = min(gmin, g))
+        end
+        (isfinite(gmin) && gmin * c * s < 8.0) && break   # resolved-band entry
+        h = -min(c * s, s - s_m)
+        Wd = radau_step_d(s, Wd, h)
+        s += h
+        nstep += 1
+        nstep > 5_000_000 && error("march_boundary: step limit exceeded (s=$s)")
+    end
+    W = ComplexF64.(Wd)
+
+    # --- Phase 2: resolved band (Float64, carried frame + purges). ---
+    Φ = 0.0
+    F = Matrix{ComplexF64}(undef, 6, 2)
+    have_F = false
+    while s > s_m + 1e-12 * S
+        # Backward-growing exponential rates from the frozen-coefficient
+        # spectrum (Re λ < 0 ⇒ grows as s decreases); rates ≲ 10/s are the
+        # algebraic power pair, not exponential growers.
+        λs = eigvals(Matrix(𝔐(s)))
+        gmin, gmax = Inf, 0.0
+        for λ in λs
+            g = -real(λ)
+            if g > 10.0 / s
+                gmin = min(gmin, g)
+                gmax = max(gmax, g)
+            end
+        end
+        resolved = isfinite(gmin) && gmin * c * s < 8.0
+        h = resolved ? -min(c * s, ztrk / gmax, s - s_m) : -min(c * s, s - s_m)
+        if resolved && !have_F
+            F = decaying_pair(params, Q, θ, s; growth=growth, seed=seed, dp_res=dp_res)
+            W .-= F * (F' * W)
+            have_F = true
+            Φ = 0.0
+        end
+        W = radau_step(s, W, h)
+        if have_F
+            F = radau_step(s, F, h)
+            F = Matrix(qr(F).Q)
+            Φ += gmax * abs(h)
+            if Φ > purge_budget
+                W .-= F * (F' * W)
+                Φ = 0.0
+            end
+        end
+        s += h
+        nstep += 1
+        nstep > 5_000_000 && error("march_boundary: step limit exceeded (s=$s)")
+    end
+    # The boundary basis must span the traces AT s_m of the solutions that
+    # decay AT INFINITY — which, after the mode families morph around
+    # s ~ |Q|, is NOT the same subspace as the locally-decaying pair at s_m.
+    # The carried frame F is exactly the marched ∞-decaying pair (it enters
+    # the resolved band before the morphing region and tracks the subspace
+    # through it), so it IS the correct E-basis; re-extracting locally at
+    # s_m substitutes the wrong subspace and corrupts Δ at O(1) on
+    # extreme-coefficient surfaces. Local extraction remains only as the
+    # fallback when the march never entered the resolved band (no morphing
+    # crossed under resolution — the subspaces then coincide).
+    if have_F
+        W .-= F * (F' * W)
+    else
+        F = purge!(W, s_m)
+    end
+    return W[:, 1], W[:, 2], F
+end
+
+"""
+    boundary_basis(params, Q, cache, θ, S; growth=30.0)
+        -> (Us, Ub, E, cond_est)
+
+Boundary data at s = S directly from the inps series (no march; used when
+the series is already valid at the BVP edge): `Us`, `Ub` are the small/large
+power-like solutions as 6-component states in the inps normalization, `E`
+the numerically-generated decaying pair, `cond_est` the condition number of
+the column-normalized basis.
+"""
+function boundary_basis(params::GGJParameters, Q::ComplexF64,
+    cache::InnerAsymptoticsCache, θ::Float64, S::Float64;
+    growth::Float64=30.0, seed::Int=20260707, dp_res::Float64=20.0)
+    xS = cis(θ) * S
+    ua, dua = physical_ua_dua(cache, xS)
+    Ub = ComplexF64[ua[1, 1], ua[2, 1], ua[3, 1], dua[1, 1], dua[2, 1], dua[3, 1]]
+    Us = ComplexF64[ua[1, 2], ua[2, 2], ua[3, 2], dua[1, 2], dua[2, 2], dua[3, 2]]
+    E = decaying_pair(params, Q, θ, S; growth=growth, seed=seed, dp_res=dp_res)
+
+    Bmat = hcat(Us ./ norm(Us), Ub ./ norm(Ub), E)
+    cond_est = cond(Bmat)
+    return Us, Ub, E, cond_est
+end
+# solve.jl
+#
+# Global spectral-element collocation solve of the GGJ inner-layer BVP on
+# the rotated ray, with the matching data Δ as a bordered unknown.
+#
+# For each parity (deltac isol = 1 "odd", 2 "even"):
+#
+#   unknowns : v at all global nodes (6 components each), plus (Δ, c₁, c₂)
+#   equations: dv/ds = e^{iθ} M(e^{iθ}s) v   collocated at the p non-left
+#              Lobatto nodes of every cell (right-biased collocation —
+#              Radau-IIA-like, stable for the dormant stiff modes);
+#              3 parity conditions at s = 0;
+#              6 matching conditions at s = S:
+#                  v(S) − Δ·Ub − c₁E₁ − c₂E₂ = Us .
+#
+# One sparse LU per parity per mesh; residual-driven bisection refinement.
+# The exponential dichotomy never enters: no quantity is ever propagated
+# across the layer, and the boundary condition splits the dichotomy exactly.
+
+"""
+    RaySolveResult
+
+Output of [`solve_ray`](@ref). `Δ` is the rescaled matching data in the
+deltac output ordering (index 1 ↔ even-BC solution, 2 ↔ odd-BC solution,
+i.e. the same swap deltac_solve applies); `Δraw` is (isol=1, isol=2) in raw
+inps normalization before `rescale_delta`.
+"""
+struct RaySolveResult
+    Δ::SVector{2,ComplexF64}
+    Δraw::SVector{2,ComplexF64}
+    cexp::SMatrix{2,2,ComplexF64,4}   # decaying-pair coefficients (col = isol)
+    Q::ComplexF64
+    θ::Float64
+    S::Float64
+    series_ok::Bool                    # pick_smax reached its tolerance
+    bc_cond::Float64                   # boundary-basis condition number
+    breaks::Vector{Float64}
+    p::Int
+    sols::Matrix{ComplexF64}           # ndof × 2 nodal solution (isol = 1, 2)
+    resid::Float64                     # final max relative collocation residual
+    δΔ_refine::Float64                 # |Δ change| in last refinement round (rel)
+    nrounds::Int
+end
+
+# -----------------------------------------------------------------------
+# Assembly.
+# -----------------------------------------------------------------------
+
+function _assemble_base(params::GGJParameters, Q::ComplexF64, θ::Float64,
+    breaks::Vector{Float64}, p::Int,
+    t::Vector{Float64}, D::Matrix{Float64},
+    Us::Vector{ComplexF64}, Ub::Vector{ComplexF64}, E::Matrix{ComplexF64})
+    N = length(breaks) - 1
+    Ng = N * p + 1
+    ndof = 6 * Ng + 3
+    ph = cis(θ)
+    βb = norm(Ub)
+
+    nnz_est = N * p * 6 * (p + 7) + 24 + 3
+    Ic = Vector{Int}()
+    Jc = Vector{Int}()
+    Vc = Vector{ComplexF64}()
+    sizehint!(Ic, nnz_est)
+    sizehint!(Jc, nnz_est)
+    sizehint!(Vc, nnz_est)
+
+    @inbounds for c in 1:N
+        a, b = breaks[c], breaks[c+1]
+        Dscale = 2 / (b - a)
+        for j in 1:p
+            s_j = a + (b - a) * (t[j+1] + 1) / 2
+            M = ph * ode_matrix(params, Q, ph * s_j)
+            row0 = 6 * ((c - 1) * p + (j - 1))
+            gj = (c - 1) * p + j + 1
+            for i in 1:6
+                row = row0 + i
+                for k in 0:p
+                    gk = (c - 1) * p + k + 1
+                    push!(Ic, row)
+                    push!(Jc, 6 * (gk - 1) + i)
+                    push!(Vc, Dscale * D[j+1, k+1])
+                end
+                for i2 in 1:6
+                    Mv = M[i, i2]
+                    if Mv != 0
+                        push!(Ic, row)
+                        push!(Jc, 6 * (gj - 1) + i2)
+                        push!(Vc, -Mv)
+                    end
+                end
+            end
+        end
+    end
+
+    # Matching rows at s = S.
+    rowm0 = 6 * N * p
+    @inbounds for i in 1:6
+        row = rowm0 + i
+        push!(Ic, row)
+        push!(Jc, 6 * (Ng - 1) + i)
+        push!(Vc, one(ComplexF64))
+        push!(Ic, row)
+        push!(Jc, 6 * Ng + 1)
+        push!(Vc, -Ub[i] / βb)
+        push!(Ic, row)
+        push!(Jc, 6 * Ng + 2)
+        push!(Vc, -E[i, 1])
+        push!(Ic, row)
+        push!(Jc, 6 * Ng + 3)
+        push!(Vc, -E[i, 2])
+    end
+
+    rhs = zeros(ComplexF64, ndof)
+    rhs[rowm0.+(1:6)] .= Us
+
+    return Ic, Jc, Vc, rhs, βb, ndof, Ng
+end
+
+"""
+    _solve_parities(Ic, Jc, Vc, rhs, ndof, Ng, N, p) -> (sol_odd, sol_even)
+
+Solve BOTH parity systems from one factorization. The two matrices differ
+only in the 3 parity rows at s = 0 (odd: unit entries at components (4,2,3);
+even: at (1,5,6)), i.e. A_even = A_odd + Σₖ e_{rₖ}(e_{c2ₖ} − e_{c1ₖ})ᵀ — a
+rank-3 update. One sparse LU (of the row-equilibrated odd system) plus a
+Woodbury correction (3 extra triangular solves + a 3×3 solve) replaces the
+second factorization, which dominates the linear-algebra cost.
+"""
+function _solve_parities(Ic::Vector{Int}, Jc::Vector{Int}, Vc::Vector{ComplexF64},
+    rhs::Vector{ComplexF64}, ndof::Int, Ng::Int, N::Int, p::Int)
+    rowp0 = 6 * N * p + 6
+    c1 = parity_rows(1)
+    c2 = parity_rows(2)
+    I2 = copy(Ic)
+    J2 = copy(Jc)
+    V2 = copy(Vc)
+    for k in 1:3
+        push!(I2, rowp0 + k)
+        push!(J2, c1[k])
+        push!(V2, one(ComplexF64))
+    end
+    A = sparse(I2, J2, V2, ndof, ndof)
+
+    # Row equilibration (as in _solve_parity; the parity rows have unit
+    # entries in both variants, so one scaling serves both systems).
+    rmax = zeros(Float64, ndof)
+    rows = rowvals(A)
+    vals = nonzeros(A)
+    for col in 1:ndof
+        for idx in nzrange(A, col)
+            r = rows[idx]
+            a = abs(vals[idx])
+            a > rmax[r] && (rmax[r] = a)
+        end
+    end
+    @inbounds for i in 1:ndof
+        rmax[i] = rmax[i] > 0 ? 1 / rmax[i] : 1.0
+    end
+    Dr = Diagonal(rmax)
+    F = lu(Dr * A)
+    x1 = F \ (Dr * rhs)
+
+    # Woodbury: (DrA + U Vᵀ)⁻¹ b = x1 − A⁻¹U (I₃ + Vᵀ A⁻¹U)⁻¹ Vᵀ x1,
+    # with U[:,k] = Dr[rₖ,rₖ]·e_{rₖ} and Vᵀ[k,:] = (e_{c2ₖ} − e_{c1ₖ})ᵀ.
+    U = zeros(ComplexF64, ndof, 3)
+    for k in 1:3
+        U[rowp0+k, k] = rmax[rowp0+k]
+    end
+    AinvU = F \ U
+    Vx1 = ComplexF64[x1[c2[k]] - x1[c1[k]] for k in 1:3]
+    S = Matrix{ComplexF64}(I, 3, 3)
+    for k in 1:3, j in 1:3
+        S[k, j] += AinvU[c2[k], j] - AinvU[c1[k], j]
+    end
+    x2 = x1 - AinvU * (S \ Vx1)
+    # UMFPACK numeric factorizations live in C malloc, invisible to the GC's
+    # heap accounting — under the refinement loop, lazily-finalized multi-GB
+    # factorizations pile up (observed: tens of GB on 24k-cell meshes). Free
+    # eagerly now that both solutions are extracted.
+    finalize(F)
+    return x1, x2
+end
+
+function _solve_parity(Ic::Vector{Int}, Jc::Vector{Int}, Vc::Vector{ComplexF64},
+    rhs::Vector{ComplexF64}, ndof::Int, Ng::Int, N::Int, p::Int, isol::Int;
+    leftcomps::AbstractVector{Int}=parity_rows(isol),
+    leftvals::AbstractVector{ComplexF64}=zeros(ComplexF64, 3))
+    rowp0 = 6 * N * p + 6
+    I2 = copy(Ic)
+    J2 = copy(Jc)
+    V2 = copy(Vc)
+    rhs = copy(rhs)
+    for k in eachindex(leftcomps)
+        push!(I2, rowp0 + k)
+        push!(J2, leftcomps[k])      # node 1 → columns 1..6
+        push!(V2, one(ComplexF64))
+        rhs[rowp0+k] = leftvals[k]
+    end
+    A = sparse(I2, J2, V2, ndof, ndof)
+
+    # Row equilibration.
+    rmax = zeros(Float64, ndof)
+    rows = rowvals(A)  # CSC: iterate columns
+    vals = nonzeros(A)
+    for col in 1:ndof
+        for idx in nzrange(A, col)
+            r = rows[idx]
+            a = abs(vals[idx])
+            a > rmax[r] && (rmax[r] = a)
+        end
+    end
+    @inbounds for i in 1:ndof
+        rmax[i] = rmax[i] > 0 ? 1 / rmax[i] : 1.0
+    end
+    Dr = Diagonal(rmax)
+    F = lu(Dr * A)
+    x = F \ (Dr * rhs)
+    finalize(F)   # eager UMFPACK free — see _solve_parities
+    return x
+end
+
+# -----------------------------------------------------------------------
+# Residual estimation for adaptive refinement.
+# -----------------------------------------------------------------------
+
+function _cell_errors(params::GGJParameters, Q::ComplexF64, θ::Float64,
+    breaks::Vector{Float64}, p::Int,
+    t::Vector{Float64}, D::Matrix{Float64}, w::Vector{Float64},
+    sol::Vector{ComplexF64})
+    N = length(breaks) - 1
+    ph = cis(θ)
+    errs = zeros(Float64, N)
+
+    # Reference-cell test points: left endpoint + inter-node midpoints.
+    τs = Float64[-1.0]
+    for j in 1:p
+        push!(τs, (t[j] + t[j+1]) / 2)
+    end
+
+    vals = Matrix{ComplexF64}(undef, p + 1, 6)
+    for c in 1:N
+        a, b = breaks[c], breaks[c+1]
+        Dscale = 2 / (b - a)
+        for j in 0:p
+            g = (c - 1) * p + j + 1
+            for i in 1:6
+                vals[j+1, i] = sol[6*(g-1)+i]
+            end
+        end
+        dvals = (Dscale .* D) * vals
+
+        emax = 0.0
+        for τ in τs
+            vτ = barycentric_eval(t, w, vals, τ)
+            dvτ = barycentric_eval(t, w, dvals, τ)
+            s_τ = a + (b - a) * (τ + 1) / 2
+            Mv = (ph * ode_matrix(params, Q, ph * s_τ)) * vτ
+            rnorm = 0.0
+            scale = 1e-300
+            for i in 1:6
+                rnorm = max(rnorm, abs(dvτ[i] - Mv[i]))
+                scale = max(scale, abs(dvτ[i]), abs(Mv[i]))
+            end
+            emax = max(emax, rnorm / scale)
+        end
+        errs[c] = emax
+    end
+    return errs
+end
+
+# -----------------------------------------------------------------------
+# Driver.
+# -----------------------------------------------------------------------
+
+"""
+    solve_ray(params::GGJParameters, Q; θ=angle(Q)/4, kmax=12, p=12,
+              S=nothing, smax_tol=1e-9, growth=30.0,
+              refine_tol=1e-8, max_rounds=8, max_cells=6000,
+              verbose=false) -> RaySolveResult
+
+Solve the GGJ inner-layer matching problem on the rotated ray x = e^{iθ}s
+and return the parity matching data. `Δ` follows the deltac output
+convention (swap + `rescale_delta`) so it is directly comparable to the
+GPEC_julia Galerkin solver.
+
+Keyword highlights:
+  - `θ`     : contour angle; default arg(Q)/4 makes the parabolic exponent
+              exactly real. θ = 0 reproduces a real-axis solve.
+  - `S`     : matching radius; default from the inps series residual
+              criterion along the ray (`pick_smax`).
+  - `p`     : polynomial order per spectral element.
+  - `refine_tol`/`max_rounds`: residual-driven bisection refinement control.
+"""
+function solve_ray(params::GGJParameters, Q::ComplexF64;
+    θ::Float64=angle(Q) / 4, kmax::Int=12, p::Int=12,
+    S::Union{Nothing,Float64}=nothing, smax_tol::Float64=1e-9,
+    growth::Float64=30.0, refine_tol::Float64=1e-8,
+    march_rtol::Float64=1e-9, sm_fac::Float64=1.0, cfl::Float64=0.4,
+    max_rounds::Int=8, max_cells::Int=6000, seed::Int=20260707,
+    ztrk::Float64=0.3, purge_budget::Float64=2.5, dp_res::Float64=20.0,
+    handoff::Union{Nothing,NTuple{2,Vector{ComplexF64}}}=nothing,
+    boundary::Union{Nothing,Tuple{Vector{ComplexF64},Vector{ComplexF64},Matrix{ComplexF64}}}=nothing,
+    verbose::Bool=false)
+
+    cache = build_asymptotics(params, Q; kmax=kmax)
+    series_ok = true
+    if S === nothing
+        Sv, cache, series_ok = pick_smax(params, Q; θ=θ, eps=smax_tol, kmax=kmax, cache=cache)
+        series_ok || @warn "pick_smax: series residual target $smax_tol not reached; using best point S=$Sv"
+    else
+        Sv = S
+    end
+
+    # Decide the BVP outer radius: the layer physics (exponential content
+    # above machine floor) ends at s_d; beyond that the solution is pure
+    # power law and is handled by the projected quotient MARCH, not by
+    # collocation. For small |Q| (s_d ≳ S) no march is needed.
+    sq = sqrt(Q)
+    ρ = real(cis(2θ) / sq)
+    # BVP outer radius s_m: where every exponential solution family has
+    # decayed ≥ 45 e-folds from the origin, computed from the actual
+    # frozen-coefficient spectrum (the asymptotic ρs rate is wrong near the
+    # origin for extreme-coefficient surfaces, where the Υ-family rate
+    # √|Q(G+KF)| dominates and kills exponential content within s ~ O(1)).
+    s_m = let acc = 0.0, sprev = 0.0, out = Sv
+        for sk in exp10.(range(-2, log10(Sv), length=400))
+            gdec = Inf
+            for λ in eigvals(Matrix(cis(θ) * ode_matrix(params, Q, cis(θ) * sk)))
+                g = -real(λ)                    # backward-growth = forward-decay rate
+                g > 10.0 / sk && (gdec = min(gdec, g))
+            end
+            isfinite(gdec) && (acc += gdec * (sk - sprev))
+            sprev = sk
+            if acc >= 45.0
+                out = sk
+                break
+            end
+        end
+        sm_fac * max(out, 2.0)
+    end
+
+    local Us, Ub, E, bc_cond, Sbvp
+    if boundary !== nothing
+        # diagnostic: externally supplied boundary data (Us, Ub, E) at s_m,
+        # e.g. from an extended-precision march
+        Us, Ub, E = boundary
+        Sbvp = s_m
+        bc_cond = cond(hcat(Us ./ norm(Us), Ub ./ norm(Ub), E))
+    elseif s_m < 0.98 * Sv
+        local UsS, UbS
+        if handoff === nothing
+            ua, dua = physical_ua_dua(cache, cis(θ) * Sv)
+            UbS = ComplexF64[ua[1, 1], ua[2, 1], ua[3, 1], dua[1, 1], dua[2, 1], dua[3, 1]]
+            UsS = ComplexF64[ua[1, 2], ua[2, 2], ua[3, 2], dua[1, 2], dua[2, 2], dua[3, 2]]
+        else
+            # diagnostic: externally supplied series hand-off (Us, Ub) at
+            # x = e^{iθ}S, e.g. evaluated in extended precision
+            UsS, UbS = handoff
+        end
+        Us, Ub, E = march_boundary(params, Q, θ, Sv, s_m, UsS, UbS;
+            rtol=march_rtol, growth=growth, seed=seed,
+            ztrk=ztrk, purge_budget=purge_budget, dp_res=dp_res)
+        Sbvp = s_m
+        bc_cond = cond(hcat(Us ./ norm(Us), Ub ./ norm(Ub), E))
+        verbose && @printf("  march: S=%.4g → s_m=%.4g, bc_cond=%.1e\n", Sv, s_m, bc_cond)
+    else
+        Us, Ub, E, bc_cond = boundary_basis(params, Q, cache, θ, Sv;
+            growth=growth, seed=seed, dp_res=dp_res)
+        Sbvp = Sv
+    end
+    bc_cond < 1e12 || @warn "boundary basis poorly conditioned (cond ≈ $bc_cond)"
+
+    t, D = cheblobatto(p)
+    w = barycentric_weights(t)
+    breaks = initial_breaks(params, Q, θ, Sbvp, p; cfl=cfl)
+
+    Δraw = MVector{2,ComplexF64}(0, 0)
+    Δprev = MVector{2,ComplexF64}(NaN, NaN)
+    cexp = zeros(MMatrix{2,2,ComplexF64})
+    local sols
+    resid = Inf
+    resid_prev = Inf
+    nstall = 0
+    δΔ = Inf
+    round_used = 0
+
+    for round in 0:max_rounds
+        round_used = round
+        N = length(breaks) - 1
+        Ic, Jc, Vc, rhs, βb, ndof, Ng = _assemble_base(params, Q, θ, breaks, p, t, D, Us, Ub, E)
+        sols = Matrix{ComplexF64}(undef, ndof, 2)
+        errs = zeros(Float64, N)
+        sol_odd, sol_even = _solve_parities(Ic, Jc, Vc, rhs, ndof, Ng, N, p)
+        for (isol, sol) in ((1, sol_odd), (2, sol_even))
+            sols[:, isol] = sol
+            Δraw[isol] = sol[6*Ng+1] / βb
+            cexp[1, isol] = sol[6*Ng+2]
+            cexp[2, isol] = sol[6*Ng+3]
+            errs = max.(errs, _cell_errors(params, Q, θ, breaks, p, t, D, w, sol))
+        end
+        resid = maximum(errs)
+        δΔ = maximum(abs.(Δraw .- Δprev) ./ (abs.(Δraw) .+ 1e-300))
+        Δprev .= Δraw
+        verbose && @printf("  round %d: %d cells, resid %.2e, Δraw[1] = %.6e %+.6eim, δΔ = %.1e\n",
+            round, N, resid, real(Δraw[1]), imag(Δraw[1]), δΔ)
+
+        (resid < refine_tol || round == max_rounds) && break
+        # Stagnation guard: on extreme-coefficient surfaces the residual
+        # bottoms out at the roundoff floor ε·‖M‖ (~10⁻⁸ ≫ refine_tol);
+        # refining past that plateau burns cells without informing Δ. Two
+        # consecutive rounds with < 2× improvement ⇒ converged-to-plateau.
+        nstall = resid > 0.5 * resid_prev ? nstall + 1 : 0
+        resid_prev = resid
+        if nstall >= 2
+            verbose && @printf("  refinement stagnated at resid %.2e (roundoff plateau) — stopping\n", resid)
+            break
+        end
+        marked = findall(>(refine_tol), errs)
+        isempty(marked) && break
+        if N + length(marked) > max_cells
+            @warn "solve_ray: cell budget ($max_cells) reached with resid=$resid"
+            break
+        end
+        newbreaks = Float64[breaks[1]]
+        for c in 1:N
+            if c in marked
+                push!(newbreaks, (breaks[c] + breaks[c+1]) / 2)
+            end
+            push!(newbreaks, breaks[c+1])
+        end
+        breaks = newbreaks
+    end
+
+    Δswapped = SVector{2,ComplexF64}(Δraw[2], Δraw[1])   # deltac_solve convention
+    Δ = rescale_delta(Δswapped, params)
+
+    return RaySolveResult(Δ, SVector(Δraw), SMatrix(cexp), Q, θ, Sv, series_ok,
+        bc_cond, breaks, p, sols, resid, δΔ, round_used)
+end
+
+solve_ray(params::GGJParameters, Q::Number; kwargs...) =
+    solve_ray(params, ComplexF64(Q); kwargs...)
+
+"""
+    solve_inner(::GGJModel{:ray}, params::GGJParameters, γ::Number; kwargs...)
+        -> SVector{2,ComplexF64}
+
+Solve the GGJ inner-layer matching problem with the rotated-ray collocation
+backend and return the parity-projected matching data `(Δ₁, Δ₂)` in the same
+convention as the `:shooting` and `:galerkin` backends (deltac output swap +
+`X₀^{2√(−D_I)}` physical rescale applied), directly interchangeable with
+them. Converts γ via `inner_Q` and forwards all keywords to
+[`solve_ray`](@ref); use `solve_ray` directly when the full
+[`RaySolveResult`](@ref) (raw Δ, contour, mesh, nodal solution, residuals)
+is wanted.
+
+Preferred backend for |Q| ≳ 1 and for Q near the imaginary axis (rotation /
+resistivity scans): validated to |Q| = 500 on the imaginary axis, where the
+other backends fail. At |Q| ≪ 1 all three backends agree; `:ray` remains
+accurate but `:galerkin` may be faster for very small |Q| real Q.
+"""
+function solve_inner(::GGJModel{:ray}, params::GGJParameters, γ::Number; kwargs...)
+    res = solve_ray(params, inner_Q(params, γ); kwargs...)
+    return res.Δ
+end
+
+"""
+    evaluate_solution(res::RaySolveResult, s::Real; isol=1) -> SVector{6,ComplexF64}
+
+Evaluate the solved state v = (Ψ, Ξ, Υ, Ψ', Ξ', Υ') at ray parameter `s`
+(diagnostic; barycentric interpolation on the containing cell).
+"""
+function evaluate_solution(res::RaySolveResult, s::Real; isol::Int=1)
+    breaks = res.breaks
+    p = res.p
+    (breaks[1] <= s <= breaks[end]) || throw(ArgumentError("s=$s outside [0, $(breaks[end])]"))
+    c = clamp(searchsortedlast(breaks, s), 1, length(breaks) - 1)
+    a, b = breaks[c], breaks[c+1]
+    t, _ = cheblobatto(p)
+    w = barycentric_weights(t)
+    vals = Matrix{ComplexF64}(undef, p + 1, 6)
+    for j in 0:p
+        g = (c - 1) * p + j + 1
+        for i in 1:6
+            vals[j+1, i] = res.sols[6*(g-1)+i, isol]
+        end
+    end
+    τ = 2 * (s - a) / (b - a) - 1
+    return SVector{6,ComplexF64}(barycentric_eval(t, w, vals, τ)...)
+end
+
+"""
+    solution_profile(res::RaySolveResult; npc=8) -> (; s, x, Ψ, Ξ, Υ)
+
+Reconstruct the physical inner-layer fields (Ψ, Ξ, Υ) on the solution
+contour, `npc` points per cell over the collocation domain [0, s_m].
+Output layout mirrors the legacy `solve_inner_profile`: each field is an
+`npts × 2` matrix with columns (isol=1 "odd", isol=2 "even"); `x = e^{iθ}s`
+is the (complex) layer coordinate. For real Q (θ = 0) this is the real-axis
+profile directly; for rotated solves the fields live on the ray — analytic
+continuation back to real x is a separate (unstable) problem, by design.
+"""
+function solution_profile(res::RaySolveResult; npc::Int=8)
+    breaks = res.breaks
+    p = res.p
+    t, _ = cheblobatto(p)
+    w = barycentric_weights(t)
+    N = length(breaks) - 1
+    npts = N * npc + 1
+    s = Vector{Float64}(undef, npts)
+    Ψ = Matrix{ComplexF64}(undef, npts, 2)
+    Ξ = Matrix{ComplexF64}(undef, npts, 2)
+    Υ = Matrix{ComplexF64}(undef, npts, 2)
+    vals = Matrix{ComplexF64}(undef, p + 1, 6)
+    k = 0
+    for c in 1:N
+        a, b = breaks[c], breaks[c+1]
+        for m in 0:(npc-1)
+            k += 1
+            τ = -1 + 2 * m / npc
+            s[k] = a + (b - a) * (τ + 1) / 2
+            for isol in 1:2
+                for j in 0:p, i in 1:6
+                    vals[j+1, i] = res.sols[6*((c-1)*p+j)+i, isol]
+                end
+                v = barycentric_eval(t, w, vals, τ)
+                Ψ[k, isol] = v[1]
+                Ξ[k, isol] = v[2]
+                Υ[k, isol] = v[3]
+            end
+        end
+    end
+    # right endpoint
+    k += 1
+    s[k] = breaks[end]
+    for isol in 1:2
+        v = evaluate_solution(res, breaks[end]; isol=isol)
+        Ψ[k, isol] = v[1]
+        Ξ[k, isol] = v[2]
+        Υ[k, isol] = v[3]
+    end
+    return (; s=s, x=cis(res.θ) .* s, Ψ=Ψ, Ξ=Ξ, Υ=Υ)
+end
+
+"""
+    asymptotic_profile(params, res::RaySolveResult, srange) -> (; s, x, Ψ, Ξ, Υ)
+
+Evaluate the ANALYTIC representation `u_small + Δ·u_big` on the ray for
+`s ≥ res.S` (where the inps series is trusted; decaying-exponential content
+is machine-dead there). Concatenating with [`solution_profile`](@ref)
+demonstrates the seamless numeric ↔ asymptotic match — the same overlap the
+outer-region matching relies on.
+"""
+function asymptotic_profile(params::GGJParameters, res::RaySolveResult,
+    srange::AbstractVector{<:Real})
+    cache = build_asymptotics(params, res.Q; kmax=12)
+    n = length(srange)
+    Ψ = Matrix{ComplexF64}(undef, n, 2)
+    Ξ = Matrix{ComplexF64}(undef, n, 2)
+    Υ = Matrix{ComplexF64}(undef, n, 2)
+    for (k, s) in enumerate(srange)
+        ua, _ = physical_ua_dua(cache, cis(res.θ) * s)
+        for isol in 1:2
+            Δ = res.Δraw[isol]
+            Ψ[k, isol] = ua[1, 2] + Δ * ua[1, 1]
+            Ξ[k, isol] = ua[2, 2] + Δ * ua[2, 1]
+            Υ[k, isol] = ua[3, 2] + Δ * ua[3, 1]
+        end
+    end
+    return (; s=collect(float.(srange)), x=cis(res.θ) .* srange, Ψ=Ψ, Ξ=Ξ, Υ=Υ)
+end
+
+"""
+    delta_convergence(params, Q; verbose=true, refine_tol=1e-9,
+                      max_rounds=16, max_cells=12000, kwargs...)
+        -> (; Δ, spread, base, table)
+
+Convergence battery for the matching data: solve at a trusted baseline, then
+re-solve with every numerical knob perturbed on an INDEPENDENT axis, and
+report the relative change of (Δ₁, Δ₂) per knob. The knobs probe orthogonal
+error sources, so the worst-case `spread` is an honest error bar on Δ:
+
+  - `θ → 1.2θ`       : contour angle — Δ is an analytic invariant of the ray,
+                       so any drift is pure numerical error (sharpest test).
+                       OUTWARD: the natural angle maximizes clearance from the
+                       x² ≈ −Q²(G+KF) pseudo-resonance, so inward checks
+                       measure their own degraded solve, not the baseline;
+  - `p → p+4`        : spectral-element order (interior discretization);
+  - `kmax+4, smax_tol/10` : series order and matching radius S (far-field BC);
+  - `refine_tol/10`  : mesh refinement depth;
+  - `march_rtol/100` : outer-region quotient-march tolerance;
+  - `sm_fac 1→1.4`   : BVP/march handoff radius;
+  - `growth 30→45`   : decaying-pair purification e-folds.
+
+Returns the baseline `Δ`, the per-parity worst-case relative `spread`, the
+baseline result struct, and the per-knob table `(name, δΔ₁, δΔ₂)`.
+"""
+function delta_convergence(params::GGJParameters, Q::ComplexF64;
+    verbose::Bool=true, refine_tol::Float64=1e-9,
+    max_rounds::Int=16, max_cells::Int=12000, kwargs...)
+    base = solve_ray(params, Q; refine_tol=refine_tol, max_rounds=max_rounds,
+        max_cells=max_cells, kwargs...)
+    variations = [
+        ("θ → 1.2θ", (; θ=(base.θ == 0 ? 0.1 : 1.2 * base.θ))),
+        ("p → p+4", (; p=base.p + 4)),
+        ("kmax+4, S↑", (; kmax=16, smax_tol=1e-10)),
+        ("refine_tol/10", (; refine_tol=refine_tol / 10)),
+        ("march_rtol/100", (; march_rtol=1e-11)),
+        ("sm_fac 1→1.4", (; sm_fac=1.4)),
+        ("growth 30→45", (; growth=45.0)),
+    ]
+    table = NamedTuple[]
+    spread = MVector{2,Float64}(0.0, 0.0)
+    verbose && @printf("Δ convergence battery, Q = %s:\n", string(Q))
+    verbose && @printf("  baseline: Δ₁ = %+.8e %+.8eim, Δ₂ = %+.8e %+.8eim\n",
+        real(base.Δ[1]), imag(base.Δ[1]), real(base.Δ[2]), imag(base.Δ[2]))
+    for (name, kw) in variations
+        r = solve_ray(params, Q; refine_tol=refine_tol, max_rounds=max_rounds,
+            max_cells=max_cells, kwargs..., kw...)
+        d1 = abs(r.Δ[1] - base.Δ[1]) / abs(base.Δ[1])
+        d2 = abs(r.Δ[2] - base.Δ[2]) / abs(base.Δ[2])
+        spread[1] = max(spread[1], d1)
+        spread[2] = max(spread[2], d2)
+        push!(table, (; name, d1, d2))
+        verbose && @printf("  %-16s δΔ₁ = %.2e   δΔ₂ = %.2e\n", name, d1, d2)
+    end
+    verbose && @printf("  %-16s δΔ₁ = %.2e   δΔ₂ = %.2e\n", "WORST CASE", spread[1], spread[2])
+    return (; Δ=base.Δ, spread=SVector(spread), base=base, table=table)
+end
+
+delta_convergence(params::GGJParameters, Q::Number; kwargs...) =
+    delta_convergence(params, ComplexF64(Q); kwargs...)
+
+export solve_inner_ray, evaluate_solution, delta_convergence,
+    solution_profile, asymptotic_profile
diff --git a/src/InnerLayer/GGJ/RayAsymptotics.jl b/src/InnerLayer/GGJ/RayAsymptotics.jl
new file mode 100644
index 00000000..06deeaf1
--- /dev/null
+++ b/src/InnerLayer/GGJ/RayAsymptotics.jl
@@ -0,0 +1,244 @@
+# RayAsymptotics.jl
+#
+# Complex-x extensions of the inps Wasow asymptotic kernel
+# (InnerAsymptotics.jl) for the rotated-ray solver (Ray.jl).
+#
+# The inps construction (GW2020 Eqs. 3–52) depends only on Q — the contour
+# enters only through the EVALUATION point. These methods therefore reuse the
+# existing, validated `InnerAsymptoticsCache` and coefficient packing verbatim
+# and add evaluation at complex x = e^{iθ}s (principal branch of x^r, valid
+# for |arg x| < π), the plain-variable conversion used by the ray boundary
+# data, the GW2020 Eq. (54) residual along the ray, and the series-radius
+# selector `pick_smax` (the ray analog of `pick_xmax`).
+#
+# Ported from GGJRay (~/Projects/GGJ_test, docs/METHOD.md §6); the real-x
+# methods of InnerAsymptotics.jl are untouched.
+
+# Complex-argument Horner evaluator: same contract as _horner(x::Real, ...)
+# in InnerAsymptotics.jl, with x^rvec on the principal branch.
+function _horner(x::Complex, c::AbstractMatrix{ComplexF64};
+    rvec::Union{Nothing,AbstractVector{<:Real}}=nothing,
+    derivative::Bool=false)
+    nrows, ncols = size(c)
+    n = ncols - 1
+
+    y = Vector{ComplexF64}(undef, nrows)
+    @inbounds for i in 1:nrows
+        y[i] = c[i, ncols]
+    end
+    @inbounds for k in (n-1):-1:0
+        for i in 1:nrows
+            y[i] = y[i] * x + c[i, k+1]
+        end
+    end
+
+    if rvec !== nothing
+        @inbounds for i in 1:nrows
+            y[i] *= ComplexF64(x)^rvec[i]
+        end
+    end
+
+    if !derivative
+        return y, nothing
+    end
+
+    dy = similar(y)
+    if rvec !== nothing
+        @inbounds for i in 1:nrows
+            dy[i] = c[i, ncols] * (rvec[i] + n)
+        end
+        @inbounds for k in (n-1):-1:0
+            for i in 1:nrows
+                dy[i] = dy[i] * x + c[i, k+1] * (rvec[i] + k)
+            end
+        end
+        @inbounds for i in 1:nrows
+            dy[i] = dy[i] * ComplexF64(x)^rvec[i] / x
+        end
+    else
+        @inbounds for i in 1:nrows
+            dy[i] = c[i, ncols] * n
+        end
+        @inbounds for k in (n-1):-1:0
+            for i in 1:nrows
+                dy[i] = dy[i] * x + c[i, k+1] * k
+            end
+        end
+        @inbounds for i in 1:nrows
+            dy[i] = dy[i] / x
+        end
+    end
+    return y, dy
+end
+
+"""
+    evaluate_asymptotics(cache, x::Complex; derivative=true, apply_T=true) -> (U, dU)
+
+Evaluate the inps asymptotic basis at complex `x` (GW2020 Eq. 53; principal
+branch, |arg x| < π). Returns the 6×2 matrix `U` whose columns are the two
+power-like asymptotic solutions in the GW2020 Eq. (2) state convention
+`u = (xΨ, Ξ, Υ, (xΨ)'/x, Ξ'/x, Υ'/x)`, and their x-derivatives `dU` if
+requested. `apply_T=false` returns the pre-Jordan-basis representation used
+by the residual measure.
+"""
+function evaluate_asymptotics(cache::InnerAsymptoticsCache, x::Complex;
+    derivative::Bool=true, apply_T::Bool=true)
+    xfac = 1 / (x * x)
+    R = cache.R
+    rvec = SVector(-R[1] / 2, -R[1] / 2, -R[2] / 2, -R[2] / 2)
+
+    cc = _pack_y_coefs(cache)
+    dd = _pack_qp_coefs(cache)
+
+    yy, dyy = _horner(xfac, cc; rvec=rvec, derivative=derivative)
+    zz, dzz = _horner(xfac, dd; derivative=derivative)
+
+    y = SMatrix{2,2,ComplexF64}(yy[1], yy[2], yy[3], yy[4])
+    q = SMatrix{2,2,ComplexF64}(zz[1], zz[2], zz[3], zz[4])
+    p21 = SMatrix{4,2,ComplexF64}(zz[5], zz[6], zz[7], zz[8],
+        zz[9], zz[10], zz[11], zz[12])
+
+    pp_m = zeros(ComplexF64, 6, 2)
+    pp_m[1, 1] = 1
+    pp_m[2, 2] = 1
+    @inbounds for i in 1:4, j in 1:2
+        pp_m[i+2, j] = p21[i, j]
+    end
+    pp = SMatrix{6,2,ComplexF64}(pp_m)
+
+    smat = @SMatrix ComplexF64[1 0; 0 xfac]
+
+    qsy = q * smat * y
+    U = pp * qsy
+
+    if derivative
+        chain = -2 * xfac / x
+        dy_v = chain .* dyy
+        dz_v = chain .* dzz
+
+        dy = SMatrix{2,2,ComplexF64}(dy_v[1], dy_v[2], dy_v[3], dy_v[4])
+        dq = SMatrix{2,2,ComplexF64}(dz_v[1], dz_v[2], dz_v[3], dz_v[4])
+        dp21 = SMatrix{4,2,ComplexF64}(dz_v[5], dz_v[6], dz_v[7], dz_v[8],
+            dz_v[9], dz_v[10], dz_v[11], dz_v[12])
+
+        dpp_m = zeros(ComplexF64, 6, 2)
+        @inbounds for i in 1:4, j in 1:2
+            dpp_m[i+2, j] = dp21[i, j]
+        end
+        dpp = SMatrix{6,2,ComplexF64}(dpp_m)
+
+        dsmat = @SMatrix ComplexF64[0 0; 0 (-2*xfac/x)]
+        dqsy = q * smat * dy + q * dsmat * y + dq * smat * y
+        dU = pp * dqsy + dpp * qsy
+
+        if apply_T
+            U = cache.T * U
+            dU = cache.T * dU
+        end
+        return U, dU
+    else
+        if apply_T
+            U = cache.T * U
+        end
+        return U, nothing
+    end
+end
+
+"""
+    physical_ua_dua(cache, x::Number) -> (ua, dua)
+
+Convert the inps 6×2 basis at (possibly complex) `x` to physical `(Ψ, Ξ, Υ)`
+values and x-derivatives (each 3×2; column 1 = large power solution,
+column 2 = small). Same map as the deltac `inpso_get_ua/dua` convention and
+the Galerkin backend's `_physical_ua_dua`, generalized to complex x.
+"""
+function physical_ua_dua(cache::InnerAsymptoticsCache, x::Number)
+    xc = ComplexF64(x)
+    U, dU = evaluate_asymptotics(cache, xc; derivative=true, apply_T=true)
+    ua = zeros(ComplexF64, 3, 2)
+    dua = zeros(ComplexF64, 3, 2)
+    @inbounds for j in 1:2
+        ua[1, j] = U[1, j] / xc
+        ua[2, j] = U[2, j]
+        ua[3, j] = U[3, j]
+        dua[1, j] = U[4, j] - U[1, j] / (xc * xc)
+        dua[2, j] = U[5, j] * xc
+        dua[3, j] = U[6, j] * xc
+    end
+    return ua, dua
+end
+
+"""
+    asymptotic_residual(cache, x::Complex) -> SVector{2,Float64}
+
+GW2020 Eq. (54) convergence measure of the two power-like series columns at
+complex `x`: ‖dU − x·J(x)·U‖∞ / max(‖dU‖∞, ‖x·J(x)·U‖∞), per column.
+"""
+function asymptotic_residual(cache::InnerAsymptoticsCache, x::Complex)
+    U, dU = evaluate_asymptotics(cache, x; derivative=true, apply_T=false)
+
+    xfac = 1 / (x * x)
+    M = cache.J[1]
+    if cache.kmax > 0
+        M = M + xfac * cache.J[2]
+    end
+    if cache.kmax > 1
+        M = M + xfac * xfac * cache.J[3]
+    end
+
+    matvec1 = dU
+    matvec2 = -x * (M * U)
+    matvec0 = matvec1 + matvec2
+
+    delta = MVector{2,Float64}(0.0, 0.0)
+    @inbounds for j in 1:2
+        n0 = 0.0
+        n1 = 0.0
+        n2 = 0.0
+        for i in 1:6
+            n0 = max(n0, abs(matvec0[i, j]))
+            n1 = max(n1, abs(matvec1[i, j]))
+            n2 = max(n2, abs(matvec2[i, j]))
+        end
+        denom = max(n1, n2)
+        delta[j] = denom == 0 ? 0.0 : n0 / denom
+    end
+    return SVector{2,Float64}(delta)
+end
+
+"""
+    pick_smax(params, Q; θ=angle(Q)/4, eps=1e-9, kmax=12, cache=nothing,
+              slogmin=-1.0, slogmax=6.5, dslog=0.01) -> (S, cache, achieved)
+
+Ray analog of `pick_xmax`: sweep the ray parameter `s` log-uniformly and
+return the smallest `s` at which the series residual along `x = e^{iθ}s`
+drops below `eps`. If the target is never reached, returns the residual-
+minimizing `s` with `achieved = false` (caller should warn). Also returns
+the asymptotics cache for reuse.
+"""
+function pick_smax(params::GGJParameters, Q::ComplexF64;
+    θ::Float64=angle(Q) / 4, eps::Float64=1e-9, kmax::Int=12,
+    cache::Union{Nothing,InnerAsymptoticsCache}=nothing,
+    slogmin::Float64=-1.0, slogmax::Float64=6.5, dslog::Float64=0.01)
+    c = cache === nothing ? build_asymptotics(params, Q; kmax=kmax) : cache
+    ph = cis(θ)
+    best_s = 10.0^slogmin
+    best_r = Inf
+    slog = slogmin
+    while slog <= slogmax
+        s = 10.0^slog
+        r = maximum(asymptotic_residual(c, ph * s))
+        if r < eps
+            return s, c, true
+        end
+        if r < best_r
+            best_r = r
+            best_s = s
+        end
+        slog += dslog
+    end
+    return best_s, c, false
+end
+
+pick_smax(params::GGJParameters, Q::Number; kwargs...) =
+    pick_smax(params, ComplexF64(Q); kwargs...)
diff --git a/src/InnerLayer/InnerLayer.jl b/src/InnerLayer/InnerLayer.jl
index 934e8f83..03b195ce 100644
--- a/src/InnerLayer/InnerLayer.jl
+++ b/src/InnerLayer/InnerLayer.jl
@@ -17,6 +17,8 @@ include("GGJ/GGJ.jl")
 import .GGJ: GGJModel, GGJParameters, build_asymptotics, evaluate_asymptotics, pick_xmax
 import .GGJ: InnerAsymptoticsCache, mercier_di, mercier_dr, inner_Q, rescale_delta
 import .GGJ: glasser_wang_2020_eq55, solve_inner_converged  # solve_inner_converged: experimental, not exported (reachable as a qualified call only)
+import .GGJ: solve_ray, RaySolveResult, pick_smax, physical_ua_dua
+import .GGJ: delta_convergence, solution_profile, asymptotic_profile, q4_surface_benchmark
 # SLAYER imports go here
 
 export InnerLayerModel, solve_inner
@@ -24,6 +26,8 @@ export GGJ, GGJModel, GGJParameters
 export build_asymptotics, evaluate_asymptotics, pick_xmax, InnerAsymptoticsCache
 export mercier_di, mercier_dr, inner_Q, rescale_delta
 export glasser_wang_2020_eq55
+export solve_ray, RaySolveResult, pick_smax, physical_ua_dua
+export delta_convergence, solution_profile, asymptotic_profile, q4_surface_benchmark
 
 # SLAYER exports go here
 
diff --git a/test/runtests_innerlayer.jl b/test/runtests_innerlayer.jl
index 7ee36d06..f5e3f7cc 100644
--- a/test/runtests_innerlayer.jl
+++ b/test/runtests_innerlayer.jl
@@ -49,3 +49,95 @@ const GGJ = IL.GGJ
         @test abs(imag(Δ[2])) < 1e-3 * abs(Δ[2])
     end
 end
+
+@testset "InnerLayer GGJ :ray backend (rotated-contour collocation)" begin
+    # Ported 2026-07-08 from the standalone GGJRay package (validated there:
+    # manufactured Δ* to 3e-14, Fortran rmatch pins, 96-equilibrium robustness
+    # scans to Q = 500i). These tests pin the port, not the method.
+    p = IL.glasser_wang_2020_eq55()
+
+    @testset "agrees with :galerkin at the paper point Q = 0.1234" begin
+        γ = 0.1234 * GGJ.q0(p)
+        Δ = IL.solve_inner(IL.GGJModel(; solver=:ray), p, γ)
+        # Same Fortran-cross-checked pins as the Galerkin testset above.
+        @test real(Δ[1]) ≈ 3.698368e4 rtol = 1e-3
+        @test real(Δ[2]) ≈ 14.759721 rtol = 1e-3
+        @test abs(imag(Δ[1])) < 1e-3 * abs(Δ[1])
+        @test abs(imag(Δ[2])) < 1e-3 * abs(Δ[2])
+        # :ray is the default backend.
+        @test IL.GGJModel() === IL.GGJModel{:ray}()
+    end
+
+    @testset "q4 physical benchmark at Q = 500i (regime beyond :galerkin)" begin
+        q4 = IL.q4_surface_benchmark()
+        γ = 500.0im * GGJ.q0(q4)
+        Δ = IL.solve_inner(IL.GGJModel(), q4, γ)
+        # Pins from the standalone GGJRay validation (post ε_mix fix;
+        # S-invariant to 3e-4 / 7e-9 and θ-stable there).
+        @test Δ[1] ≈ 2.4720608737 + 13.354123514im rtol = 1e-4
+        @test Δ[2] ≈ 0.13749694953 + 0.74275468725im rtol = 1e-4
+
+        # Δ is an analytic invariant of the contour angle: the outward
+        # θ-check drift is a direct numerical error measurement.
+        Q = GGJ.inner_Q(q4, γ)
+        r2 = IL.solve_ray(q4, Q; θ=1.2 * angle(Q) / 4)
+        @test abs(r2.Δ[2] - Δ[2]) / abs(Δ[2]) < 1e-5
+        @test abs(r2.Δ[1] - Δ[1]) / abs(Δ[1]) < 1e-3
+    end
+end
+
+@testset "InnerLayer GGJ :ray internal machinery" begin
+    p = IL.q4_surface_benchmark()
+
+    @testset "cheblobatto nodes and differentiation matrix" begin
+        t, D = GGJ.cheblobatto(8)
+        @test length(t) == 9
+        @test issorted(t)                         # ascending, per the reflected convention
+        @test t[1] ≈ -1 && t[end] ≈ 1
+        @test D * ones(9) ≈ zeros(9) atol = 1e-10  # d/dt of a constant is zero
+        @test D * t ≈ ones(9) atol = 1e-10         # d/dt of the linear function is one
+    end
+
+    @testset "ode_matrix: ordinary point and type-generic build" begin
+        Q = 5.0im
+        # x = 0 is an ordinary point: the coefficient matrix is finite there.
+        M0 = GGJ.ode_matrix(p, Q, 0.0)
+        @test all(isfinite, M0)
+        @test M0[1, 4] == 1 && M0[2, 5] == 1 && M0[3, 6] == 1   # v' = (Ψ',Ξ',Υ') block
+        # The extended-precision build agrees with the Float64 build.
+        Md = GGJ.ode_matrix(Complex{GGJ.Double64}, p, Q, 0.3)
+        @test ComplexF64.(Md) ≈ GGJ.ode_matrix(p, Q, 0.3) rtol = 1e-12
+    end
+
+    @testset "parity_rows match the deltac boundary convention" begin
+        @test GGJ.parity_rows(1) == [4, 2, 3]     # odd:  Ψ'(0)=Ξ(0)=Υ(0)=0
+        @test GGJ.parity_rows(2) == [1, 5, 6]     # even: Ψ(0)=Ξ'(0)=Υ'(0)=0
+    end
+
+    @testset "decaying_pair is an orthonormal 6×2 frame" begin
+        Q = 5.0im
+        θ = angle(Q) / 4
+        E = GGJ.decaying_pair(p, Q, θ, 60.0)
+        @test size(E) == (6, 2)
+        @test all(isfinite, E)
+        @test E' * E ≈ I(2) atol = 1e-10          # columns orthonormal
+    end
+
+    @testset "profile diagnostics reconstruct finite fields" begin
+        res = IL.solve_ray(p, 5.0im)
+        prof = IL.solution_profile(res; npc=6)
+        @test size(prof.Ψ, 2) == 2 && length(prof.s) == size(prof.Ψ, 1)
+        @test all(isfinite, prof.Ψ) && all(isfinite, prof.Ξ) && all(isfinite, prof.Υ)
+        # Analytic tail evaluates on the trusted series radius.
+        asy = IL.asymptotic_profile(p, res, [res.S, 2 * res.S])
+        @test all(isfinite, asy.Ψ)
+    end
+
+    @testset "delta_convergence: small spread, consistent with solve_inner" begin
+        Q = 5.0im
+        conv = IL.delta_convergence(p, Q; verbose=false)
+        Δ = IL.solve_inner(IL.GGJModel(; solver=:ray), p, Q * GGJ.q0(p))
+        @test conv.Δ ≈ Δ rtol = 1e-6              # baseline == the plain solve
+        @test maximum(conv.spread) < 1e-4         # honest error bar is small here
+    end
+end