diff --git a/docs/development/references.md b/docs/development/references.md
index 9dbb3c411..4aa06c67a 100644
--- 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
index 06a117a30..51a401514 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -21,7 +21,7 @@ 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 |
 |---|---|---|---|
@@ -29,12 +29,16 @@ foreach(r -> save_to_hdf5("gpec.h5", r), (pe, kf, ef, tear))
 | 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/stability.md b/docs/src/stability.md
index aa86df2ec..74875de15 100644
--- a/docs/src/stability.md
+++ b/docs/src/stability.md
@@ -141,6 +141,15 @@ 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.
+
 Local stability (Mercier, resistive interchange and ballooning) is a separate stage on the same equilibrium; see [Local Stability](ballooning.md).
 
 ## API reference
@@ -149,5 +158,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
index 38b4cdd8b..9a280c199 100644
--- 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
new file mode 100644
index 000000000..d3351bbb5
--- /dev/null
+++ 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
diff --git a/src/Deck.jl b/src/Deck.jl
index ea0f4b17c..1cd61eb50 100644
--- a/src/Deck.jl
+++ b/src/Deck.jl
@@ -78,7 +78,7 @@ const _RETIRED_SECTIONS = Dict("DEBUG" => "[DEBUG] is retired: output_benchmark_
 
 # Every section a deck may hold, besides the analytic-equilibrium sections.
 const _DECK_SECTIONS = ("Run", "Equilibrium", "Wall", "LocalStability", "Ballooning", "ForceFreeStates", "Forward", "Riccati",
-    "Galerkin", "Match", "GGJ", "Tearing", "SLAYER", "ForcingTerms", "PerturbedEquilibrium", "KineticForces", "ErrorFields")
+    "Galerkin", "Match", "GGJ", "Tearing", "SLAYER", "GrowthRate", "ForcingTerms", "PerturbedEquilibrium", "KineticForces", "ErrorFields")
 const _INTEGRATORS = ("Forward" => Forward, "Riccati" => Riccati, "Galerkin" => Galerkin)
 
 # Stop the run when a parsed gpec.toml section sets keys that no longer exist, so the deck is updated.
@@ -377,6 +377,16 @@ function tearing_from_deck(inputs::AbstractDict)
     return (; problem, model=inner_model == "SLAYER" ? slayer_model(inputs) : ggj_model(inputs), config)
 end
 
+"""
+    growth_rate_from_deck(inputs) -> Union{Nothing,GrowthRateConfig}
+
+The `GrowthRateConfig` of `[GrowthRate]`; `nothing` without that section.
+"""
+function growth_rate_from_deck(inputs::AbstractDict)
+    haskey(inputs, "GrowthRate") || return nothing
+    return GrowthRateConfig(; section_keywords(inputs["GrowthRate"], "GrowthRate", _struct_fields(GrowthRateConfig))...)
+end
+
 """
     perturbed_equilibrium_from_deck(inputs) -> Union{Nothing,NamedTuple}
 
@@ -447,7 +457,7 @@ function parse_deck(inputs::AbstractDict, path::AbstractString)
     end
     has_ffs = haskey(inputs, "ForceFreeStates")
     if !has_ffs
-        for section in ("Wall", "Forward", "Riccati", "Galerkin", "Match", "GGJ", "Tearing", "SLAYER", "LocalStability", "Ballooning",
+        for section in ("Wall", "Forward", "Riccati", "Galerkin", "Match", "GGJ", "Tearing", "SLAYER", "GrowthRate", "LocalStability", "Ballooning",
             "ForcingTerms", "PerturbedEquilibrium", "KineticForces", "ErrorFields")
             haskey(inputs, section) && error("[$section] is set but there is no [ForceFreeStates] section, so the run stops after the equilibrium")
         end
@@ -462,6 +472,15 @@ function parse_deck(inputs::AbstractDict, path::AbstractString)
         error("[KineticForces] is set but neither [PerturbedEquilibrium] (the torque) nor [ForceFreeStates] kinetic = \"KineticForces\" uses it")
     (haskey(inputs, "KineticForces") || kinetic_source) && isempty(get(get(inputs, "Equilibrium", Dict{String,Any}()), "kinetic_file", "")) &&
         error("this run needs kinetic profiles ([KineticForces], or kinetic = \"KineticForces\"): set kinetic_file in [Equilibrium]")
+    growth_rate = growth_rate_from_deck(inputs)
+    # The growth-rate solve needs the plain ideal solve as its marginal one and the free-boundary energy of a Forward solve.
+    if growth_rate !== nothing
+        (; keywords, algorithm) = force_balance
+        haskey(inputs, "Match") && error("[GrowthRate] needs the ideal solve as its marginal one: remove the [Match] section")
+        get(keywords, :kinetic, nothing) === nothing || error("[GrowthRate] needs an ideal force-balance solve: remove kinetic from [ForceFreeStates]")
+        get(keywords, :vac_flag, false) || error("[GrowthRate] needs the free-boundary energy: set [ForceFreeStates] vac_flag = true")
+        algorithm isa Forward || error("[GrowthRate] needs integrator = \"Forward\" in [ForceFreeStates]")
+    end
     match = match_from_deck(inputs)
     # The MatchProblem guards, checked before the force-free-states solve rather than after it.
     if match !== nothing
@@ -483,7 +502,7 @@ function parse_deck(inputs::AbstractDict, path::AbstractString)
         error("[GGJ] is set but neither [Match] nor [Tearing] inner_model = \"GGJ\" uses it")
     haskey(inputs, "SLAYER") && tearing_model != "SLAYER" && !(has_ffs && get(inputs["ForceFreeStates"], "psilim_from_layer_overlap", false) === true) &&
         error("[SLAYER] is set but neither [Tearing] inner_model = \"SLAYER\" nor [ForceFreeStates] psilim_from_layer_overlap = true uses it")
-    return (; run=run_settings(inputs, path), force_balance, local_stability=local_stability_from_deck(inputs), match, tearing,
+    return (; run=run_settings(inputs, path), force_balance, local_stability=local_stability_from_deck(inputs), match, tearing, growth_rate,
         kinetic_forces=haskey(inputs, "KineticForces") ? kinetic_forces_control(inputs) : nothing, perturbed_equilibrium, forcing,
         error_fields=error_fields_from_deck(inputs, path))
 end
@@ -578,6 +597,7 @@ function deck_inputs(deck::NamedTuple, config::Equilibrium.EquilibriumConfig, an
         merge!(inputs, _tearing_sections(filter(p -> p.first != :critical_resonant_field, problem), get(problem, :critical_resonant_field, nothing),
             model, config))
     end
+    deck.growth_rate === nothing || (inputs["GrowthRate"] = _deck_table(deck.growth_rate))
     deck.kinetic_forces === nothing || merge!(inputs, deck_sections(deck.kinetic_forces))
     pe = deck.perturbed_equilibrium
     pe === nothing || (inputs["PerturbedEquilibrium"] = _deck_table(merge(pe.response, pe.drive)))
@@ -625,6 +645,7 @@ deck_sections(ls::LocalStabilityResult) = Dict{String,Any}("LocalStability" => _
 deck_sections(m::ForceFreeStates.MatchResult) =
     merge(deck_sections(m.model), Dict{String,Any}("Match" => merge(_deck_table(m.problem, (:ffs,)), Dict{String,Any}("inner_model" => "GGJ"))))
 deck_sections(kf::KineticForcesResult) = deck_sections(kf.control)
+deck_sections(g::IdealGrowthRateResult) = g.config === nothing ? Dict{String,Any}() : Dict{String,Any}("GrowthRate" => _deck_table(g.config))
 deck_sections(t::TearingResult) = t.model === nothing ? Dict{String,Any}() : _tearing_sections(t.problem, t.problem.critical_resonant_field, t.model, t.config)
 deck_sections(ef::ErrorFieldsResult) = _error_fields_sections(ef.settings, ef.tolerances === nothing ? "" : _RECORDED_TOLERANCES)
 
diff --git a/src/ForceFreeStates/CoreTypes.jl b/src/ForceFreeStates/CoreTypes.jl
index 1cb794f65..89eedba73 100644
--- a/src/ForceFreeStates/CoreTypes.jl
+++ b/src/ForceFreeStates/CoreTypes.jl
@@ -37,6 +37,8 @@ Solve-time scratch state for a force-balance solve: the resolved mode space, sin
   - `q1lim::Float64` - Safety factor derivative at psilim
   - `wall_settings::Union{Nothing,Vacuum.WallSettings}` - Wall shape settings for vacuum calculations
   - `verbose::Bool` - Log solve progress; set from the `verbose` keyword of `solve`
+  - `count_conjugate_points::Bool` - Run the exact conjugate-point count of the forward integration (`OdeState.nzero_conj`); set by
+    `growth_rate` for its solves, off otherwise
 """
 @kwdef mutable struct ForceFreeStatesInternal <: ModeSpace
     mlow::Int = 0
@@ -59,6 +61,7 @@ Solve-time scratch state for a force-balance solve: the resolved mode space, sin
     q1lim::Float64 = 0.0
     wall_settings::Union{Nothing,Vacuum.WallSettings} = nothing
     verbose::Bool = false
