d-burg · Nov 7, 2026 · Nov 7, 2026 · Nov 7, 2026 · Nov 7, 2026 · Nov 7, 2026
diff --git a/docs/development/references.md b/docs/development/references.md
@@ -34,6 +34,25 @@ The ForceFreeStates module (ideal MHD stability analysis) implements methods fro
   - Published: Physics of Plasmas **25**, 032507 (2018)
   - Describes: Riccati method for ideal MHD eigenvalue problem
 
+- **Arnold (1967)**: "On a characteristic class entering into conditions of quantization"
+  - Published: Functional Analysis and Its Applications **1**, 1 (1967)
+  - Describes: the Maslov index of a path of Lagrangian subspaces, counted through the eigenvalues of a unitary (Cayley) image; basis of the exact conjugate-point count `nzero_conj` in `FixedBoundaryStability.jl`
+
+- **Robbin & Salamon (1993)**: "The Maslov index for paths"
+  - Published: Topology **32**, 827 (1993)
+  - Describes: the Maslov index as a signed count of crossings of a Lagrangian path with a reference subspace, which is what the Cayley eigenphases crossing π compute step by step
+
+- **Wilson, Snyder, Huysmans & Miller (2002)**: "Numerical studies of edge localized instabilities in tokamaks"
+  - Published: Physics of Plasmas **9**, 1277 (2002)
+  - Describes: the ELITE formulation; Eq. 9 is the perpendicular (incompressible) inertia that `PerpendicularInertia` implements; the growth rate is found by shooting with a secant iteration on an edge matching eigenvalue
+
+- **Snyder, Wilson et al. (2002)**: "Edge localized modes and the pedestal: A model based on coupled peeling–ballooning modes"
+  - Published: Physics of Plasmas **9**, 2037 (2002)
+  - Describes: ELITE's peeling–ballooning application and the γ > ω_*pi/2 diamagnetic-stabilization test that `IdealGrowthRateResult.omega_star_pi` serves
+
+- **Lunniss (2016)**: "Arbitrary-n ELITE" (PhD thesis, University of York; etheses.whiterose.ac.uk/16638)
+  - Describes: ELITE's radial finite-difference sweep (App. A.1), the edge matching matrix T(γ²) and the fictitious eigenvalue λ(γ²) driven to zero (§2.5.8), and the arbitrary-n inertia options; the comparison point for the count-based bracket of `growth_rate`, which ELITE does not have
+
 ## PerturbedEquilibrium Module
 
 The PerturbedEquilibrium module implements GPEC-style perturbed equilibrium calculations from:

diff --git a/docs/src/api.md b/docs/src/api.md
@@ -21,20 +21,24 @@ foreach(r -> save_to_hdf5("gpec.h5", r), (pe, kf, ef, tear))
 ```
 
 `solve(eq, alg; kwargs...)` is sugar for the problem form `solve(ForceBalanceProblem(eq; kwargs...), alg)`: the problem names what is solved (the force-balance equation on this equilibrium, with this mode range, wall and closure) and the algorithm names how.
-`main` translates each deck section into the problem, algorithm or control it names, calls these functions for every section present, saves each result into `[Run] output_file`, and returns `(; equil, local_stability, ffs, match, resp, pe, kinetic, error_fields, tearing)`, each `nothing` when its section is absent.
+`main` translates each deck section into the problem, algorithm or control it names, calls these functions for every section present, saves each result into `[Run] output_file`, and returns `(; equil, local_stability, ffs, match, growth_rate, resp, pe, kinetic, error_fields, tearing)`, each `nothing` when its section is absent.
 
 | Stage | Deck section | Call | Result |
 |---|---|---|---|
 | Equilibrium | `[Equilibrium]` (or `[Solovev]`, `[LargeAspectRatio]`, `[TJAnalytic]`) | `PlasmaEquilibrium(file; kwargs...)`, `PlasmaEquilibrium(config)`, `PlasmaEquilibrium(dd)` | `PlasmaEquilibrium` |
 | ForceFreeStates | `[ForceFreeStates]` + `[Forward]`, `[Riccati]` or `[Galerkin]` | `solve(ForceBalanceProblem(eq; ...), Forward() \| Riccati() \| Galerkin())` | `ForceFreeStatesResult` |
 | LocalStability | `[LocalStability]` + `[Ballooning]` | `solve(LocalStabilityProblem(eq; ...), Ballooning(; ...))` | `LocalStabilityResult` |
 | Match | `[Match]` + `[GGJ]` | `solve(MatchProblem(ffs; ...), GGJModel())` | `ForceFreeStatesResult` with `closure == :matched` |
+| GrowthRate | `[GrowthRate]` | `growth_rate(ForceBalanceProblem(eq; ...), Forward(), GrowthRateConfig(; ...); marginal=ffs)` | `IdealGrowthRateResult` |
 | Drive | `[ForcingTerms]` | `RMPField(ft_ctrl)` (a `ForcingTermsControl`), `RMPField(file)`, `RMPField(modes)` | `RMPField` |
 | PerturbedEquilibrium | `[PerturbedEquilibrium]` | `plasma_response(ffs; ...)`, then `perturbed_equilibrium(resp, rmp; ...)` | `PlasmaResponse`, `PerturbedEquilibriumState` |
 | KineticForces | `[KineticForces]` | `kinetic_torque(pe, KineticForcesControl(; ...))` | `KineticForcesResult` |
 | ErrorFields | `[ErrorFields]` | `error_fields(resp, pe.coil_sets, CoilConfig(...); ...)` | `ErrorFieldsResult` |
 | Tearing | `[Tearing]` + `[SLAYER]` or `[GGJ]` | `solve(TearingProblem(ffs; ...), SLAYERModel() \| GGJModel(); config=TearingConfig(...))` | `TearingResult` |
 
+`[GrowthRate]` runs the ideal growth-rate solve after the force-free-states stage; its keys are the `GrowthRateConfig` fields, and its `rho`, left at NaN, derives the mass density from the equilibrium's kinetic profiles.
+It needs `[ForceFreeStates] vac_flag = true` and `integrator = "Forward"`, and cannot share a deck with `[Match]`, since its marginal solve is the plain ideal one.
+
 The [Workflow](workflow.md) page describes what each stage computes; the module pages document the keywords, which are the deck keys of the matching section.
 
 ## Inspecting a result

diff --git a/docs/src/citations.md b/docs/src/citations.md
@@ -68,6 +68,22 @@ Reformulates the DCON eigenvalue problem as a Riccati matrix ODE, enabling paral
 
 ---
 
+> P. B. Snyder, R. J. Groebner, A. W. Leonard, T. H. Osborne, and H. R. Wilson, "Development and validation of a predictive model for the pedestal height,"
+> *Physics of Plasmas* **16**, 056118 (2009).
+> DOI: [10.1063/1.3122146](https://doi.org/10.1063/1.3122146)
+
+The EPED1 pedestal model. Source of the diamagnetic threshold γ > ω_*pi/2 with an effective ω_*pi linear in n, against which the ideal growth rate from `growth_rate` is compared (see the thresholds section of the stability page).
+
+---
+
+> P. B. Snyder, R. J. Groebner, J. W. Hughes, T. H. Osborne, M. Beurskens, A. W. Leonard, H. R. Wilson, and X. Q. Xu, "A first-principles predictive model of the pedestal height and width: development, testing and ITER optimization with the EPED model,"
+> *Nuclear Fusion* **51**, 103016 (2011).
+> DOI: [10.1088/0029-5515/51/10/103016](https://doi.org/10.1088/0029-5515/51/10/103016)
+
+The EPED1.6 model. Replaces the linear-in-n diamagnetic threshold with a two-segment fit to two-fluid calculations that rolls over at high n.
+
+---
+
 ## PerturbedEquilibrium
 
 > J.-K. Park, A. H. Boozer, and A. H. Glasser, "Computation of three-dimensional tokamak and spherical torus equilibria,"

diff --git a/docs/src/stability.md b/docs/src/stability.md
@@ -141,6 +141,38 @@ The perpendicular diffusivities that set the layer widths are described with the
 - Both integrators use the adaptive ninth-order `Vern9` solver at `eulerlagrange_tolerance`; the Riccati basis is recovered by renormalizing whenever the solution exceeds `renorm_threshold`.
 - The two crossing conventions differ slightly, so Riccati and Forward free-boundary energies agree to within a few percent, not to machine precision.
 
+### Counting fixed-boundary modes
+
+Two counts of the conjugate points (fixed-boundary modes fitting inside the ψ interval) come out of the forward integration.
+
+- `nzero` is DCON's criterion: a sign change of the smallest-magnitude eigenvalue of W⁻¹ = U₁U₂⁻¹ between stored steps, confirmed at the interpolated zero. It is cheap and is what the `crit` scan in the output records, but it misses a crossing whenever a stiff harmonic holds a tiny eigenvalue of W⁻¹ near zero at every stored step, and it can report crossings where the frame has lost its Hermitian structure.
+- `nzero_conj` is the Maslov index of the solution frame (Arnold 1967; Robbin & Salamon 1993), computed at every accepted step and exact for a Lagrangian frame. The Cayley transform Ψ = (U₁ − i s U₂)(U₁ + i s U₂)⁻¹ of the column-normalized frame is unitary when U₁†U₂ is Hermitian, has an eigenvalue −1 exactly where det U₁ = 0 (a conjugate point) and +1 where U₂ is singular (a pole of W⁻¹, harmless), and is invariant under U → U·T, so the Gaussian reductions do not disturb it. Each eigenphase is followed through its eigenvector from one step to the next and counts one conjugate point when its shortest arc crosses π. The phases are seeded one accepted step into each chunk, because at the axis start and after an ideal crossing U₁ has vanishing columns whose phases sit exactly at π and belong to the interval's end, not its interior; a step in which a phase rotates by more than π/4 is bisected on the integrator's dense output until the arcs are short, since an eigenvalue of W⁻¹ can pass zero within one step near a conjugate point. The scale s is set from the frame one step into the chunk and never changed within it: a change of s moves every phase but none through π, yet a large s squashes the phases toward π until a pole of W⁻¹ becomes indistinguishable from a zero crossing, so rescaling is not a substitute for bisection.
+
+`nconj_flags` counts the steps where the count cannot be trusted: a rotation above π/2 after bisection, or a unitarity defect max||λ|−1| of Ψ above 1e-3, which means the integrated frame is no longer Lagrangian to working precision (the Solovev regression fixture shows this after its q = 2 crossing, where `nzero` also misfires). `maslov_defect` records the largest defect seen. The count costs two matrix inversions and one eigendecomposition per accepted step, roughly 10–15 % of a solve, so it is off by default; `growth_rate` switches it on for all of its solves and uses it, with the sign of the least free-boundary energy, as the number of modes faster than the trial γ, which is what makes its bracket monotone. Checked on the TJ analytic equilibrium at `qa = 2.9, p0 = 0.8`, where a fixed-boundary mode at γ² ≈ 1.2e8 is invisible to `nzero`: `nzero_conj` is 1 below that eigenvalue and 0 above, and the Morse index `nzero_conj` + (negative free-boundary eigenvalues) is non-increasing in γ² and continuous across the pole, as it must be.
+
+### Diamagnetic stabilization thresholds
+
+`growth_rate` returns the ideal-MHD γ and the peak ion diamagnetic frequency `omega_star_pi` (the maximum of |ω_*pi| over the integration domain, at the mode's n). It applies no stability threshold. Pedestal models call a peeling–ballooning mode unstable only when γ exceeds half of an effective diamagnetic frequency, and the choice of that frequency moves the stability boundary, so every boundary quoted from these numbers must name the convention it used.
+
+| Convention | Unstable when | Ratio in terms of the outputs |
+|:--|:--|:--|
+| Ideal | γ > 0 | none |
+| EPED1 ([Snyder (2009)](citations.md#ForceFreeStates)) | γ > ω_eff/2 with ω_eff half the pedestal maximum of ω_*pi, so ω_eff ∝ n | R₁ = `gamma / (omega_star_pi/4)` > 1 |
+| EPED1.6 ([Snyder (2011)](citations.md#ForceFreeStates)) | γ > ω_eff/2 with ω_eff from a two-segment fit to two-fluid calculations: ∝ n up to n_crit, a weaker slope above | R₁.₆ = R₁/√2 for n ≤ n_crit, and R₁/(√2 (n_crit/n + 0.168 (n − n_crit)/n)) above; unstable when > 1 |
+
+R₁ is the quantity ELITE prints as γ/(ω_*pi,max/4). The EPED1.6 row is written in the form used by the stability-map post-processing of the OMFIT ELITE module, which defaults to n_crit = 6 and is what the conventional pedestal stability diagrams (pressure gradient α against edge current, contour at 1, maximum over n) show.
+
+The motivation for the second fit is that ω_*pi varies strongly across a pedestal while a high-n mode is radially narrow, so the frequency the mode feels grows more slowly than n; a threshold linear in n over-stabilizes high n.
+
+Points to keep in view when using either ratio:
+
+- The two ratios are equal at n ≈ 1.54 n_crit (n ≈ 9 for n_crit = 6). Below that the EPED1.6 form is the more stable one, by √2 at n ≤ n_crit; above it the less stable, by a factor 1.92 at n = 25 and 2.64 at n = 50. The boundary of a map taken as the maximum over n therefore moves, and so does the most-unstable n.
+- With the EPED1.6 form the maximum over n is often set by the highest n computed, so the n range is part of the result.
+- n_crit is a model parameter, not a constant: the value 6 is a default, and a published DIII-D analysis scales it with the edge safety factor (n_crit ≈ 21.7/q at the peak pressure gradient).
+- The factor √2 and the slope 0.168 are quoted from that post-processing; check them against the fit in [Snyder (2011)](citations.md#ForceFreeStates) before relying on them outside a like-for-like comparison.
+- `omega_star_pi` comes from the attached kinetic profiles with n_i ∝ n_e and Z_i = 1, and its peak is taken over `[psilow, psilim]`; a truncated domain can cut off part of the pedestal gradient.
+- A comparison of two codes is independent of the convention as long as both growth rates go through the same one; the boundary location is not.
+
 Local stability (Mercier, resistive interchange and ballooning) is a separate stage on the same equilibrium; see [Local Stability](ballooning.md).
 
 ## API reference
@@ -149,5 +181,5 @@ The Galerkin solver's functions are documented on the [Galerkin Solver](galerkin
 
 ```@autodocs
 Modules = [GeneralizedPerturbedEquilibrium.ForceFreeStates]