+    count_conjugate_points::Bool = false
     """
     Inter-surface Δ' matrix, `msing × msing`, PEST3 convention: the odd-parity projection of
     `delta_prime_raw`, set by `compute_delta_prime_matrix!` (Riccati path only).
diff --git a/src/ForceFreeStates/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl
index d40700059..52ad96a6e 100644
--- a/src/ForceFreeStates/EulerLagrange.jl
+++ b/src/ForceFreeStates/EulerLagrange.jl
@@ -67,6 +67,19 @@ records needed to undo it, and singular-surface asymptotic coefficients.
   - `ising_start::Int` - Index of the last singular surface inside the starting ψ.
   - `needs_crossing::Bool` - Whether a rational surface must be crossed after the current region.
   - `nzero::Int` - Number of detected zero crossings of the critical eigenvalue.
+  - `count_conj::Bool` - Whether the exact conjugate-point count runs (forward Hamiltonian runs that request it: no kinetic source, or a
+    `PerpendicularInertia` one); off by default since it costs two LU factorizations and one eigenvalue solve per accepted step.
+  - `nconj::Int` - Signed running count of Cayley eigenphases that wrapped through π (conjugate points, fixed-boundary modes inside the
+    ψ interval), reseeded after each ideal crossing so it counts per Newcomb interval.
+  - `nzero_conj::Int` - Number of conjugate points, `abs(nconj)` at the end of the run; −1 when the count did not run.
+  - `nconj_flags::Int` - Accepted steps whose count is unreliable (an eigenphase rotation above π/2, or a unitarity defect above 1e-3 that
+    shows the frame is not Lagrangian to working precision); a nonzero value means `nzero_conj` may be wrong.
+  - `maslov_theta::Vector{Float64}` - Cayley eigenphases at the last counted step, aligned with the columns of `maslov_vecs`.
+  - `maslov_vecs::Matrix{ComplexF64}` - Cayley eigenvectors at the last counted step, which identify each eigenphase across a step.
+  - `maslov_pending::Bool` - Whether the next accepted step seeds the eigenphases and the Cayley scale instead of advancing the count.
+  - `maslov_s::Float64` - Scale s of the Cayley transform (U₁ − i s U₂)(U₁ + i s U₂)⁻¹, set to ‖U₁‖/‖U₂‖ of the column-normalized frame one step
+    into each chunk.
+  - `maslov_defect::Float64` - Largest unitarity defect max| |λ| − 1 | of the Cayley eigenvalues over the run.
   - `new::Bool` - Whether the next norm call records fresh reference norms `unorm0`.
   - `unorm::Vector{Float64}` - Current solution-column norms.
   - `unorm0::Vector{Float64}` - Reference norms since the last reduction.
@@ -111,6 +124,17 @@ records needed to undo it, and singular-surface asymptotic coefficients.
     needs_crossing::Bool = false
     nzero::Int = 0
 
+    # Exact conjugate-point count
+    count_conj::Bool = false
+    nconj::Int = 0
+    nzero_conj::Int = -1
+    nconj_flags::Int = 0
+    maslov_theta::Vector{Float64} = zeros(Float64, numpert_total)
+    maslov_vecs::Matrix{ComplexF64} = zeros(ComplexF64, numpert_total, numpert_total)
+    maslov_pending::Bool = false
+    maslov_s::Float64 = 1.0
+    maslov_defect::Float64 = 0.0
+
     # Used for Gaussian reduction
     new::Bool = true
     unorm::Vector{Float64} = zeros(Float64, numpert_total)
@@ -285,6 +309,7 @@ function forward_eulerlagrange_integration(ctrl::ForceFreeStatesControl, equil::
 
     odet = OdeState(intr.numpert_total, ctrl.algorithm.numsteps_init, ctrl.algorithm.numunorms_init, intr.msing)
     initialize_el_at_axis!(odet, ctrl, mats, equil.profiles, intr)
+    odet.count_conj = intr.count_conjugate_points && kinetic_source(ctrl.kinetic) isa Union{Nothing,PerpendicularInertia}
 
     chunks = chunk_el_integration_bounds(odet, ctrl, intr)
 
@@ -293,6 +318,8 @@ function forward_eulerlagrange_integration(ctrl::ForceFreeStatesControl, equil::
     end
 
     for chunk in chunks
+        # Reseeds the count one step into the chunk: at the Frobenius start and after an ideal crossing U₁ has vanishing columns
+        odet.count_conj && maslov_seed!(odet)
         integrate_el_region!(odet, ctrl, equil, mats, intr, chunk)
         if verbose
             @info "   ψ = $((@sprintf "%.3f" odet.psifac)),  q = $((@sprintf "%.3f" odet.q)),  steps = $(odet.total_steps)"
@@ -318,6 +345,7 @@ function forward_eulerlagrange_integration(ctrl::ForceFreeStatesControl, equil::
         @info "Evaluating fixed-boundary stability criterion"
     end
     odet.nzero = evaluate_stability_criterion!(odet, equil.profiles)
+    odet.count_conj && (odet.nzero_conj = abs(odet.nconj))
 
     # Undo Gaussian reduction to get true solution vectors
     transform_u!(odet, intr)
@@ -663,6 +691,7 @@ function cross_kinetic_singular_surf!(
     odet.psifac = ksurf.psifac + dpsi
     sing_der!(du2, odet.u, params, odet.psifac)
     odet.u .+= (du1 .+ du2) .* dpsi
+    odet.count_conj && maslov_advance!(odet, odet.u)   # compared with the phases of the chunk's last accepted step
 
     odet.q = equil.profiles.q_spline(odet.psifac; hint=odet.spline_hint)
     store_ode_data!(odet, odet.psifac, odet.u)
@@ -746,6 +775,8 @@ function integrate_el_region!(
 
         compute_solution_norms!(integrator.u, odet, ctrl, intr, false)
         column_abstol!(integrator.opts.abstol, integrator.u, integrator.opts.reltol)
+        # The dense output lets a step with a fast-moving eigenphase be bisected for the conjugate-point count
+        odet.count_conj && maslov_advance!(odet, integrator.u; refine=τ -> integrator(integrator.tprev + τ * (integrator.t - integrator.tprev)))
 
         # The first step is always saved, even when q_range ≈ 0
         near_start = abs(odet.q - q_start) < near_q_frac * q_range || steps_in_segment[] == 1
diff --git a/src/ForceFreeStates/FixedBoundaryStability.jl b/src/ForceFreeStates/FixedBoundaryStability.jl
index 29896ac3d..d557fc04e 100644
--- a/src/ForceFreeStates/FixedBoundaryStability.jl
+++ b/src/ForceFreeStates/FixedBoundaryStability.jl
@@ -106,3 +106,126 @@ smallest magnitude. Also flags whether W⁻¹ was non-Hermitian beyond tolerance
     crit = evals[findmin(abs, evals)[2]]
     return crit, nonherm
 end
+
+"""
+    cayley_eigen!(θ, V, u1, u2, s) -> defect
+
+Eigen-decompose the Cayley transform Ψ = (U₁ − i s U₂)(U₁ + i s U₂)⁻¹ of the solution frame `[U₁; U₂]`,
+writing the eigenphases in (−π, π] to `θ` and the eigenvectors to the columns of `V`, and return the
+unitarity defect max| |λ| − 1 |. For a Lagrangian frame Ψ is unitary; an eigenvalue −1 occurs exactly where
+det U₁ = 0, a conjugate point, while an eigenvalue +1 marks a singular U₂, a harmless pole of W⁻¹ (Arnold
+1967, Funct. Anal. Appl. 1, 1; Robbin & Salamon 1993, Topology 32, 827). Ψ is invariant under U → U·T for
+any invertible T, so Gaussian reduction leaves it unchanged and the frame can be passed with normalized
+columns (see `normalized_frame`). When U₁†U₂ drifts from Hermitian, Ψ is not exactly unitary; the
+eigenpairs are used regardless and the defect reports the drift.
+"""
+function cayley_eigen!(θ::AbstractVector{Float64}, V::AbstractMatrix{ComplexF64}, u1::AbstractMatrix{ComplexF64}, u2::AbstractMatrix{ComplexF64}, s::Real)
+    psi = (u1 .- (im * s) .* u2) / (u1 .+ (im * s) .* u2)
+    E = eigen!(psi)
+    defect = 0.0
+    for (j, l) in enumerate(E.values)
+        defect = max(defect, abs(abs(l) - 1))
+        a = angle(l)
+        θ[j] = a <= -π ? Float64(π) : a
+    end
+    V .= E.vectors
+    return defect
+end
+
+"""
+    normalized_frame(u) -> (u1, u2)
+
+Copies of the two halves of the solution frame `u = [U₁; U₂]` with each column of the stacked frame scaled
+to unit norm, so that the Cayley inverse is conditioned by the frame's directions and not by the column
+norms, which span many decades between Gaussian reductions.
+"""
+function normalized_frame(u::AbstractArray{ComplexF64,3})
+    u1 = copy(@view u[:, :, 1])
+    u2 = copy(@view u[:, :, 2])
+    for j in axes(u1, 2)
+        c = hypot(norm(@view u1[:, j]), norm(@view u2[:, j]))
+        c > 0 && (u1[:, j] ./= c; u2[:, j] ./= c)
+    end
+    return u1, u2
+end
+
+"""
+    maslov_match(θold, Vold, θnew, Vnew) -> (w, maxrot, overlap)
+
+Follow each eigenphase from one step to the next through the eigenvectors, assigning old to new by greatest
+overlap |⟨v_old, v_new⟩| first until every pair is matched (the eigenvectors of a unitary Ψ are orthonormal,
+so a clear step gives overlaps near 1; within a cluster of nearly equal phases the assignment is arbitrary
+but harmless, since swapped partners rotate by the same small arcs). Each phase moves along the shortest arc
+δ = angle(exp(i(θnew − θold))); return the signed number of arcs that cross π (a conjugate point, +1 upward
+and −1 downward), the largest |δ|, and the smallest overlap used.
+"""
+function maslov_match(θold::AbstractVector{Float64}, Vold::AbstractMatrix{ComplexF64}, θnew::AbstractVector{Float64}, Vnew::AbstractMatrix{ComplexF64})
+    O = abs.(Vold' * Vnew)
+    N = length(θnew)
+    partner = zeros(Int, N)
+    overlap = 1.0
+    for _ in 1:N
+        o, idx = findmax(O)
+        j, k = Tuple(idx)
+        partner[k] = j
+        overlap = min(overlap, o)
+        O[j, :] .= -1.0
+        O[:, k] .= -1.0
+    end
+    w, maxrot = 0, 0.0
+    for k in 1:N
+        d = angle(cis(θnew[k] - θold[partner[k]]))
+        maxrot = max(maxrot, abs(d))
+        target = θold[partner[k]] + d
+        w += target > π ? 1 : target <= -π ? -1 : 0
+    end
+    return w, maxrot, overlap
+end
+
+"""
+    maslov_seed!(odet)
+
+Mark the conjugate-point count for reseeding at the next accepted step. At a chunk start U₁ has vanishing
+columns (every column at the axis, the zeroed resonant column after an ideal crossing), so eigenphases sit
+exactly at π there; seeding one step later counts conjugate points in the open interval, as Newcomb's
+theorem states it, and sets the Cayley scale s = ‖U₁‖/‖U₂‖ from a nonsingular frame.
+"""
+function maslov_seed!(odet::OdeState)
+    odet.maslov_pending = true
+    return odet
+end
+
+"""
+    maslov_advance!(odet, u; refine=nothing, depth=0)
+
+Advance the conjugate-point count from the stored eigenpairs to those of frame `u`: add the signed number
+of eigenphases that crossed π to `odet.nconj`. A phase that rotated by more than π/4 in the step leaves its
+direction of travel in doubt (an eigenvalue of W⁻¹ passing zero or infinity quickly), so when `refine(τ)`
+can supply the frame at fraction τ of the step (the integrator's dense output), the step is bisected, up to
+six times, until every arc is short. The step is counted in `odet.nconj_flags` when a rotation above π/2
+remains or the unitarity defect exceeds 1e-3 (the frame is not Lagrangian to working precision, so its
+eigenphases are unreliable). A pending seed stores the eigenpairs and the scale instead.
+"""
+function maslov_advance!(odet::OdeState, u::AbstractArray{ComplexF64,3}; refine=nothing, depth::Int=0)
+    u1, u2 = normalized_frame(u)
+    if odet.maslov_pending
+        n1, n2 = norm(u1), norm(u2)
+        odet.maslov_s = n1 > 0 && n2 > 0 ? n1 / n2 : 1.0
+        odet.maslov_defect = max(odet.maslov_defect, cayley_eigen!(odet.maslov_theta, odet.maslov_vecs, u1, u2, odet.maslov_s))
+        odet.maslov_pending = false
+        return odet
+    end
+    θnew, Vnew = similar(odet.maslov_theta), similar(odet.maslov_vecs)
+    defect = cayley_eigen!(θnew, Vnew, u1, u2, odet.maslov_s)
+    w, maxrot, _ = maslov_match(odet.maslov_theta, odet.maslov_vecs, θnew, Vnew)
+    if maxrot > π / 4 && refine !== nothing && depth < 6
+        maslov_advance!(odet, refine(0.5); refine=τ -> refine(0.5 * τ), depth=depth + 1)
+        return maslov_advance!(odet, u; refine=τ -> refine(0.5 + 0.5 * τ), depth=depth + 1)
+    end
+    odet.maslov_defect = max(odet.maslov_defect, defect)
+    (maxrot > π / 2 || defect > 1e-3) && (odet.nconj_flags += 1)
+    odet.nconj += w
+    odet.maslov_theta .= θnew
+    odet.maslov_vecs .= Vnew
+    return odet
+end
diff --git a/src/ForceFreeStates/ForceBalance.jl b/src/ForceFreeStates/ForceBalance.jl
index f4db1a240..57542f7cc 100644
--- a/src/ForceFreeStates/ForceBalance.jl
+++ b/src/ForceFreeStates/ForceBalance.jl
@@ -92,9 +92,12 @@ ffs = solve(prob, Riccati())
 ffs = solve(eq, Riccati(); nn=1, vac_flag=true)   # equivalent one-line form
 ```
 """
-function CommonSolve.solve(prob::ForceBalanceProblem, alg::AbstractIntegrator; verbose::Bool=false)
+CommonSolve.solve(prob::ForceBalanceProblem, alg::AbstractIntegrator; verbose::Bool=false) = _solve_force_balance(prob, alg; verbose)
+
+# `solve` with the internal switch for the exact conjugate-point count, which growth_rate turns on for its solves.
+function _solve_force_balance(prob::ForceBalanceProblem, alg::AbstractIntegrator; verbose::Bool=false, count_conjugate_points::Bool=false)
     total_start = time()
-    ctrl, intr, equil, metric, mats, layer_overlap = _prepare_solve(prob, alg; verbose)
+    ctrl, intr, equil, metric, mats, layer_overlap = _prepare_solve(prob, alg; verbose, count_conjugate_points)
     result = run_force_free_states(ctrl, equil, mats, intr, metric; cut_solution=alg isa Galerkin && alg.cut_solution)
     result = _with(result; problem=prob, algorithm=alg, layer_overlap=layer_overlap)
 
@@ -123,7 +126,7 @@ end
 
 # Everything solve does before the integration: the control, the mode space, the (re-formed) equilibrium, the metric,
 # the matrices, the layer-overlap scan and the kinetic matrices.
-function _prepare_solve(prob::ForceBalanceProblem, alg::AbstractIntegrator; verbose::Bool=false)
+function _prepare_solve(prob::ForceBalanceProblem, alg::AbstractIntegrator; verbose::Bool=false, count_conjugate_points::Bool=false)
     equil = prob.equil
     ctrl_kwargs = copy(prob.ctrl_kwargs)
     ctrl = ForceFreeStatesControl(; algorithm=alg, ctrl_kwargs...)
@@ -133,9 +136,9 @@ function _prepare_solve(prob::ForceBalanceProblem, alg::AbstractIntegrator; verb
         error("kinetic_factor scales the kinetic matrices, but no kinetic source is set: pass kinetic=... " *
               "([ForceFreeStates] kinetic = \"KineticForces\")")
     # A kinetic model (not supplied or synthetic matrices) is evaluated on the equilibrium's kinetic profiles.
-    ctrl.kinetic === nothing || ctrl.kinetic isa Union{KineticMatrices,TestKineticMatrices} || equil.kinetic !== nothing ||
+    ctrl.kinetic === nothing || ctrl.kinetic isa Union{KineticMatrices,TestKineticMatrices,PerpendicularInertia} || equil.kinetic !== nothing ||
         error("the KineticForces kinetic source needs kinetic profiles on the equilibrium: set kinetic_file or call attach_kinetic_profiles!")
-    ctrl.kinetic isa Union{Nothing,KineticMatrices,TestKineticMatrices} || Equilibrium.has_ntv_fields(equil.kinetic) ||
+    ctrl.kinetic isa Union{Nothing,KineticMatrices,TestKineticMatrices,PerpendicularInertia} || Equilibrium.has_ntv_fields(equil.kinetic) ||
         error("the KineticForces kinetic source needs n_e, T_i, T_e and omega_E in the kinetic profiles; the kinetic file lacks at least one")
     ctrl.kinetic !== nothing && !(alg isa Forward) &&
         error("$(nameof(typeof(alg)))() does not support kinetic runs; use Forward().")
@@ -148,7 +151,7 @@ function _prepare_solve(prob::ForceBalanceProblem, alg::AbstractIntegrator; verb
               "nodes; regenerate the file in double precision for converged Δ′." file = eq_config.eq_filename
     end
 
-    intr = ForceFreeStatesInternal(; verbose)
+    intr = ForceFreeStatesInternal(; verbose, count_conjugate_points)
     intr.wall_settings = prob.wall
 
     resolve_mode_space!(intr, ctrl)
@@ -272,6 +275,9 @@ function prepare_force_free_states!(
     elseif ctrl.kinetic isa TestKineticMatrices
         intr.verbose && @info "Computing kinetic matrices (synthetic test matrices, scale $(ctrl.kinetic.scale))"
         km = evaluate_kinetic_matrices(ctrl, equil, mats, intr, metric)
+    elseif ctrl.kinetic isa PerpendicularInertia
+        intr.verbose && @info "Computing perpendicular-inertia matrices (γ² = $(ctrl.kinetic_factor) s⁻²)"
+        km = evaluate_kinetic_matrices(ctrl, equil, mats, intr, metric)
     elseif ctrl.kinetic !== nothing
         intr.verbose && @info "Computing kinetic matrices (KineticForces model)"
         calculated_cb, axis_psi_c = kinetic_model_source(ctrl.kinetic, equil, intr)
diff --git a/src/ForceFreeStates/ForceFreeStates.jl b/src/ForceFreeStates/ForceFreeStates.jl
index f1374dae3..7c08fe8ea 100644
--- a/src/ForceFreeStates/ForceFreeStates.jl
+++ b/src/ForceFreeStates/ForceFreeStates.jl
@@ -72,6 +72,8 @@ include("InputSnapshot.jl")
 include("Output.jl")
 include("Run.jl")
 include("ForceBalance.jl")
+include("Inertia.jl")
+include("InertiaOutput.jl")
 
 # Small tolerance and root-finder iteration cap used throughout ForceFreeStates
 global eps = 1e-10
diff --git a/src/ForceFreeStates/Inertia.jl b/src/ForceFreeStates/Inertia.jl
new file mode 100644
index 000000000..91a669c13
--- /dev/null
+++ b/src/ForceFreeStates/Inertia.jl
@@ -0,0 +1,358 @@
+"""
+    PerpendicularInertia(; rho)
+
+Perpendicular plasma inertia for the force-balance equations, the ELITE inertia model: the kinetic
+energy K = ½∫ρ|ξ_⊥|² dV with the parallel inertia dropped (incompressible displacement). With the
+normal displacement Ξ = ξ·∇ψ and the in-surface variable ξ_s = (ξ×B)·e_ψ of the energy principle
+[Glasser 2016 Eqs. 8-9], ξ×B = ξ_s∇ψ + Ξu with u = χ′(∇ζ − q∇θ), so |ξ_⊥|² = |ξ×B|²/B² is algebraic
+in the two retained variables and K adds the Hermitian bands ⟨μ₀ρ|∇ψ|²/B²⟩ to A, ⟨μ₀ρ u·∇ψ/B²⟩ to
+C and ⟨μ₀ρ|u|²/B²⟩ to H of [Glasser 2016 Eq. A6], each scaled by γ² for a mode ∝ e^{γt}
+(δW + γ²K is stationary).
+
+Pass it as the `kinetic` source of a [`ForceBalanceProblem`](@ref) with `kinetic_factor = γ²` in
+s⁻², or let [`growth_rate`](@ref) find the γ² at which the free-boundary energy vanishes. The F, K,
+G reduction then takes the kinetic (non-Q-factorized) path, on which the rational surfaces are
+regular points for γ² > 0 and no ideal crossings are needed.
+
+## Fields
+
+  - `rho` - Mass density in kg/m³: a number, or a function of normalized flux ψ_N.
+"""
+@kwdef struct PerpendicularInertia
+    rho::Union{Float64,Function}
+end
+
+_rho_at(src::PerpendicularInertia, psi::Real) = src.rho isa Function ? Float64(src.rho(psi)) : src.rho
+
+"""
+    PerpendicularInertia(equil::Equilibrium.PlasmaEquilibrium)
+
+The inertia of the equilibrium's attached kinetic profiles: ρ(ψ_N) = μ m_p n_e(ψ_N), with μ the mean ion mass per
+electron of the main ions and the impurity at the local Z_eff, the density the SLAYER layer uses. Errors when no
+kinetic profiles are attached (set `kinetic_file` or call `attach_kinetic_profiles!`).
+"""
+function PerpendicularInertia(equil::Equilibrium.PlasmaEquilibrium)
+    data = InnerLayer._layer_kinetic_data(equil, "mass density ρ")
+    _require_profiles(data, "the inertia density ρ = μ m_p n_e")
+    loaded = InnerLayer.SLAYER._slayer_profiles(data, nothing)
+    mu = InnerLayer.SLAYER._layer_mass_per_electron(data, InnerLayer.SLAYER._layer_zeff(loaded, data))
+    n_e = loaded.profiles.n_e
+    return PerpendicularInertia(; rho=psi -> (mu isa Real ? mu : mu(psi)) * Utilities.M_P * n_e(psi))
+end
+
+"""
+    ion_diamagnetic_frequency(equil::Equilibrium.PlasmaEquilibrium, n::Integer; z_i=1) -> ψ_N -> ω_*pi
+
+The ion diamagnetic frequency in rad/s for toroidal mode number `n`, ω_*pi = n (dp_i/dΨ)/(Z_i e n_i) with Ψ the
+poloidal flux per radian, from the attached kinetic profiles with n_i ∝ n_e and p_i = n_i T_i, the layer
+convention (positive for the usual outward-decreasing profiles). The EPED1 criterion counts a mode as unstable
+when γ > ω_*pi/2 at the pedestal.
+"""
+function ion_diamagnetic_frequency(equil::Equilibrium.PlasmaEquilibrium, n::Integer; z_i::Real=1)
+    data = InnerLayer._layer_kinetic_data(equil, "ion diamagnetic frequency ω_*pi")
+    _require_profiles(data, "ω_*pi")
+    p = InnerLayer.SLAYER._slayer_profiles(data, nothing).profiles
+    scale = -n / (z_i * equil.psio)
+    return psi -> scale * (p.T_i(psi) * p.n_e(psi; deriv=DerivOp(1)) / p.n_e(psi) + p.T_i(psi; deriv=DerivOp(1)))
+end
+
+# The shared kinetic-profile loader needs n_e, T_e and T_i even where only some of them enter.
+_has_profiles(data) = data !== nothing && data.n_e !== nothing && data.T_e !== nothing && data.T_i !== nothing
+function _require_profiles(data, what::AbstractString)
+    _has_profiles(data) || throw(ArgumentError("$what needs kinetic data with n_e, T_e and T_i (the shared profile loader reads all three)"))
+    return data
+end
+
+# Toroidal mode number of the least-stable free-boundary mode of `r` (its lowest n without free-boundary energies).
+_least_stable_n(r::ForceFreeStatesResult) = r.free_boundary === nothing ? r.nlow : r.nlow + r.free_boundary.n_tor_idx[1]
+
+# Peak |ω_*pi| over the integration domain of `r` at its least-stable mode's n; NaN without the kinetic profiles.
+function _peak_omega_star_pi(r::ForceFreeStatesResult)
+    _has_profiles(r.equil.kinetic) || return NaN
+    omega = ion_diamagnetic_frequency(r.equil, _least_stable_n(r))
+    psis = filter(x -> r.psilow <= x <= r.psilim, r.metric.xs)
+    return isempty(psis) ? NaN : maximum(abs ∘ omega, psis)
+end
+
+"""
+    inertia_kinetic_matrices(src::PerpendicularInertia, equil, intr, metric) -> (kw_flat, kt_flat)
+
+The perpendicular-inertia bands on the metric's (ψ, θ) grid in the (ψ, mode pair, component) layout
+of [`KineticMatrices`](@ref): components A, C and H carry ⟨μ₀ρ|∇ψ|²/B²⟩, ⟨μ₀ρ u·∇ψ/B²⟩ and
+⟨μ₀ρ|u|²/B²⟩ per unit γ², the others and the torque matrices are zero. The metric stores
+J·g_ij with g_ij = e_i·e_j/J² [Glasser 2016 Eq. A1], so the contravariant products are cofactors:
+|∇ψ|² = g22g33 − g23², ∇ψ·∇θ = g23g31 − g12g33, ∇ψ·∇ζ = g12g23 − g22g31, |∇θ|² = g11g33 − g31²,
+|∇ζ|² = g11g22 − g12², ∇θ·∇ζ = g31g12 − g11g23 (all in the stored J·g units), and
+B² = χ′²(g22 + 2q g23 + q²g33)/J from B = (χ′/J)(e_θ + q e_ζ) [Glasser 2016 Eq. 5].
+"""
+function inertia_kinetic_matrices(src::PerpendicularInertia, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal, metric::MetricData)
+    mpsi, mtheta, mpert = metric.mpsi, metric.mtheta, intr.mpert
+    chi1 = 2π * equil.psio
+    integrands = zeros(Float64, mpsi, mtheta, 3)
+    for ipsi in 1:mpsi
+        q = equil.profiles.q_spline(metric.xs[ipsi])
+        rho = _rho_at(src, metric.xs[ipsi])
+        isfinite(rho) && rho > 0 || error("PerpendicularInertia: ρ(ψ_N = $(metric.xs[ipsi])) = $rho; the mass density must be finite and positive")
+        mu0rho = Utilities.MU_0 * rho
+        for jtheta in 1:mtheta
+            g11, g22, g33, g23, g31, g12, jac = @view metric.fs[ipsi, jtheta, 1:7]
+            gradpsi2 = g22 * g33 - g23^2
+            gradpsi_theta = g23 * g31 - g12 * g33
+            gradpsi_zeta = g12 * g23 - g22 * g31
+            gradtheta2 = g11 * g33 - g31^2
+            gradzeta2 = g11 * g22 - g12^2
+            gradtheta_zeta = g31 * g12 - g11 * g23
+            b2 = chi1^2 * (g22 + 2q * g23 + q^2 * g33) / jac
+            weight = jac * mu0rho / b2
+            integrands[ipsi, jtheta, 1] = weight * gradpsi2
+            integrands[ipsi, jtheta, 2] = weight * chi1 * (gradpsi_zeta - q * gradpsi_theta)
+            integrands[ipsi, jtheta, 3] = weight * chi1^2 * (gradzeta2 - 2q * gradtheta_zeta + q^2 * gradtheta2)
+        end
+    end
+    fc = Utilities.FourierCoefficients(metric.xs, metric.ys, integrands, mpert - 1)
+
+    np = intr.numpert_total
+    kw_flat = zeros(ComplexF64, mpsi, np^2, 6)
+    kt_flat = zeros(ComplexF64, mpsi, np^2, 6)
+    band = zeros(ComplexF64, 2mpert - 1)
+    for ipsi in 1:mpsi, (qty, ic) in ((1, 1), (2, 3), (3, 6))
+        # Same band layout as the ideal matrices: coefficient m at mpert−m, conjugates fill the other half
+        for m in 0:(mpert-1)
+            band[mpert-m] = Utilities.get_complex_coeff(fc, ipsi, m, qty)
+        end
+        for k in 1:(mpert-1)
+            band[k+mpert] = conj(band[mpert-k])
+        end
+        for ipert_n in 1:intr.npert, ipert_m in 1:mpert, jpert_m in 1:mpert
+            ipert = ipert_m + (ipert_n - 1) * mpert
+            jpert = jpert_m + (ipert_n - 1) * mpert
+            kw_flat[ipsi, ipert+(jpert-1)*np, ic] = band[jpert_m-ipert_m+mpert]
+        end
+    end
+    return kw_flat, kt_flat
+end
+
+"""
+    GrowthRateConfig(; rho=NaN, gamma_sq_init=1e9, gamma_sq_min=1e4, gamma_sq_max=1e16, rtol=1e-3, ftol=1e-6, maxiter=30)
+
+The `[GrowthRate]` deck section: its presence runs [`growth_rate`](@ref) after the force-balance stage, and the
+same settings drive the three-argument `growth_rate` method.
+
+## Fields
+
+  - `rho` - Uniform mass density in kg/m³; NaN (the default) derives ρ(ψ_N) from the equilibrium's kinetic profiles.
+  - `gamma_sq_init` - First γ² in s⁻² tried when bracketing the root.
+  - `gamma_sq_min` - γ² floor in s⁻² below which the mode counts as marginal.
+  - `gamma_sq_max` - γ² ceiling in s⁻² beyond which bracketing stops with an error.
+  - `rtol` - Relative width of the γ² bracket at which the root iteration stops.
+  - `ftol` - Fraction of the marginal |λ_min| below which |λ_min| counts as zero.
+  - `maxiter` - Maximum root-iteration steps after bracketing.
+"""
+@kwdef struct GrowthRateConfig
+    rho::Float64 = NaN
+    gamma_sq_init::Float64 = 1e9
+    gamma_sq_min::Float64 = 1e4
+    gamma_sq_max::Float64 = 1e16
+    rtol::Float64 = 1e-3
+    ftol::Float64 = 1e-6
+    maxiter::Int = 30
+end
+
+# The inertia a deck section asks for: uniform, or from the equilibrium's kinetic profiles.
+inertia_source(config::GrowthRateConfig, equil::Equilibrium.PlasmaEquilibrium) =
+    isnan(config.rho) ? PerpendicularInertia(equil) : PerpendicularInertia(; rho=config.rho)
+
+"""
+    IdealGrowthRateResult
+
+The outcome of [`growth_rate`](@ref).
+
+## Fields
+
+  - `gamma::Float64` - Growth rate in s⁻¹ of the least-stable free-boundary mode; 0 when the plasma is marginally stable or stable.
+  - `gamma_sq::Float64` - γ² at which the least-stable free-boundary energy eigenvalue crosses zero.
+  - `converged::Bool` - Whether the root iteration met its `ftol` or `rtol` stop (always true for a stable or marginal plasma); false when
+    `maxiter` ran out, in which case `gamma` is the last iterate and `history` brackets the root, or when any solve flagged an ambiguous
+    step in its conjugate-point count.
+  - `inertia::PerpendicularInertia` - The inertia source the solves used.
+  - `config::Union{Nothing,GrowthRateConfig}` - The deck settings when the solve came from a `[GrowthRate]` section or the three-argument method; `nothing` otherwise.
+  - `omega_star_pi::Float64` - Peak ion diamagnetic frequency |ω_*pi| in rad/s over the integration domain, at the toroidal n of the
+    least-stable free-boundary mode of the solve at `gamma_sq` (of the marginal solve when stable), from the equilibrium's kinetic
+    profiles (NaN without n_e, T_e and T_i); the EPED1 criterion counts the mode as unstable when `gamma > omega_star_pi/2`. The
+    peak is taken over `[psilow, psilim]`, so a truncated domain can cut off part of a pedestal gradient.
+  - `marginal::ForceFreeStatesResult` - The ideal (γ² = 0) solve, whose least-stable energy decides stability.
+  - `result::Union{Nothing,ForceFreeStatesResult}` - The solve at `gamma_sq`, carrying the eigenmode; `nothing` when stable.
+  - `history::Vector` - Every (γ², λ_min, nzero) evaluated, in order, with λ_min the least-stable free-boundary energy at that γ² and
+    `nzero` the exact conjugate-point count of that solve.
+"""
+struct IdealGrowthRateResult
+    gamma::Float64
+    gamma_sq::Float64
+    converged::Bool
+    inertia::PerpendicularInertia
+    config::Union{Nothing,GrowthRateConfig}
+    omega_star_pi::Float64
+    marginal::ForceFreeStatesResult
+    result::Union{Nothing,ForceFreeStatesResult}
+    history::Vector{@NamedTuple{gamma_sq::Float64, lambda::Float64, nzero::Int}}
+end
+
+"""
+    growth_rate(prob::ForceBalanceProblem, alg::Forward; inertia, gamma_sq_init=1e9, gamma_sq_min=1e4, gamma_sq_max=1e16,
+                rtol=1e-3, ftol=1e-6, maxiter=30, marginal=nothing, config=nothing, verbose=false) -> IdealGrowthRateResult
+    growth_rate(prob::ForceBalanceProblem, alg::Forward, config::GrowthRateConfig; marginal=nothing, verbose=false)
+
+Growth rate of the fastest ideal mode of `prob` with the perpendicular inertia `inertia` (a
+[`PerpendicularInertia`](@ref)), ELITE's eigenvalue problem solved by shooting. The inertia enters as
+`kinetic_factor = γ²`; at each γ² the number of modes growing faster than γ is the number of conjugate
+points `nzero` plus one when the least-stable free-boundary energy λ_min is negative, and that count cannot
+rise with γ² since the kinetic energy is positive. The ideal solve gives the count at γ = 0; when it is
+zero the plasma is stable and `gamma` is 0. Otherwise the count is bracketed by decades from
+`gamma_sq_init` (down while nothing is faster, up while something still is). An interior crossing still
+present at the lower end is a fixed-boundary eigenvalue inside the bracket, where λ_min has a pole; since
+the fastest mode is free-boundary and at least as fast as any fixed-boundary one, the lower end is lifted
+above it by bisecting on the count. Then only λ_min changes sign in the bracket and a bracketed secant
+(Illinois) iteration brings it to zero, stopping when |λ_min| falls below `ftol` times the reference
+|λ_min| (the marginal value, or the bracket's negative end when the marginal energy is positive) or the
+bracket is narrower than `rtol` in γ², within `maxiter` steps (a warning and `converged = false`
+otherwise). Every evaluation is one force-balance solve. `gamma_sq_min` is the floor below which the mode
+counts as marginal (`gamma` = 0), `gamma_sq_max` the ceiling beyond which bracketing stops with an error.
+When the count bisection reaches `rtol` with a crossing still present, the fastest mode is a fixed-boundary
+one (its edge displacement vanishes) and the bracket's upper end is reported as converged.
+
+The conjugate points are counted exactly, as the eigenphases of the Cayley transform of the solution frame
+that wrap through π along ψ (its Maslov index), not by the smallest-eigenvalue detector `diagnostics.nzero`;
+a step whose phase match is ambiguous is warned about and marks the result not converged. Starting
+`gamma_sq_init` above the fastest growth rate expected is cheapest, since the upward scan needs only the sign of λ_min.
+
+A `marginal` result of the same problem (an earlier ideal `solve`) is reused instead of solving it again;
+`config` is recorded on the result and is what the three-argument method passes. The problem must carry no
+kinetic source. The inertial solves take every other control of the problem, so a two-pass ψ grid is the
+same as the marginal one; with `singfac_min > 0` the kinetic surface finder can still place crossings at
+the sharp cond(F̂) peaks that appear near rational surfaces when γ²ρ is far below any physical root, where
+only the sign of λ_min matters, and `singfac_min = 0` gives a single chunk (the rational surfaces are
+regular with inertia).
+"""
+function growth_rate(prob::ForceBalanceProblem, alg::Forward; inertia::PerpendicularInertia,
+    gamma_sq_init::Real=1e9, gamma_sq_min::Real=1e4, gamma_sq_max::Real=1e16, rtol::Real=1e-3, ftol::Real=1e-6, maxiter::Int=30,
+    marginal::Union{Nothing,ForceFreeStatesResult}=nothing, config::Union{Nothing,GrowthRateConfig}=nothing, verbose::Bool=false)
+    get(prob.ctrl_kwargs, :kinetic, nothing) === nothing ||
+        error("growth_rate needs an ideal marginal problem: remove its kinetic source, the inertia is the only kinetic term")
+    history = @NamedTuple{gamma_sq::Float64, lambda::Float64, nzero::Int}[]
+    lambda_min(r) = minimum(real, r.free_boundary.et)
+    counted_solve(p) = _solve_force_balance(p, alg; verbose, count_conjugate_points=true)
+    counted(r) = r.diagnostics !== nothing && r.diagnostics.nzero_conj >= 0
+    flagged = Ref(false)
+    function check_flags(r, gamma_sq)
+        nflag = r.diagnostics.nconj_flags
+        nflag > 0 || return
+        flagged[] = true
+        @warn "growth_rate: the conjugate-point count at γ² = $gamma_sq had $nflag ambiguous step(s); the result is marked not converged"
+    end
+    faster(h) = h.nzero + (h.lambda < 0)   # modes growing faster than γ; non-increasing in γ²
+    solves = Dict{Float64,ForceFreeStatesResult}()
+    function finish(gamma_sq, converged)
+        r = get(solves, gamma_sq, nothing)
+        return IdealGrowthRateResult(sqrt(gamma_sq), gamma_sq, converged && !flagged[], inertia, config, _peak_omega_star_pi(something(r, marginal)), marginal, r, history)
+    end
+
+    if marginal !== nothing && !counted(marginal)
+        verbose && @info "growth_rate: the supplied marginal solve has no conjugate-point count; solving it again with the count"
+        marginal = nothing
+    end
+    marginal = marginal === nothing ? counted_solve(prob) : marginal
+    marginal.free_boundary === nothing && error("growth_rate needs the free-boundary energies: set vac_flag=true")
+    check_flags(marginal, 0.0)
+    push!(history, (gamma_sq=0.0, lambda=lambda_min(marginal), nzero=marginal.diagnostics.nzero_conj))
+    verbose && @info @sprintf("growth_rate: marginal λ_min = %.6e (nzero = %d)", history[end].lambda, history[end].nzero)
+    faster(history[1]) == 0 && return finish(0.0, true)
+
+    kwargs = copy(prob.ctrl_kwargs)
+    kwargs[:kinetic] = inertia
+    kwargs[:kinetic_factor] = Float64(gamma_sq_init)
+    kwargs[:kinetic] = kinetic_matrices(ForceBalanceProblem(prob.equil, prob.wall, kwargs), alg; verbose)   # the bands once, reused at every γ²
+    seen = Dict{Float64,Int}()   # γ² → history index; each γ² is one solve
+    function evaluate(gamma_sq)
+        gamma_sq = Float64(gamma_sq)
+        haskey(seen, gamma_sq) && return history[seen[gamma_sq]]
+        kw = copy(kwargs)
+        kw[:kinetic_factor] = gamma_sq
+        r = counted_solve(ForceBalanceProblem(prob.equil, prob.wall, kw))
+        solves[gamma_sq] = r
+        check_flags(r, gamma_sq)
+        push!(history, (gamma_sq=gamma_sq, lambda=lambda_min(r), nzero=r.diagnostics.nzero_conj))
+        seen[gamma_sq] = length(history)
+        verbose && @info @sprintf("growth_rate: γ² = %.4e  λ_min = %.6e  nzero = %d", gamma_sq, history[end].lambda, history[end].nzero)
+        return history[end]
+    end
+
+    # Bracket the count by decades: down while nothing is faster, up while something still is
+    hi = Float64(gamma_sq_init)
+    h_hi = evaluate(hi)
+    while faster(h_hi) == 0
+        hi /= 10
+        if hi < gamma_sq_min
+            verbose && @info @sprintf("growth_rate: nothing grows faster than γ² = %.2e, below gamma_sq_min: marginal", 10hi)
+            return finish(0.0, true)
+        end
+        h_hi = evaluate(hi)
+    end
+    lo, h_lo = hi, h_hi
+    while faster(h_hi) > 0
+        lo, h_lo = hi, h_hi
+        hi *= 10
+        hi > gamma_sq_max && error("growth_rate: a mode still grows at γ² = $lo; raise gamma_sq_max or check the inertia density")
+        h_hi = evaluate(hi)
+    end
+    # The fastest root lies in (lo, hi]. A crossing still present at lo is a fixed-boundary eigenvalue inside the bracket, where
+    # λ_min has a pole; the fastest mode is free-boundary and at least as fast (interlacing), so bisect on the count, which is
+    # monotone, until lo is above that eigenvalue and only the free-boundary energy changes sign in the bracket
+    while h_lo.nzero > 0 && hi - lo > rtol * hi
+        mid = sqrt(lo * hi)
+        h_mid = evaluate(mid)
+        if faster(h_mid) == 0
+            hi, h_hi = mid, h_mid
+        else
+            lo, h_lo = mid, h_mid
+        end
+    end
+    if h_lo.nzero > 0   # the root sits on the pole to within rtol: the fastest mode has no edge displacement
+        verbose && @info @sprintf("growth_rate: the fastest mode is fixed-boundary to within rtol, γ² in [%.6e, %.6e]", lo, hi)
+        return finish(hi, true)
+    end
+
+    # Illinois iteration on the sign change; the retained end's λ is halved when it is kept twice running
+    scale = history[1].lambda < 0 ? abs(history[1].lambda) : abs(h_lo.lambda)
+    flo, fhi = h_lo.lambda, h_hi.lambda
+    gamma_sq, side, converged = hi, 0, false
+    for _ in 1:maxiter
+        gamma_sq = (lo * fhi - hi * flo) / (fhi - flo)
+        h = evaluate(gamma_sq)
+        h.nzero == 0 || error("growth_rate: an interior zero crossing appeared inside the bracket at γ² = $gamma_sq; the free-boundary energy is not monotone there")
+        lam = h.lambda
+        if abs(lam) <= ftol * scale
+            converged = true
+            break
+        end
+        if lam < 0
+            lo, flo = gamma_sq, lam
+            side == -1 && (fhi /= 2)
+            side = -1
+        else
+            hi, fhi = gamma_sq, lam
+            side == 1 && (flo /= 2)
+            side = 1
+        end
+        if hi - lo <= rtol * hi
+            converged = true
+            break
+        end
+    end
+    converged || @warn "growth_rate: no convergence in $maxiter steps; γ² lies in [$lo, $hi] and the last iterate $gamma_sq is reported"
+    return finish(gamma_sq, converged)
+end
+
+function growth_rate(prob::ForceBalanceProblem, alg::Forward, config::GrowthRateConfig; marginal::Union{Nothing,ForceFreeStatesResult}=nothing, verbose::Bool=false)
+    return growth_rate(prob, alg; inertia=inertia_source(config, prob.equil), config.gamma_sq_init, config.gamma_sq_min, config.gamma_sq_max,
+        config.rtol, config.ftol, config.maxiter, marginal, config, verbose)
+end
diff --git a/src/ForceFreeStates/InertiaOutput.jl b/src/ForceFreeStates/InertiaOutput.jl
new file mode 100644
index 000000000..9612e60c8
--- /dev/null
+++ b/src/ForceFreeStates/InertiaOutput.jl
@@ -0,0 +1,44 @@
+# Writer for an `IdealGrowthRateResult`. Layout under the caller's parent group (inputs live only under Input/):
+#
+#   ForceFreeStates/GrowthRate/
+#   ├── gamma, gamma_sq, converged, omega_star_pi -- growth rate, its square, the stop flag and the EPED1 ω_*pi reference
+#   ├── psi, rho                       -- inertia density on the marginal solve's ψ_N grid
+#   └── History/                       -- every (γ², λ_min, nzero) evaluated by the root search
+
+const GROWTH_RATE_H5_ANNOTATIONS = [
+    "gamma" => (; long_name="growth rate γ of the least-stable free-boundary mode with perpendicular inertia (0 when stable)", units="s^-1"),
+    "gamma_sq" => (; long_name="γ² at which the least-stable free-boundary energy eigenvalue crosses zero (0 when stable)", units="s^-2"),
+    "converged" => (; long_name="1 when the root iteration met its stop criterion, 0 when maxiter ran out (gamma is then the last iterate)"),
+    "omega_star_pi" => (; long_name="peak ion diamagnetic frequency |ω_*pi| over the integration domain at the lowest n (NaN without kinetic T_i)",
+        units="rad/s"),
+    "psi" => (; long_name="normalized poloidal flux ψ_N grid of the marginal solve", dims=("psi",), scale="psi"),
+    "rho" => (; long_name="perpendicular-inertia mass density ρ", units="kg/m^3", dims=("psi",), attach=(1 => "psi",)),
+    "History/gamma_sq" => (; long_name="γ² evaluated by the growth-rate root search, in order", units="s^-2", dims=("iteration",)),
+    "History/lambda_min" => (; long_name="power-normalized least-stable free-boundary energy at each evaluated γ² (per unit ⟨|ξ|²⟩)", dims=("iteration",)),
+    "History/nzero" => (; long_name="number of conjugate points (fixed-boundary modes inside the ψ interval) at each evaluated γ²", dims=("iteration",))
+]
+
+"""
+    save_to_hdf5(parent::Union{HDF5.File,HDF5.Group}, result::IdealGrowthRateResult)
+
+Write `result` into a `ForceFreeStates/GrowthRate/` subgroup of `parent`, replacing any existing one, and record its
+inputs in `Input/`. The growth rate, its square and the peak ω_*pi are scalars; the inertia density is sampled on
+the marginal solve's ψ_N grid; `History/` holds every (γ², λ_min, nzero) the root search evaluated.
+"""
+function save_to_hdf5(parent::Union{HDF5.File,HDF5.Group}, result::IdealGrowthRateResult)
+    ffs = haskey(parent, "ForceFreeStates") ? parent["ForceFreeStates"] : create_group(parent, "ForceFreeStates")
+    g = Utilities.HDF5Annotations.replace_group!(ffs, "GrowthRate")
+    g["gamma"] = result.gamma
+    g["gamma_sq"] = result.gamma_sq
+    g["converged"] = Int(result.converged)
+    g["omega_star_pi"] = result.omega_star_pi
+    psi = collect(Float64, result.marginal.metric.xs)
+    g["psi"] = psi
+    g["rho"] = [_rho_at(result.inertia, x) for x in psi]
+    g["History/gamma_sq"] = [h.gamma_sq for h in result.history]
+    g["History/lambda_min"] = [h.lambda for h in result.history]
+    g["History/nzero"] = [h.nzero for h in result.history]
+    Utilities.HDF5Annotations.annotate!(g, GROWTH_RATE_H5_ANNOTATIONS)
+    result.config === nothing || Utilities.record_inputs!(parent, result)   # the API form has no deck section to record
+    return g
+end
diff --git a/src/ForceFreeStates/Kinetic.jl b/src/ForceFreeStates/Kinetic.jl
index 51351e556..b1b831a24 100644
--- a/src/ForceFreeStates/Kinetic.jl
+++ b/src/ForceFreeStates/Kinetic.jl
@@ -142,6 +142,9 @@ function evaluate_kinetic_matrices(
     if ctrl.kinetic isa TestKineticMatrices
         kw_flat, kt_flat = fixed_kinetic_matrices(intr.mpert, intr.numpert_total, mpsi, ctrl.kinetic.scale, intr.mlow, mats, xs)
         return matrices(collect(xs), kw_flat, kt_flat)
+    elseif ctrl.kinetic isa PerpendicularInertia
+        kw_flat, kt_flat = inertia_kinetic_matrices(ctrl.kinetic, equil, intr, metric)
+        return matrices(collect(xs), kw_flat, kt_flat)
     end
     isnothing(calculated_source) && error(
         "The KineticForces kinetic source requires the KineticForces callback: solve the problem, or pass " *
diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl
index bda2b7c8d..504423aac 100755
--- a/src/GeneralizedPerturbedEquilibrium.jl
+++ b/src/GeneralizedPerturbedEquilibrium.jl
@@ -71,7 +71,9 @@ include("Rerun.jl")
 
 # Scripting-API surface: the integrator selectors, the published result, the equilibrium
 # constructor and the forcing description, re-exported so a user needs one `using`.
-using .ForceFreeStates: AbstractIntegrator, Forward, Riccati, Galerkin, TestKineticMatrices, KineticMatrices, ForceFreeStatesResult
+using .ForceFreeStates:
+    AbstractIntegrator, Forward, Riccati, Galerkin, TestKineticMatrices, KineticMatrices, ForceFreeStatesResult, PerpendicularInertia, IdealGrowthRateResult, growth_rate,
+    GrowthRateConfig, ion_diamagnetic_frequency
 using .ForceFreeStates: MatchProblem
 using .InnerLayer: GGJModel, SLAYERModel
 using .Tearing: TearingProblem, TearingConfig, TearingResult, CriticalResonantFieldControl
@@ -140,10 +142,10 @@ equilibrium and writes nothing. The raw forcing file contents are captured in
 `Input/RawInputs/ForcingTerms/`, and `preloaded_coil_sets` lets the rerun path inject the coil
 geometry stored in `Input/RawInputs/Coils/`.
 
-Returns `(; equil, local_stability, ffs, match, resp, pe, kinetic, error_fields, tearing)`, each
+Returns `(; equil, local_stability, ffs, match, growth_rate, resp, pe, kinetic, error_fields, tearing)`, each
 `nothing` when its section is absent. `equil` is the equilibrium the stages used (the solve's, on
 the two-pass grid), `ffs` the result the later stages read (the matched one when the deck has
-`[Match]`), and `match` its `MatchResult`.
+`[Match]`), `match` its `MatchResult`, and `growth_rate` the `IdealGrowthRateResult` of `[GrowthRate]`.
 """
 function main_from_inputs(
     inputs::Dict{String,Any},
@@ -158,7 +160,7 @@ function main_from_inputs(
     runtimes = Vector{Pair{String,Float64}}()
     deck = parse_deck(inputs, path)
     (; output_file, verbose) = deck.run
-    local_stability = ffs = match = resp = pe = kinetic = error_fields = tearing = nothing
+    local_stability = ffs = match = growth = resp = pe = kinetic = error_fields = tearing = nothing
 
     equil = _timed_stage(runtimes, "equilibrium", "Equilibrium") do
         Equilibrium.PlasmaEquilibrium(eq_config, additional_input)
@@ -169,7 +171,9 @@ function main_from_inputs(
 
     if deck.force_balance !== nothing
         ffs = _timed_stage(runtimes, "force_free_states", "Force-Free States") do
-            solve(ForceBalanceProblem(equil; deck.force_balance.keywords...), deck.force_balance.algorithm; verbose)
+            # growth_rate reuses this solve as its marginal, so it carries the conjugate-point count when a [GrowthRate] section follows
+            ForceFreeStates._solve_force_balance(ForceBalanceProblem(equil; deck.force_balance.keywords...), deck.force_balance.algorithm; verbose,
+                count_conjugate_points=deck.growth_rate !== nothing)
         end
         # Every later stage, local stability included, sees the equilibrium the solve formed (re-formed on the two-pass grid).
         equil = ffs.equil
@@ -188,6 +192,12 @@ function main_from_inputs(
             save_to_hdf5(target, ffs; git_version, local_stability)
         end
         @info "Results written to $output_file"
+        if deck.growth_rate !== nothing
+            growth = _timed_stage(runtimes, "growth_rate", "GrowthRate") do
+                growth_rate(ffs.problem, deck.force_balance.algorithm, deck.growth_rate; marginal=ffs, verbose)   # the problem the marginal solve came from
+            end
+            save_to_hdf5(output_file, growth)
+        end
     end
 
     if deck.perturbed_equilibrium !== nothing
@@ -229,7 +239,7 @@ function main_from_inputs(
     ffs === nothing || _write_runtimes!(output_file, runtimes)
     @info "\n$_BANNER\n  GPEC completed successfully in $(@sprintf("%.3f", total_dt)) s\n$_BANNER"
 
-    return (; equil, local_stability, ffs, match, resp, pe, kinetic, error_fields, tearing)
+    return (; equil, local_stability, ffs, match, growth_rate=growth, resp, pe, kinetic, error_fields, tearing)
 end
 
 """
@@ -305,6 +315,7 @@ export main, write_imas
 export solve, plasma_response, perturbed_equilibrium, kinetic_torque, error_fields
 export KineticProfileData, IonSpecies, attach_kinetic_profiles!
 export PlasmaEquilibrium, ForceBalanceProblem, Forward, Riccati, Galerkin, TestKineticMatrices, KineticMatrices, kinetic_matrices, ForceFreeStatesResult, RMPField, PlasmaResponse
+export PerpendicularInertia, IdealGrowthRateResult, GrowthRateConfig, growth_rate, ion_diamagnetic_frequency
 export MatchProblem, TearingProblem, TearingConfig, TearingResult, GGJModel, SLAYERModel
 export LocalStabilityProblem, Ballooning, LocalStabilityResult
 export KineticForcesControl, KineticForcesResult, ErrorFieldsResult
diff --git a/src/Utilities/HDF5Schema.jl b/src/Utilities/HDF5Schema.jl
index f17b53511..f66e0fca2 100644
--- a/src/Utilities/HDF5Schema.jl
+++ b/src/Utilities/HDF5Schema.jl
@@ -289,6 +289,7 @@ const RUNTIME_H5_ANNOTATIONS = [
     "Info/Runtimes/equilibrium" => (; long_name="wall-clock time of the Equilibrium construction stage", units="s"),
     "Info/Runtimes/local_stability" => (; long_name="wall-clock time of the LocalStability stage (Mercier, resistive interchange and ballooning scans)", units="s"),
     "Info/Runtimes/force_free_states" => (; long_name="wall-clock time of the ForceFreeStates stability stage", units="s"),
+    "Info/Runtimes/growth_rate" => (; long_name="wall-clock time of the ForceFreeStates growth-rate stage (inertial γ² root search)", units="s"),
     "Info/Runtimes/match" => (; long_name="wall-clock time of the [Match] inner-layer matching stage", units="s"),
     "Info/Runtimes/tearing" => (; long_name="wall-clock time of the Tearing stage (SLAYER inner-layer solve)", units="s"),
     "Info/Runtimes/perturbed_equilibrium" => (; long_name="wall-clock time of the PerturbedEquilibrium stage", units="s"),
diff --git a/test/runtests_inertia.jl b/test/runtests_inertia.jl
new file mode 100644
index 000000000..368dac2f3
--- /dev/null
+++ b/test/runtests_inertia.jl
@@ -0,0 +1,217 @@
+using TOML
+using LinearAlgebra
+using HDF5
+include("h5_metadata_check.jl")
+
+# Perpendicular inertia (ELITE's kinetic energy K = ½∫ρ|ξ_⊥|² dV) as a kinetic source of the force-balance
+# equations: the A, C, H bands it adds, the B² identity they rely on, the γ²-scaled solves and growth_rate.
+@testset "perpendicular inertia" begin
+    GPEC = GeneralizedPerturbedEquilibrium
+    FFS = GPEC.ForceFreeStates
+    template = joinpath(@__DIR__, "test_data", "regression_tj_ideal_example")
+    deck = TOML.parsefile(joinpath(template, "gpec.toml"))
+    equil = GPEC.Equilibrium.PlasmaEquilibrium(GPEC.equilibrium_config(deck, template), GPEC.analytic_config(deck))
+    wall = GPEC.Vacuum.parse_wall(get(deck, "Wall", nothing))
+    ffs_kwargs = Dict{Symbol,Any}(p.first => p.second for p in GPEC.force_balance_keywords(deck, template) if !(p.first in (:nn, :wall)))
+    ffs_kwargs[:vac_flag] = true
+    forward = Forward(; eulerlagrange_tolerance=1e-8)
+
+    mlow, mhigh = -3, 4
+    mpert = mhigh - mlow + 1
+    metric = FFS.make_metric(equil, mpert)
+    intr = FFS.ForceFreeStatesInternal(; mpert, numpert_total=mpert, mlow, mhigh, nlow=1, nhigh=1, npert=1)
+    ipsi_mid = cld(metric.mpsi, 2)
+    bands(kw, ic, ipsi) = reshape(kw[ipsi, :, ic], mpert, mpert)
+
+    @testset "A, C, H bands" begin
+        kw, kt = FFS.inertia_kinetic_matrices(FFS.PerpendicularInertia(; rho=1e-7), equil, intr, metric)
+        @test size(kw) == (metric.mpsi, mpert^2, 6)
+        @test all(iszero, kt)
+        @test all(iszero, kw[:, :, [2, 4, 5]])
+        for ic in (1, 3, 6)
+            m = bands(kw, ic, ipsi_mid)
+            @test norm(m - m') <= 1e-12 * norm(m)
+            @test norm(m) > 0
+        end
+        @test all(>(0), eigvals(Hermitian(bands(kw, 1, ipsi_mid))))
+        @test all(>(0), eigvals(Hermitian(bands(kw, 6, ipsi_mid))))
+
+        kw2, _ = FFS.inertia_kinetic_matrices(FFS.PerpendicularInertia(; rho=2e-7), equil, intr, metric)
+        @test kw2 ≈ 2kw rtol = 1e-12
+
+        # A ψ-dependent density scales each surface's bands by ρ(ψ)/ρ₀ and so lowers them toward the edge
+        kwf, _ = FFS.inertia_kinetic_matrices(FFS.PerpendicularInertia(; rho=psi -> 1e-7 * (1 - psi)), equil, intr, metric)
+        for ipsi in (2, ipsi_mid, metric.mpsi - 1), ic in (1, 3, 6)
+            @test bands(kwf, ic, ipsi) ≈ (1 - metric.xs[ipsi]) * bands(kw, ic, ipsi) rtol = 1e-12
+        end
+        diagA = [real(kwf[ipsi, 1, 1]) / real(kw[ipsi, 1, 1]) for ipsi in 1:metric.mpsi]
+        @test issorted(diagA; rev=true)
+    end
+
+    @testset "B² from the metric matches the equilibrium field" begin
+        chi1 = 2π * equil.psio
+        bnodes = equil.eqfun_B.nodal_derivs.partials
+        @test size(bnodes)[2:3] == (metric.mpsi, metric.mtheta)
+        for ipsi in (3, ipsi_mid, metric.mpsi - 2), jtheta in (1, cld(metric.mtheta, 3), cld(metric.mtheta, 2))
+            g11, g22, g33, g23, g31, g12, jac = metric.fs[ipsi, jtheta, 1:7]
+            q = equil.profiles.q_spline(metric.xs[ipsi])
+            @test chi1^2 * (g22 + 2q * g23 + q^2 * g33) / jac ≈ bnodes[1, ipsi, jtheta]^2 rtol = 1e-3
+        end
+    end
+
+    lambda_min(r) = minimum(real, r.free_boundary.et)
+    ideal = solve(equil, forward; nn=1, wall, ffs_kwargs...)
+    inertial_kwargs = merge(ffs_kwargs, Dict{Symbol,Any}(:kinetic => FFS.PerpendicularInertia(; rho=1e-7), :singfac_min => 0.0))
+
+    @testset "inertial solves" begin
+        @test equil.kinetic === nothing
+        weak = solve(equil, forward; nn=1, wall, inertial_kwargs..., kinetic_factor=1e4)
+        @test weak.closure === :regular
+        @test weak.kinetic.kmsing == 0
+
+        # The inertial layers at the rational surfaces add a term linear in γ, so λ_min approaches the
+        # ideal value as √γ²: a decade in γ² shrinks the offset by √10, and the offset at γ² = 1e2 is ~1.5e-4
+        # of λ (measured 6.5e-5 absolute on this fixture).
+        weaker = solve(equil, forward; nn=1, wall, inertial_kwargs..., kinetic_factor=1e2)
+        offset_weak = lambda_min(weak) - lambda_min(ideal)
+        offset_weaker = lambda_min(weaker) - lambda_min(ideal)
+        @test 0 < offset_weaker < 3e-4 * lambda_min(ideal)
+        @test 8 < offset_weak / offset_weaker < 12
+
+        strong = solve(equil, forward; nn=1, wall, inertial_kwargs..., kinetic_factor=1e10)
+        @test lambda_min(strong) > lambda_min(weak)
+
+        # The kinetic crossing machinery (singfac_min > 0, det F̄ near-zeros) and the single chunk agree
+        crossed = solve(equil, forward; nn=1, wall, merge(inertial_kwargs, Dict{Symbol,Any}(:singfac_min => 1e-4))..., kinetic_factor=1e4)
+        @test crossed.kinetic.kmsing > 0
+        @test lambda_min(crossed) ≈ lambda_min(weak) rtol = 1e-5
+    end
+
+    @testset "growth_rate" begin
+        @test_throws ErrorException solve(equil, forward; nn=1, wall, ffs_kwargs..., kinetic_factor=2.0)
+
+        # The fixture is ideal-stable, so growth_rate stops after the marginal solve
+        @test lambda_min(ideal) > 0
+        prob = ForceBalanceProblem(equil; nn=1, wall, ffs_kwargs...)
+        gr = FFS.growth_rate(prob, forward; inertia=FFS.PerpendicularInertia(; rho=1e-7))
+        @test gr.gamma == 0
+        @test gr.gamma_sq == 0
+        @test gr.result === nothing
+        @test length(gr.history) == 1
+        @test gr.history[1].lambda ≈ lambda_min(ideal) rtol = 1e-12
+    end
+
+    @testset "rho and omega_star_pi from kinetic profiles" begin
+        ktemplate = joinpath(@__DIR__, "test_data", "regression_tj_kinetic_calculated")
+        kdeck = TOML.parsefile(joinpath(ktemplate, "gpec.toml"))
+        kequil = GPEC.Equilibrium.PlasmaEquilibrium(GPEC.equilibrium_config(kdeck, ktemplate), GPEC.analytic_config(kdeck))
+        @test kequil.kinetic !== nothing
+        src = FFS.PerpendicularInertia(kequil)
+        # The density the SLAYER layer uses at a rational surface
+        @test FFS._rho_at(src, 0.5) ≈ only(GPEC.InnerLayer.layer_rho(kequil, [(psifac=0.5,)])) rtol = 1e-12
+        @test FFS._rho_at(src, 0.3) > 0
+        # ω_*pi = n (dp_i/dΨ)/(Z e n_i) with p_i = n_e T_i, against a finite difference of the profile splines
+        p = GPEC.InnerLayer.SLAYER._slayer_profiles(kequil.kinetic, nothing).profiles
+        p_i(ψ) = p.n_e(ψ) * p.T_i(ψ)
+        h = 1e-4
+        for n in (1, 3)
+            omega = FFS.ion_diamagnetic_frequency(kequil, n)
+            fd = -n / kequil.psio * (p_i(0.5 + h) - p_i(0.5 - h)) / (2h) / p.n_e(0.5)
+            @test omega(0.5) ≈ fd rtol = 1e-6
+        end
+        @test FFS.ion_diamagnetic_frequency(kequil, 2)(0.5) ≈ 2 * FFS.ion_diamagnetic_frequency(kequil, 1)(0.5)
+        @test FFS.ion_diamagnetic_frequency(kequil, 1)(0.5) > 0   # outward-decreasing p_i: the layer sign convention
+        @test_throws ArgumentError FFS.ion_diamagnetic_frequency(equil, 1)   # the ideal fixture carries no kinetic data
+        @test_throws ArgumentError FFS.PerpendicularInertia(equil)
+
+        # The kinetic TJ fixture is ideal-stable too: γ = 0 with the domain's peak ω_*pi reported
+        kkw = Dict{Symbol,Any}(q.first => q.second for q in GPEC.force_balance_keywords(kdeck, ktemplate) if !(q.first in (:nn, :wall, :kinetic, :kinetic_factor)))
+        kkw[:vac_flag] = true
+        kprob = ForceBalanceProblem(kequil; nn=1, wall=GPEC.Vacuum.parse_wall(get(kdeck, "Wall", nothing)), kkw...)
+        gr = FFS.growth_rate(kprob, forward; inertia=src)
+        @test gr.gamma == 0 && isfinite(gr.omega_star_pi) && gr.omega_star_pi > 0
+        omega1 = FFS.ion_diamagnetic_frequency(kequil, 1)
+        @test gr.omega_star_pi ≈ maximum(abs ∘ omega1, filter(x -> gr.marginal.psilow <= x <= gr.marginal.psilim, gr.marginal.metric.xs))
+        # The config form reuses a marginal solve and derives the same inertia
+        cfg = FFS.GrowthRateConfig()
+        gr2 = FFS.growth_rate(kprob, forward, cfg; marginal=gr.marginal)
+        @test gr2.gamma == gr.gamma && gr2.omega_star_pi == gr.omega_star_pi && gr2.marginal === gr.marginal
+        @test gr.config === nothing && gr.inertia === src && gr2.config === cfg && isnan(gr2.config.rho)
+        @test FFS._rho_at(FFS.inertia_source(cfg, kequil), 0.4) == FFS._rho_at(src, 0.4)
+        @test FFS._rho_at(FFS.inertia_source(FFS.GrowthRateConfig(; rho=2.0), kequil), 0.4) == 2.0
+        # The HDF5 writer: the group, its annotations and the recorded [GrowthRate] inputs
+        mktempdir() do dir
+            path = joinpath(dir, "gpec.h5")
+            FFS.save_to_hdf5(path, gr2)
+            h5open(path, "r") do f
+                @test isempty(_collect_metadata_violations(f))
+                @test read(f["ForceFreeStates/GrowthRate/gamma"]) == 0.0
+                @test read(f["ForceFreeStates/GrowthRate/converged"]) == 1
+                @test read(f["ForceFreeStates/GrowthRate/omega_star_pi"]) == gr2.omega_star_pi
+                @test length(read(f["ForceFreeStates/GrowthRate/History/lambda_min"])) == 1
+                @test read(f["ForceFreeStates/GrowthRate/rho"]) ≈ [FFS._rho_at(src, x) for x in read(f["ForceFreeStates/GrowthRate/psi"])]
+                @test haskey(TOML.parse(read(f["Input/gpec_toml_raw"])), "GrowthRate")
+            end
+        end
+    end
+
+    @testset "root search on a free-boundary-unstable fixture" begin
+        # q_a = 2.9 makes the TJ fixture external-kink unstable at n = 1
+        kink_deck = merge(deck, Dict{String,Any}("TJAnalytic" => merge(deck["TJAnalytic"], Dict{String,Any}("qa" => 2.9))))
+        kink_equil = GPEC.Equilibrium.PlasmaEquilibrium(GPEC.equilibrium_config(kink_deck, template), GPEC.analytic_config(kink_deck))
+        kink_prob = ForceBalanceProblem(kink_equil; nn=1, wall, ffs_kwargs...)
+        inertia = FFS.PerpendicularInertia(; rho=1e-7)
+        gr = FFS.growth_rate(kink_prob, forward; inertia)
+        @test gr.gamma > 0 && gr.converged && gr.result !== nothing
+        @test all(h.nzero == 0 for h in gr.history)
+        @test abs(minimum(real, gr.result.free_boundary.et)) <= 1e-6 * abs(gr.history[1].lambda)
+        @test gr.result.closure == :regular && isnan(gr.omega_star_pi)
+        # Running out of iterations is reported, not hidden
+        short = @test_logs (:warn, r"no convergence in 1 steps") match_mode = :any FFS.growth_rate(kink_prob, forward; inertia, maxiter=1, marginal=gr.marginal)
+        @test !short.converged && short.gamma > 0 && short.result !== nothing
+        # A root below gamma_sq_min counts as marginal: γ = 0, no inertial solve reported
+        floor = FFS.growth_rate(kink_prob, forward; inertia, gamma_sq_min=1e12, marginal=gr.marginal)
+        @test floor.gamma == 0 && floor.gamma_sq == 0 && floor.result === nothing && floor.converged
+        @test all(h.lambda > 0 for h in floor.history[2:end])
+        # The marginal problem must be ideal, and the density positive
+        @test_throws "ideal marginal problem" FFS.growth_rate(ForceBalanceProblem(kink_equil; nn=1, wall, ffs_kwargs..., kinetic=GPEC.KineticForcesControl()), forward; inertia)
+        @test_throws "finite and positive" FFS.inertia_kinetic_matrices(FFS.PerpendicularInertia(; rho=0.0), equil, intr, metric)
+    end
+
+    @testset "a fixed-boundary mode inside the bracket" begin
+        # p0 = 0.8 adds a fixed-boundary mode at γ² ≈ 1.2e8 (ρ = 1e-7) below the external kink at ≈ 4.9e9, a pole of λ_min
+        # between them. The crossing counter watches the smallest-magnitude eigenvalue of W⁻¹, which a stiff harmonic holds
+        # near zero on this fixture, so it misses the conjugate point; until the count is exact these pin the wrong answers.
+        hp_deck = merge(deck, Dict{String,Any}("TJAnalytic" => merge(deck["TJAnalytic"], Dict{String,Any}("qa" => 2.9, "p0" => 0.8))))
+        hp_equil = GPEC.Equilibrium.PlasmaEquilibrium(GPEC.equilibrium_config(hp_deck, template), GPEC.analytic_config(hp_deck))
+        inertia = FFS.PerpendicularInertia(; rho=1e-7)
+        # Starting below the pole, only the count can carry the search past it to the fastest mode
+        low_start = FFS.growth_rate(ForceBalanceProblem(hp_equil; nn=1, wall, ffs_kwargs...), forward; inertia, gamma_sq_init=1e8)
+        @test low_start.converged
+        @test low_start.gamma_sq > 4e9
+        # A close wall makes the marginal free-boundary energy positive, so only the count says the plasma is unstable, and the
+        # bracket [1e8, 1e9] then holds the pole: the count bisection must lift its lower end before the sign iteration
+        close_wall = GPEC.Vacuum.ConformalWall(; standoff_fraction=0.05)
+        walled = FFS.growth_rate(ForceBalanceProblem(hp_equil; nn=1, wall=close_wall, ffs_kwargs...), forward; inertia, gamma_sq_init=1e8)
+        @test walled.gamma > 0
+        @test any(h.nzero > 0 for h in walled.history)
+    end
+
+    @testset "kinetic data without T_e" begin
+        ktemplate = joinpath(@__DIR__, "test_data", "regression_tj_kinetic_calculated")
+        kdeck = TOML.parsefile(joinpath(ktemplate, "gpec.toml"))
+        kequil = GPEC.Equilibrium.PlasmaEquilibrium(GPEC.equilibrium_config(kdeck, ktemplate), GPEC.analytic_config(kdeck))
+        data = kequil.kinetic
+        GPEC.attach_kinetic_profiles!(kequil, GPEC.Equilibrium.KineticProfileData(; psi=data.psi, n_e=data.n_e, T_i=data.T_i))
+        @test_throws ArgumentError FFS.PerpendicularInertia(kequil)
+        @test_throws ArgumentError FFS.ion_diamagnetic_frequency(kequil, 1)
+        kkw = Dict{Symbol,Any}(q.first => q.second for q in GPEC.force_balance_keywords(kdeck, ktemplate) if !(q.first in (:nn, :wall, :kinetic, :kinetic_factor)))
+        kkw[:vac_flag] = true
+        gr = FFS.growth_rate(
+            ForceBalanceProblem(kequil; nn=1, wall=GPEC.Vacuum.parse_wall(get(kdeck, "Wall", nothing)), kkw...),
+            forward;
+            inertia=FFS.PerpendicularInertia(; rho=1e-7)
+        )
+        @test gr.gamma == 0 && isnan(gr.omega_star_pi)
+    end
+end
diff --git a/test/runtests_solve_api.jl b/test/runtests_solve_api.jl
index 2c17a146b..085ba3bad 100644
--- a/test/runtests_solve_api.jl
+++ b/test/runtests_solve_api.jl
@@ -561,6 +561,84 @@ include("display_check.jl")
             end
         end
 
+        @testset "deck growth rate under [GrowthRate]" begin
+            config = GrowthRateConfig(; rho=1e-7)
+            growth_section = Dict{String,Any}("GrowthRate" => Dict{String,Any}("rho" => 1e-7))
+            # The fixture is ideal-stable, so both paths stop at the marginal solve with γ = 0.
+            mktempdir() do dir
+                ffs = solve(equil, forward; nn=1, wall=wall, ffs_kwargs...)
+                api = growth_rate(ForceBalanceProblem(equil; nn=1, wall=wall, ffs_kwargs...), forward, config; marginal=ffs)
+                reference = GPEC.main([_stage_deck(dir; extra_sections=growth_section)])
+
+                @test reference.growth_rate isa IdealGrowthRateResult
+                @test reference.growth_rate.gamma == api.gamma == 0
+                @test length(reference.growth_rate.history) == length(api.history)
+                @test [h.gamma_sq for h in reference.growth_rate.history] ≈ [h.gamma_sq for h in api.history] rtol = 1e-10
+                @test [h.lambda for h in reference.growth_rate.history] ≈ [h.lambda for h in api.history] rtol = 1e-10
+                @test isequal(reference.growth_rate.omega_star_pi, api.omega_star_pi)
+                @test reference.growth_rate.config == api.config
+                h5file = joinpath(dir, "gpec.h5")
+                @test h5open(f -> haskey(f, "Info/Runtimes/growth_rate"), h5file)
+                @test h5open(f -> haskey(f, "ForceFreeStates/GrowthRate/gamma"), h5file)
+                # The recorded input carries the [GrowthRate] section.
+                recorded = TOML.parse(h5open(f -> read(f, "Input/gpec_toml_raw"), h5file))
+                @test haskey(recorded, "GrowthRate")
+                @test recorded["GrowthRate"]["rho"] == 1e-7
+            end
+
+            # Lowering q_a to 2.9 makes the fixture free-boundary unstable, so the root find runs and γ > 0 is compared.
+            kink_deck = merge(deck, Dict{String,Any}("TJAnalytic" => merge(deck["TJAnalytic"], Dict{String,Any}("qa" => 2.9))))
+            mktempdir() do dir
+                kink_equil = GPEC.Equilibrium.PlasmaEquilibrium(GPEC.equilibrium_config(kink_deck, template), GPEC.analytic_config(kink_deck))
+                ffs = solve(kink_equil, forward; nn=1, wall=wall, ffs_kwargs...)
+                @test minimum(real, ffs.free_boundary.et) < 0
+                api = growth_rate(ForceBalanceProblem(kink_equil; nn=1, wall=wall, ffs_kwargs...), forward, config; marginal=ffs)
+                reference = GPEC.main([_stage_deck(dir; extra_sections=merge(growth_section, Dict{String,Any}("TJAnalytic" => kink_deck["TJAnalytic"])))])
+
+                @test api.gamma > 0
+                @test reference.growth_rate.gamma ≈ api.gamma rtol = 1e-12
+                @test length(reference.growth_rate.history) == length(api.history) > 1
+                @test [h.gamma_sq for h in reference.growth_rate.history] ≈ [h.gamma_sq for h in api.history] rtol = 1e-10
+                @test [h.lambda for h in reference.growth_rate.history] ≈ [h.lambda for h in api.history] rtol = 1e-10
+                @test [h.nzero for h in reference.growth_rate.history] == [h.nzero for h in api.history]
+            end
+
+            parse_staged(dir; kwargs...) = (deck_dir = _stage_deck(dir; kwargs...); GPEC.parse_deck(TOML.parsefile(joinpath(deck_dir, "gpec.toml")), deck_dir))
+            # On the two-pass auto grid the deck must hand growth_rate the problem its marginal solve came from, so both paths re-form from pass 1.
+            auto_overrides = Dict{String,Any}("grid_type" => "auto", "mpsi" => 0, "psi_accuracy" => 3e-2)
+            mktempdir() do dir
+                auto_deck = merge(kink_deck, Dict{String,Any}("Equilibrium" => merge(kink_deck["Equilibrium"], auto_overrides)))
+                pass1 = GPEC.Equilibrium.PlasmaEquilibrium(GPEC.equilibrium_config(auto_deck, template), GPEC.analytic_config(auto_deck))
+                api = growth_rate(ForceBalanceProblem(pass1; nn=1, wall=wall, ffs_kwargs...), forward, config)
+                reference =
+                    GPEC.main([_stage_deck(dir; eq_overrides=auto_overrides, extra_sections=merge(growth_section, Dict{String,Any}("TJAnalytic" => kink_deck["TJAnalytic"])))])
+                @test api.gamma > 0
+                @test reference.growth_rate.gamma ≈ api.gamma rtol = 1e-10
+                @test [h.lambda for h in reference.growth_rate.history] ≈ [h.lambda for h in api.history] rtol = 1e-10
+            end
+
+            # A kinetic source in [ForceFreeStates] would make the marginal solve kinetic-MHD while the inertial ones are ideal.
+            mktempdir() do dir
+                deck_dir = _stage_deck(dir; eq_overrides=Dict{String,Any}("kinetic_file" => "kinetic.dat"), ffs_overrides=Dict{String,Any}("kinetic" => "KineticForces"),
+                    extra_sections=growth_section)
+                cp(joinpath(@__DIR__, "test_data", "regression_tj_kinetic_calculated", "kinetic.dat"), joinpath(deck_dir, "kinetic.dat"))
+                @test_throws "[GrowthRate] needs an ideal force-balance solve" GPEC.parse_deck(TOML.parsefile(joinpath(deck_dir, "gpec.toml")), deck_dir)
+            end
+
+            # The root find needs a Forward free-boundary solve as its plain ideal marginal one, so the deck stops at parse time.
+            mktempdir() do dir
+                @test_throws "[GrowthRate] needs the free-boundary energy" parse_staged(dir; ffs_overrides=Dict{String,Any}("vac_flag" => false), extra_sections=growth_section)
+            end
+            mktempdir() do dir
+                @test_throws "[GrowthRate] needs integrator" parse_staged(dir; ffs_overrides=Dict{String,Any}("integrator" => "Riccati"), extra_sections=growth_section)
+            end
+            mktempdir() do dir
+                matched = merge(growth_section, Dict{String,Any}("Match" => Dict{String,Any}("eta" => 1e-7, "rho" => 1e-7, "rotation" => 1.0)))
+                @test_throws "[GrowthRate] needs the ideal solve as its marginal one" parse_staged(dir; ffs_overrides=Dict{String,Any}("integrator" => "Riccati"),
+                    extra_sections=matched)
+            end
+        end
+
         @testset "kinetic_format" begin
             # A rates-v1 kinetic file read through the deck's [Equilibrium] keys and through the API keyword.
             mktempdir() do dir
diff --git a/test/suites.jl b/test/suites.jl
index dfe79f9e4..8374236d6 100644
--- a/test/suites.jl
+++ b/test/suites.jl
@@ -50,6 +50,7 @@ const TEST_SUITES = [
     "runtests_slayer_runner.jl" => "errorfields",
     "runtests_kinetic.jl" => "core",
     "runtests_kinetic_validity.jl" => "core",
+    "runtests_inertia.jl" => "core",
     "runtests_multiion.jl" => "core",
     "runtests_fullruns.jl" => "api",
     "runtests_torque_response.jl" => "api",