-Pages = ["ForceFreeStates.jl", "CoreTypes.jl", "Surfaces/Types.jl", "Riccati/Types.jl", "Matching/DeltaPrime.jl", "Result.jl", "Surfaces/ResistEval.jl", "Matching/ResonantMatch.jl", "Matching/MatchProblem.jl", "EulerLagrange.jl", "Surfaces/Finding.jl", "Surfaces/Asymptotics.jl", "Fourfit.jl", "Kinetic.jl", "FixedKineticMatrices.jl", "FixedBoundaryStability.jl", "Utils.jl", "Free.jl", "Riccati/Propagators.jl", "Riccati/Crossings.jl", "Riccati/DeltaPrimeBVP.jl", "Riccati/Driver.jl"]
+Pages = ["ForceFreeStates.jl", "CoreTypes.jl", "Surfaces/Types.jl", "Riccati/Types.jl", "Matching/DeltaPrime.jl", "Result.jl", "Surfaces/ResistEval.jl", "Matching/ResonantMatch.jl", "Matching/MatchProblem.jl", "EulerLagrange.jl", "Surfaces/Finding.jl", "Surfaces/Asymptotics.jl", "Fourfit.jl", "Kinetic.jl", "Inertia.jl", "FixedKineticMatrices.jl", "FixedBoundaryStability.jl", "Utils.jl", "Free.jl", "Riccati/Propagators.jl", "Riccati/Crossings.jl", "Riccati/DeltaPrimeBVP.jl", "Riccati/Driver.jl"]
 ```
diff --git a/docs/src/workflow.md b/docs/src/workflow.md
@@ -21,6 +21,7 @@ Every stage after `Equilibrium` reads the force-free states, and `KineticForces`
 |---|---|---|
 | Is the equilibrium read and mapped correctly; what are q and β_N? | `[Equilibrium]` (+ `[ForceFreeStates]`, which triggers the write) | `Equilibrium/q_95`, `Equilibrium/beta_N`, `Equilibrium/Profiles/` |
 | Is the plasma ideal-kink stable, with or without a wall? | `[ForceFreeStates]` (+ `[Wall]`) | `ForceFreeStates/FreeBoundaryStability/W_freeboundary`, `…/eigenmode_energies` |
+| How fast does the least-stable free-boundary ideal mode grow? | `[GrowthRate]` (+ `[ForceFreeStates] vac_flag = true`, Forward) | `ForceFreeStates/GrowthRate/` |
 | Is it Mercier, resistive-interchange or ballooning unstable? | `[LocalStability]` (+ `[Ballooning]`) | `LocalStability/D_I`, `LocalStability/D_R`, `LocalStability/ballooning_Delta_prime` |
 | What is the tearing stability index Δ′ at each rational surface? | `[ForceFreeStates] integrator = "Riccati"` or `"Galerkin"` (+ `[Riccati]` or `[Galerkin]`) | `SingularSurfaces/Delta_prime_matrix` |
 | How does the plasma respond to an RMP or error field; do islands overlap? | `[ForcingTerms]` + `[PerturbedEquilibrium]` | `PerturbedEquilibrium/SingularCoupling/resonant_area_weighted_field`, `…/island_half_width`, `…/chirikov_parameter` |

diff --git a/regression-harness/cases/lar_reversed_shear_growth_rate_n1.toml b/regression-harness/cases/lar_reversed_shear_growth_rate_n1.toml
@@ -0,0 +1,79 @@
+# Regression case: ELITE-style ideal growth rate of the large-aspect-ratio reversed-shear n=1 example.
+# The example is free-boundary unstable at n=1, so [GrowthRate] root-finds a nonzero γ with a uniform
+# deuterium mass density (1e20 m⁻³). Tracks γ, γ², the root-find history and the ideal energy it starts from.
+# Each [quantities.*] block names an HDF5 path, how to extract it, and its noise floor.
+[case]
+name = "lar_reversed_shear_growth_rate_n1"
+description = "Large-aspect-ratio reversed-shear equilibrium, n=1, ELITE-style ideal growth rate with uniform deuterium density"
+example_dir = "examples/LAR_reversed_shear_example"
+
+[overrides]
+"GrowthRate.rho" = 3.3e-7                # Uniform mass density [kg/m³] (deuterium at 1e20 m⁻³); NaN derives ρ(ψ) from the kinetic profiles
+
+# Growth rate (γ ≈ 3.0e3 1/s); absolute floors match a 1e-10 shift in λ_min near the root
+[quantities.gamma]
+h5path = "ForceFreeStates/GrowthRate/gamma"
+type = "real_scalar"
+extract = "value"
+label = "growth rate gamma (1/s)"
+noise_threshold = 1e-5
+order = 10
+
+[quantities.gamma_sq]
+h5path = "ForceFreeStates/GrowthRate/gamma_sq"
+type = "real_scalar"
+extract = "value"
+label = "growth rate squared gamma^2 (1/s^2)"
+noise_threshold = 1e-1
+order = 11
+
+# Root-find history: every γ² evaluated, the least-stable free-boundary energy there, and the fixed-boundary zero count
+[quantities.history_gamma_sq]
+h5path = "ForceFreeStates/GrowthRate/History/gamma_sq"
+type = "real_vector"
+extract = "all_real"
+label = "root-find gamma^2 history"
+noise_threshold = 1e-1
+order = 20
+
+[quantities.history_lambda_min]
+h5path = "ForceFreeStates/GrowthRate/History/lambda_min"
+type = "real_vector"
+extract = "all_real"
+label = "root-find lambda_min history"
+noise_threshold = 1e-10
+order = 21
+
+[quantities.history_nzero]
+h5path = "ForceFreeStates/GrowthRate/History/nzero"
+type = "real_vector"
+extract = "all_real"
+label = "root-find nzero history (integer counts)"
+noise_threshold = 0
+order = 22
+
+# Ideal marginal solve the root find starts from
+[quantities.et_real]
+h5path = "ForceFreeStates/FreeBoundaryStability/eigenmode_energies"
+type = "complex_vector"
+extract = "real_first"
+label = "total energy Re(et[1])"
+noise_threshold = 1e-10
+order = 30
+
+[quantities.msing]
+h5path = "SingularSurfaces/rational_count"
+type = "int_scalar"
+extract = "value"
+label = "# singular surfaces"
+noise_threshold = 0
+order = 31
+
+# Runtime (special: not from H5)
+[quantities.runtime]
+h5path = ""
+type = "runtime"
+extract = "value"
+label = "Runtime (s)"
+noise_threshold = 0.0
+order = 999