diff --git a/Cslib.lean b/Cslib.lean
index c5bad333..2f40f693 100644
--- a/Cslib.lean
+++ b/Cslib.lean
@@ -2,6 +2,11 @@ module  -- shake: keep-all
 
 public import Cslib.Algorithms.Lean.MergeSort.MergeSort
 public import Cslib.Algorithms.Lean.TimeM
+public import Cslib.CPS.DiscreteLinearTime.AsymptoticStability
+public import Cslib.CPS.DiscreteLinearTime.Basic
+public import Cslib.CPS.DiscreteLinearTime.Cayley
+public import Cslib.CPS.DiscreteLinearTime.Controllability
+public import Cslib.CPS.DiscreteLinearTime.Reachability
 public import Cslib.Computability.Automata.Acceptors.Acceptor
 public import Cslib.Computability.Automata.Acceptors.OmegaAcceptor
 public import Cslib.Computability.Automata.DA.Basic
@@ -43,6 +48,8 @@ public import Cslib.Foundations.Semantics.FLTS.Basic
 public import Cslib.Foundations.Semantics.FLTS.FLTSToLTS
 public import Cslib.Foundations.Semantics.FLTS.LTSToFLTS
 public import Cslib.Foundations.Semantics.FLTS.Prod
+public import Cslib.Foundations.Semantics.GameSemantics.Basic
+public import Cslib.Foundations.Semantics.GameSemantics.HMLGame
 public import Cslib.Foundations.Semantics.LTS.Basic
 public import Cslib.Foundations.Semantics.LTS.Bisimulation
 public import Cslib.Foundations.Semantics.LTS.Simulation
@@ -75,6 +82,7 @@ public import Cslib.Languages.LambdaCalculus.LocallyNameless.Untyped.FullBeta
 public import Cslib.Languages.LambdaCalculus.LocallyNameless.Untyped.FullBetaConfluence
 public import Cslib.Languages.LambdaCalculus.LocallyNameless.Untyped.Properties
 public import Cslib.Languages.LambdaCalculus.Named.Untyped.Basic
+public import Cslib.Logics.HennessyMilnerLogic.Basic
 public import Cslib.Logics.LinearLogic.CLL.Basic
 public import Cslib.Logics.LinearLogic.CLL.CutElimination
 public import Cslib.Logics.LinearLogic.CLL.EtaExpansion
diff --git a/Cslib/CPS/DiscreteLinearTime/AsymptoticStability.lean b/Cslib/CPS/DiscreteLinearTime/AsymptoticStability.lean
new file mode 100644
index 00000000..64325d85
--- /dev/null
+++ b/Cslib/CPS/DiscreteLinearTime/AsymptoticStability.lean
@@ -0,0 +1,228 @@
+
+/-
+Copyright (c) 2026 Bashar Hamade. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Bashar Hamade
+-/
+
+module
+
+public import Cslib.Init
+
+
+public import Mathlib.Analysis.Normed.Group.Basic
+public import Mathlib.Analysis.Normed.Module.RCLike.Real
+public import Mathlib.Analysis.Normed.Operator.ContinuousLinearMap
+public import Mathlib.Analysis.Complex.Order
+public import Mathlib.Analysis.Normed.Algebra.Spectrum
+public import Mathlib.Order.Basic
+public import Mathlib.Analysis.Normed.Algebra.GelfandFormula
+public import Cslib.CPS.DiscreteLinearTime.Basic
+public import Cslib.CPS.DiscreteLinearTime.Cayley
+
+@[expose] public section
+
+/-!
+# Asymptotic Stability Analysis
+
+This module analyzes the asymptotic stability of discrete-time linear systems.
+It uses the Gelfand spectral radius formula to show that if the
+spectral radius of the system matrix is less than one,
+the state converges to zero under zero input.
+
+## Main Theorems
+* `asymptotic_stability_discrete`: The main stability theorem.
+* `gelfand_le_one_when_spectral_radius_le_one`: Relates spectral radius to Gelfand limit.
+
+## References
+https://en.wikipedia.org/wiki/Lyapunov_stability
+https://www.cds.caltech.edu/~murray/books/AM08/pdf/fbs-public_24Jul2020.pdf
+-/
+
+
+
+
+open scoped ComplexOrder
+
+variable {σ ι : Type*}
+variable [NormedAddCommGroup σ] [NormedSpace ℂ σ]
+variable [NormedAddCommGroup ι] [NormedSpace ℂ ι]
+
+
+
+
+/-- System evolution function. -/
+noncomputable def system_evolution (sys : DiscreteLinearSystemState σ ι) : ℕ → σ
+  | 0 => sys.x₀
+  | k + 1 => sys.a (system_evolution sys k) + sys.B (sys.u k)
+
+/-- Discrete State space representation property. -/
+def state_system_equation (sys : DiscreteLinearSystemState σ ι) : Prop :=
+  ∀ k : ℕ, sys.x (k + 1) = sys.a (sys.x k) + sys.B (sys.u k)
+
+
+
+
+
+
+
+/-- Lemma: State evolution under zero input. -/
+lemma state_evolution_zero_input (sys : DiscreteLinearSystemState σ ι)
+ (h_init : sys.x 0 = sys.x₀)
+ (h_state : state_system_equation sys)
+ (h_zero_input : sys.u = zero_input) :
+ ∀ k, sys.x k = (sys.a ^ k) sys.x₀ := by
+  intro k
+  induction k with
+ | zero =>
+   simp [pow_zero, h_init]
+ | succ k ih =>
+   have h1 : sys.x (k + 1) = sys.a (sys.x k) + sys.B (sys.u k) := h_state k
+   rw [ih] at h1
+   rw [h_zero_input] at h1
+   unfold zero_input at h1
+   simp only [ContinuousLinearMap.map_zero, add_zero] at h1
+   rw [h1]
+   rw [←ContinuousLinearMap.comp_apply]
+   congr 1
+   rw [system_power_multiplication_flopped]
+
+
+
+
+
+
+/-- Bound on state norm. -/
+theorem bound_x_norm
+    (sys : DiscreteLinearSystemState σ ι)
+    (hx : ∀ k, sys.x k = (sys.a ^ k) sys.x₀) (N : ℕ) :
+    ∀ k, N < k → ‖sys.x k‖ ≤ ‖sys.a ^ k‖ * ‖sys.x₀‖ := by
+  intro k hk
+  rw [hx k]
+  exact ContinuousLinearMap.le_opNorm (sys.a ^ k) sys.x₀
+
+/-- Spectral radius strictly less than one. -/
+def spectral_radius_less_than_one (a : σ →L[ℂ] σ) : Prop :=
+  spectralRadius ℂ a < 1
+
+/-- If the spectral radius of `a` is less than one, then the
+Gelfand formula limits to less than one. -/
+theorem gelfand_le_one_when_spectral_radius_le_one
+    [CompleteSpace σ] (a : σ →L[ℂ] σ) :
+    spectral_radius_less_than_one a →
+    Filter.limsup (fun (n : ℕ) => (‖a ^ n‖₊ : ENNReal) ^ (1 / n : ℝ)) Filter.atTop < 1 := by
+  intro h_spectral
+  unfold spectral_radius_less_than_one at h_spectral
+  have h_gelfand := spectrum.limsup_pow_nnnorm_pow_one_div_le_spectralRadius a
+  convert lt_of_le_of_lt h_gelfand h_spectral
+
+
+
+/-- If the Gelfand limit is < 1, then eventually the root-term is bounded by some `r < 1`. -/
+theorem gelfand_eventually_bounded [CompleteSpace σ]
+    (a : σ →L[ℂ] σ) (h : Filter.limsup (fun n :
+    ℕ ↦ (‖a ^ n‖₊ : ENNReal) ^ (1 / n : ℝ)) Filter.atTop < 1) :
+    ∃ (r : ENNReal) (N : ℕ), 0 < r ∧ r < 1
+    ∧ ∀ (k : ℕ), N < k → (‖a ^ k‖₊ : ENNReal) ^ (1 / k : ℝ) < r :=
+by
+  obtain ⟨r, h_limsup_lt_r, h_r_lt_one⟩ := exists_between h
+  have r_pos : 0 < r := lt_of_le_of_lt (zero_le _) h_limsup_lt_r
+  have eventually_lt := Filter.eventually_lt_of_limsup_lt h_limsup_lt_r
+  rw [Filter.eventually_atTop] at eventually_lt
+  obtain ⟨N, hN⟩ := eventually_lt
+  use r, N
+  refine ⟨r_pos, h_r_lt_one, fun k hk => hN k (Nat.le_of_lt hk)⟩
+
+
+/-- If terms are eventually bounded by `r^k` with `r < 1`, the sequence tends to zero. -/
+theorem power_to_zero [CompleteSpace σ] (sys : DiscreteLinearSystemState σ ι)
+  (r : ENNReal) (N : ℕ)
+  (r_pos : 0 < r)
+  (r_lt_one : r < 1)
+  (h_power : ∀ (k : ℕ), N < k → (‖sys.a ^ k‖₊ : ENNReal) < r ^ k) :
+  Filter.Tendsto (fun k => ‖sys.a ^ k‖) Filter.atTop (nhds 0) := by
+  have r_real_zero : Filter.Tendsto (fun n => (r ^ n).toReal) Filter.atTop (nhds 0) := by
+    simp only [ENNReal.toReal_pow, tendsto_pow_atTop_nhds_zero_iff, ENNReal.abs_toReal]
+    cases r with
+    | top => simp
+    | coe x =>
+      simp only [ENNReal.coe_toReal, NNReal.coe_lt_one]
+      exact ENNReal.coe_lt_one_iff.mp r_lt_one
+  rw [Metric.tendsto_atTop]
+  intro ε ε_pos
+  obtain ⟨N₁, hN₁⟩ := Metric.tendsto_atTop.mp r_real_zero ε ε_pos
+  use max N N₁ + 1
+  intro k hk
+  have hkN : N < k := lt_of_le_of_lt (le_max_left N N₁) (Nat.lt_of_succ_le hk)
+  have hkN₁ : N₁ ≤ k := le_trans (le_max_right N N₁) (Nat.le_of_succ_le hk)
+  have h_bound := h_power k hkN
+  have h_r_small := hN₁ k hkN₁
+  simp only [dist_zero_right, norm_norm, gt_iff_lt]
+  simp only [ENNReal.toReal_pow, dist_zero_right, norm_pow, Real.norm_eq_abs,
+    ENNReal.abs_toReal] at h_r_small
+  have h_r_ennreal : (r ^ k).toReal < ε := by
+    rw [ENNReal.toReal_pow]
+    exact h_r_small
+  have h_finite : (r ^ k) ≠ ⊤ := by
+    simp only [ne_eq, ENNReal.pow_eq_top_iff]
+    intro h
+    cases h with
+    | intro h_r_top h_k_ne_zero =>
+      have : r < ⊤ := r_lt_one.trans_le le_top
+      exact ne_of_lt this h_r_top
+  calc ‖sys.a ^ k‖
+  = (‖sys.a ^ k‖₊ : ℝ) := by simp
+  _ = (↑‖sys.a ^ k‖₊ : ENNReal).toReal := by simp
+  _ < (r ^ k).toReal := (ENNReal.toReal_lt_toReal ENNReal.coe_ne_top h_finite).mpr h_bound
+  _ < ε := h_r_ennreal
+
+
+
+
+
+
+/-- Main theorem: Asymptotic Stability. If the spectral radius is < 1,
+the zero-input response converges to zero. -/
+theorem asymptotic_stability_discrete [CompleteSpace σ] (sys : DiscreteLinearSystemState σ ι)
+  (h_init : sys.x 0 = sys.x₀)
+  (h_state : state_system_equation sys)
+  (h_zero_input : sys.u = zero_input)
+  (h_spectral : spectral_radius_less_than_one sys.a) :
+  Filter.Tendsto (fun k => ‖sys.x k‖) Filter.atTop (nhds 0) := by
+  have h_gelfand := spectrum.limsup_pow_nnnorm_pow_one_div_le_spectralRadius sys.a
+  have h_gelfand_le_one : Filter.limsup (fun (n : ℕ) =>
+  (‖sys.a ^ n‖₊ : ENNReal) ^ (1 / n : ℝ)) Filter.atTop < 1 := by
+      unfold spectral_radius_less_than_one at h_spectral
+      refine lt_of_le_of_lt ?_ h_spectral
+      exact h_gelfand
+  have eventually_bounded := gelfand_eventually_bounded sys.a h_gelfand_le_one
+  obtain ⟨r, N, r_pos, r_lt_one, h_bound⟩ := eventually_bounded
+  have h_power : ∀ (k : ℕ), N < k → ↑‖sys.a ^ k‖₊ < r ^ k := by
+      intros k' hk'
+      specialize h_bound k' hk'
+      have h_k'_pos : 0 < k' := Nat.zero_lt_of_lt hk'
+      have h_inv_k' : (k' : ℝ)⁻¹ * k' = 1 := by
+        field_simp
+      have h_pow : (↑‖sys.a ^ k'‖₊ ^ (1 / k' : ℝ)) ^ (k' : ℝ) < r ^ k' := by
+        rw [← ENNReal.rpow_natCast r k']
+        exact ENNReal.rpow_lt_rpow h_bound (Nat.cast_pos.mpr h_k'_pos)
+      rw [← ENNReal.rpow_natCast, ← ENNReal.rpow_mul] at h_pow
+      simp only [one_div, ENNReal.rpow_natCast] at h_pow
+      rw [h_inv_k'] at h_pow
+      simp only [ENNReal.rpow_one] at h_pow
+      exact h_pow
+  have hx : ∀ k, sys.x k = (sys.a ^ k) sys.x₀ :=
+      state_evolution_zero_input sys h_init h_state h_zero_input
+  have pow_zero : Filter.Tendsto (fun k => ‖sys.a ^ k‖) Filter.atTop (nhds 0) :=
+      power_to_zero sys r N r_pos r_lt_one h_power
+  simp only [hx]
+  have h_norm_eq : ∀ k, ‖(sys.a ^ k) sys.x₀‖ ≤ ‖sys.a ^ k‖ * ‖sys.x₀‖ :=
+    fun k => ContinuousLinearMap.le_opNorm (sys.a ^ k) sys.x₀
+  have h_prod_zero : Filter.Tendsto (fun k => ‖sys.a ^ k‖ * ‖sys.x₀‖) Filter.atTop (nhds 0) := by
+    rw [← zero_mul ‖sys.x₀‖]
+    exact Filter.Tendsto.mul_const ‖sys.x₀‖ pow_zero
+  exact tendsto_of_tendsto_of_tendsto_of_le_of_le
+    (tendsto_const_nhds)
+    h_prod_zero
+    (fun k => norm_nonneg _)
+    h_norm_eq
diff --git a/Cslib/CPS/DiscreteLinearTime/Basic.lean b/Cslib/CPS/DiscreteLinearTime/Basic.lean
new file mode 100644
index 00000000..32bd2222
--- /dev/null
+++ b/Cslib/CPS/DiscreteLinearTime/Basic.lean
@@ -0,0 +1,109 @@
+/-
+Copyright (c) 2026 Bashar Hamade. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Bashar Hamade
+-/
+
+module
+
+public import Cslib.Init
+public import Mathlib.Analysis.Normed.Module.Basic
+public import Mathlib.Analysis.Normed.Operator.ContinuousLinearMap
+public import Mathlib.LinearAlgebra.Span.Basic
+public import Mathlib.LinearAlgebra.FiniteDimensional.Defs
+public import Mathlib.Data.Finset.Basic
+public import Mathlib.Data.Complex.Basic
+public import Mathlib.Order.CompleteLattice.Basic
+public import Mathlib.Analysis.Complex.Order
+public import Mathlib.Analysis.Normed.Field.Lemmas
+public import Mathlib.Analysis.Normed.Field.Basic
+public import Mathlib.Analysis.Complex.Exponential
+public import Mathlib.Analysis.Normed.Group.Basic
+public import Mathlib.Analysis.Normed.Module.RCLike.Real
+public import Mathlib.Analysis.Normed.Algebra.Spectrum
+public import Mathlib.Order.Basic
+public import Mathlib.LinearAlgebra.Charpoly.Basic
+public import Mathlib.LinearAlgebra.Matrix.Charpoly.LinearMap
+
+@[expose] public section
+
+open scoped ComplexOrder
+
+
+set_option linter.style.emptyLine false
+set_option linter.deprecated.module false
+
+universe u v
+
+section DiscreteLinearSystem
+
+/-!
+# Basic definitions for Discrete Linear Time Systems
+
+This module defines the state space representation of a discrete-time linear dynamical system.
+It includes the definition of the system state, the evolution function,
+and the property of satisfying the state equation.
+
+## Main Definitions
+* `DiscreteLinearSystemState`: Structure representing the system matrices (A and B),
+the current state, input, and initial state.
+* `DiscreteLinearSystemState.system_evolution`: Function computing the state at time `k`
+given an input sequence.
+* `DiscreteLinearSystemState.satisfies_state_equation`: Proposition stating that
+the sequence `x` satisfies the linear difference equation `x(k+1) = A x(k) + B u(k)`.
+
+
+## References
+
+https://en.wikipedia.org/wiki/State-space_representation
+https://www.cds.caltech.edu/~murray/books/AM08/pdf/fbs-public_24Jul2020.pdf
+
+-/
+
+
+
+variable {σ : Type u} {ι : Type v}
+variable [TopologicalSpace σ] [NormedAddCommGroup σ] [NormedSpace ℂ σ]
+variable [TopologicalSpace ι] [NormedAddCommGroup ι] [NormedSpace ℂ ι]
+variable [Inhabited ι]
+
+/-- Discrete-time linear dynamical system with state equation x(k+1) = A·x(k) + B·u(k). -/
+structure DiscreteLinearSystemState (σ : Type u) (ι : Type v)
+    [TopologicalSpace σ] [NormedAddCommGroup σ] [NormedSpace ℂ σ]
+    [TopologicalSpace ι] [NormedAddCommGroup ι] [NormedSpace ℂ ι] where
+  /-- State transition matrix (A), mapping the current state to the next state component (n×n). -/
+  a : σ →L[ℂ] σ
+  /-- Input matrix (B), mapping the current input to the next state component (n×p). -/
+  B : ι →L[ℂ] σ
+  /-- Initial state -/
+  x₀ : σ
+  /-- State sequence -/
+  x : ℕ → σ
+  /-- Input sequence -/
+  u : ℕ → ι
+
+variable {sys : DiscreteLinearSystemState σ ι}
+
+/-- System evolution function from initial state -/
+noncomputable def DiscreteLinearSystemState.system_evolution (u : ℕ → ι) : ℕ → σ
+  | 0 => sys.x₀
+  | k + 1 => sys.a (system_evolution u k) + sys.B (u k)
+
+/-- Discrete state space representation property -/
+def DiscreteLinearSystemState.satisfies_state_equation : Prop :=
+  ∀ k : ℕ, sys.x (k + 1) = sys.a (sys.x k) + sys.B (sys.u k)
+
+
+/-- Evolution from zero initial state with given input -/
+noncomputable def DiscreteLinearSystemState.evolve_from_zero
+   (u : ℕ → ι) (sys : DiscreteLinearSystemState σ ι) : ℕ → σ
+  | 0 => 0
+  | k + 1 => sys.a (evolve_from_zero u sys k) + sys.B (u k)
+
+
+
+
+/-- Zero input sequence -/
+def zero_input : ℕ → ι := fun _ => 0
+
+end DiscreteLinearSystem
diff --git a/Cslib/CPS/DiscreteLinearTime/Cayley.lean b/Cslib/CPS/DiscreteLinearTime/Cayley.lean
new file mode 100644
index 00000000..70502118
--- /dev/null
+++ b/Cslib/CPS/DiscreteLinearTime/Cayley.lean
@@ -0,0 +1,283 @@
+/-
+Copyright (c) 2026 Bashar Hamade. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors:Bashar Hamade
+-/
+
+module
+
+public import Cslib.Init
+public import Mathlib.Analysis.Normed.Module.Basic
+public import Mathlib.Analysis.Normed.Operator.ContinuousLinearMap
+public import Mathlib.LinearAlgebra.Span.Basic
+public import Mathlib.LinearAlgebra.FiniteDimensional.Defs
+public import Mathlib.Data.Finset.Basic
+public import Mathlib.Data.Complex.Basic
+public import Mathlib.Order.CompleteLattice.Basic
+public import Mathlib.Analysis.Normed.Module.Basic
+public import Mathlib.Data.Complex.Basic
+public import Mathlib.Analysis.Complex.Order
+public import Mathlib.Analysis.Normed.Field.Lemmas
+public import Mathlib.Analysis.Normed.Field.Basic
+public import Mathlib.Analysis.Complex.Exponential
+
+public import Mathlib.Analysis.Normed.Group.Basic
+public import Mathlib.Analysis.Normed.Module.RCLike.Real
+
+public import Mathlib.Analysis.Normed.Operator.ContinuousLinearMap
+public import Mathlib.Analysis.Complex.Order
+public import Mathlib.Analysis.Normed.Algebra.Spectrum
+public import Mathlib.Order.Basic
+public import Mathlib.LinearAlgebra.Charpoly.Basic
+
+public import Mathlib.LinearAlgebra.Matrix.Charpoly.LinearMap
+public import Mathlib.Algebra.Algebra.Bilinear
+
+
+@[expose] public section
+
+open scoped ComplexOrder
+
+
+/-!
+# Cayley-Hamilton implications for Controllability
+
+Auxiliary results based on the Cayley-Hamilton theorem for proving controllability results.
+
+
+## References
+https://www.statslab.cam.ac.uk/~rrw1/oc/L08.pdf
+
+-/
+
+set_option linter.style.emptyLine false
+set_option linter.deprecated.module false
+
+variable {σ ι : Type*}
+variable [TopologicalSpace σ] [NormedAddCommGroup σ] [NormedSpace ℂ σ]
+variable [TopologicalSpace ι] [NormedAddCommGroup ι] [NormedSpace ℂ ι]
+variable [Inhabited ι]
+
+
+
+/-- Power definition: A^(k+1) = A^k * A. -/
+lemma system_power_multiplication (a : σ →L[ℂ] σ) (k : ℕ) :
+    a ^ (k + 1) = (a ^ k).comp a := by
+  induction k with
+  | zero =>
+    simp only [zero_add, pow_one, pow_zero]
+    exact ContinuousLinearMap.id_comp a
+
+  | succ k ih =>
+    rw [pow_succ]
+    rw [ih]
+    rfl
+
+/-- Power definition commutation: A^(k+1) = A * A^k. -/
+lemma system_power_multiplication_flopped (a : σ →L[ℂ] σ) (k : ℕ) :
+    a ^ (k + 1) = a.comp (a^k) := by
+  induction k with
+  | zero =>
+    simp only [zero_add, pow_one, pow_zero]
+    exact ContinuousLinearMap.id_comp a
+
+  | succ k ih =>
+    rw [pow_succ]
+    rw [ih]
+
+    simp only [←ContinuousLinearMap.mul_def]
+    rw [mul_assoc]
+
+    congr 1
+
+
+
+/-- Commutativity helper: A (A^i B v) = A^(i+1) B v. -/
+lemma state_power_comm [FiniteDimensional ℂ σ]
+    (a : σ →L[ℂ] σ) (B : ι →L[ℂ] σ) (n : ℕ) (i : Fin n) (v : ι) :
+    a.toLinearMap ((a ^ i.val) (B v)) = (a ^ (i.val + 1)) (B v) := by
+
+    simp only [system_power_multiplication_flopped]
+    rfl
+
+
+
+/-- The degree of the characteristic polynomial equals the dimension of the space. -/
+lemma deg_eq_dim_of_space [FiniteDimensional ℂ σ]
+    (a : σ →L[ℂ] σ) (n : ℕ) (h_dim : Module.finrank ℂ σ = n) :
+    a.toLinearMap.charpoly.natDegree = n := by
+  rw [← h_dim]
+  exact LinearMap.charpoly_natDegree a.toLinearMap
+
+/-- By Cayley-Hamilton, A^n can be expressed as a linear combination of lower powers. -/
+lemma state_repr_as_lower_powers [FiniteDimensional ℂ σ]
+    (a : σ →L[ℂ] σ) (n : ℕ) (h_dim : Module.finrank ℂ σ = n) :
+    ∃ (c : Fin n → ℂ), a.toLinearMap ^ n = ∑ j : Fin n, c j • (a.toLinearMap ^ j.val) := by
+
+  let p := a.toLinearMap.charpoly
+
+  have ch : Polynomial.aeval a.toLinearMap p = 0 := LinearMap.aeval_self_charpoly _
+  have deg_p : p.natDegree = n := by rw [← h_dim]; exact LinearMap.charpoly_natDegree _
+  have monic_p : p.Monic := LinearMap.charpoly_monic _
+  have lead_coeff : p.coeff n = 1 := by rw [← deg_p]; exact monic_p.leadingCoeff
+
+  use fun j => -(p.coeff j.val)
+  simp only [neg_smul]
+
+
+  rw [Finset.sum_neg_distrib]
+
+  have expand : Polynomial.aeval a.toLinearMap p =
+      a.toLinearMap ^ n + ∑ j : Fin n, p.coeff j.val • a.toLinearMap ^ j.val := by
+    rw [Polynomial.aeval_eq_sum_range, deg_p, Finset.sum_range_succ]
+    rw [lead_coeff, one_smul, add_comm]
+    congr 1
+    exact (Fin.sum_univ_eq_sum_range (fun i => p.coeff i • a.toLinearMap ^ i) n).symm
+
+  rw [expand] at ch
+  exact eq_neg_of_add_eq_zero_left ch
+
+
+/-- The span of the controllability vectors is invariant under A. -/
+lemma controllabilityColumnSpace_invariant [FiniteDimensional ℂ σ]
+    (a : σ →L[ℂ] σ) (B : ι →L[ℂ] σ) (n : ℕ) (h_dim : Module.finrank ℂ σ = n) :
+    Submodule.map a.toLinearMap
+    (Submodule.span ℂ (⋃ i : Fin n, Set.range (fun v => (a ^ i.val) (B v)))) ≤
+    Submodule.span ℂ (⋃ i : Fin n, Set.range (fun v => (a ^ i.val) (B v))):= by
+
+
+
+
+
+  rw [Submodule.map_span_le]
+
+
+  intro x hx
+
+
+  simp only [Set.mem_iUnion] at hx
+  obtain ⟨i, hx⟩ := hx
+  simp only [Set.mem_range] at hx
+  obtain ⟨v, rfl⟩ := hx
+
+
+
+
+  by_cases h : i.val + 1 < n
+  · have : (a.toLinearMap ((a ^ i.val) (B v))) = (a ^ (i.val + 1)) (B v) := by
+      apply state_power_comm a B n i v
+
+    rw [this]
+    apply Submodule.subset_span
+    simp only [Set.mem_iUnion]
+    use ⟨i.val + 1, h⟩
+    simp only [Set.mem_range]
+    use v
+
+  · push_neg at h
+
+
+    have ch := LinearMap.aeval_self_charpoly a.toLinearMap
+
+
+    have deg_ch : a.toLinearMap.charpoly.natDegree = n := by
+      apply deg_eq_dim_of_space a  n h_dim
+
+
+    have i_eq : i.val = n - 1 := by
+      have : i.val < n := i.prop
+      omega
+
+
+
+
+
+    have : ∃ (c : Fin n → ℂ), a.toLinearMap ^ n = ∑ j : Fin n, c j • (a.toLinearMap ^ j.val) := by
+
+
+      apply state_repr_as_lower_powers a  n h_dim
+
+    obtain ⟨c, hc⟩ := this
+
+
+    have h_apply : (a ^ n) (B v) = ∑ j : Fin n, c j • ((a ^ j.val) (B v)) := by
+      have pow_eq : ∀ k x, (a ^ k) x = (a.toLinearMap ^ k) x := fun k x => by
+        induction k with
+        | zero => simp [pow_zero]
+        | succ k ih =>
+
+          simp [pow_succ', ih]
+
+
+
+      rw [pow_eq n]
+      rw [congrFun (congrArg DFunLike.coe hc) (B v)]
+      simp only [LinearMap.sum_apply, LinearMap.smul_apply]
+      congr 1
+      ext j
+      rw [← pow_eq j.val]
+
+
+    rw [i_eq]
+    have : a.toLinearMap ((a ^ (n - 1)) (B v)) = (a ^ n) (B v) := by
+      have hn : n = n - 1 + 1 := by omega
+      conv_rhs => rw [hn, pow_succ']
+      rfl
+
+    rw [this, h_apply]
+
+
+    apply Finset.sum_induction _ (· ∈ Submodule.span ℂ _)
+    · exact fun _ _ => Submodule.add_mem _
+    · exact Submodule.zero_mem _
+    · intro j _
+      apply Submodule.smul_mem
+      apply Submodule.subset_span
+      simp only [Set.mem_iUnion]
+      use j
+      simp only [Set.mem_range]
+      use v
+
+/-- Main Lemma: Higher powers A^j (j ≥ n) lie in the span of the first n powers. -/
+theorem cayley_hamilton_controllability' [FiniteDimensional ℂ σ]
+    (a : σ →L[ℂ] σ) (B : ι →L[ℂ] σ) (n : ℕ) (hn : n > 0)
+    (h_dim : Module.finrank ℂ σ = n) :
+    ∀ j ≥ n, ∀ v : ι, (a ^ j) (B v)
+    ∈ Submodule.span ℂ (⋃ i : Fin n, Set.range (fun v => (a ^ i.val) (B v))) := by
+  intro j hj v
+  induction j using Nat.strong_induction_on with
+  | _ j ih =>
+    by_cases hjn : j < n
+    · apply Submodule.subset_span
+      simp only [Set.mem_iUnion, Set.mem_range]
+
+
+      exact ⟨⟨j, hjn⟩, v, rfl⟩
+
+    · push_neg at hjn
+      by_cases hj_zero : j = 0
+      · omega
+      · have hj_pos : j > 0 := Nat.pos_of_ne_zero hj_zero
+        have : j = (j - 1) + 1 := by
+          omega
+
+        rw [this, pow_succ', ContinuousLinearMap.mul_apply]
+        apply controllabilityColumnSpace_invariant a B n h_dim
+        apply Submodule.mem_map_of_mem
+        have h_pred_ge : j - 1 ≥ n ∨ j - 1 < n := by
+          by_cases h : j - 1 ≥ n
+          · left; exact h
+          · right; push_neg at h; exact h
+
+
+        cases h_pred_ge with
+        | inl h_ge =>
+          apply ih
+          · omega
+          · exact h_ge
+        | inr h_lt =>
+
+
+          apply Submodule.subset_span
+          simp only [Set.mem_iUnion, Set.mem_range]
+          exact ⟨⟨j - 1, h_lt⟩, v, rfl⟩
diff --git a/Cslib/CPS/DiscreteLinearTime/Controllability.lean b/Cslib/CPS/DiscreteLinearTime/Controllability.lean
new file mode 100644
index 00000000..caf5e3d0
--- /dev/null
+++ b/Cslib/CPS/DiscreteLinearTime/Controllability.lean
@@ -0,0 +1,340 @@
+/-
+Copyright (c) 2026 Bashar Hamade. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Bashar Hamade
+-/
+
+module
+
+public import Cslib.Init
+public import Cslib.CPS.DiscreteLinearTime.Basic
+public import Cslib.CPS.DiscreteLinearTime.Reachability
+public import Cslib.CPS.DiscreteLinearTime.Cayley
+@[expose] public section
+universe u v
+
+variable {σ : Type u} {ι : Type v}
+variable [TopologicalSpace σ] [NormedAddCommGroup σ] [NormedSpace ℂ σ]
+variable [TopologicalSpace ι] [NormedAddCommGroup ι] [NormedSpace ℂ ι]
+variable [Inhabited ι]
+variable {sys : DiscreteLinearSystemState σ ι}
+
+open BigOperators Finset DiscreteLinearSystemState
+/-!
+# Controllability Analysis
+
+This module analyzes the controllability of discrete-time linear systems.
+It introduces the controllability matrix, the controllability column space,
+and connects these concepts to reachability.
+The main result is the Kalman Controllability Rank Condition.
+
+## Main Definitions
+* `controllabilityColumn`: The i-th term `A^i B`.
+* `controllabilityMatrix`: The matrix formed by columns `[B, AB, ..., A^(n-1)B]`.
+* `controllabilityColumnSpace`: The subspace spanned by the columns of the controllability matrix.
+* `totalReachableSubmodule`: The join of all reachable subspaces.
+
+## Main Theorems
+* `evolution_eq_matrix_form`: Expresses the state after `k` steps using the controllability matrix.
+* `controllabilityColumnSpace_stabilizes`: The span of the controllability matrix stabilizes
+after `n` steps (where `n` is the dimension of the state space), due to Cayley-Hamilton.
+* `full_finrank_equivalent_to_reachability`: A system is reachable if and only
+if the controllability matrix has full rank.
+
+
+## References
+https://www.statslab.cam.ac.uk/~rrw1/oc/L08.pdf
+https://en.wikipedia.org/wiki/Controllability
+
+-/
+
+
+/-- The i-th column block of the controllability matrix: Aⁱ B -/
+def controllabilityColumn (a : σ →L[ℂ] σ) (B : ι →L[ℂ] σ) (i : ℕ) : ι →L[ℂ] σ :=
+  (a ^ i).comp B
+
+/-- The controllability matrix up to step `n`, viewed as a function
+from index `i` to the operator `A^i B`. -/
+def controllabilityMatrix (a : σ →L[ℂ] σ) (B : ι →L[ℂ] σ) (n : ℕ) : Fin n → (ι →L[ℂ] σ) :=
+  fun i => (a ^ (i : ℕ)).comp B
+
+/-- The subspace spanned by the columns of the controllability matrix up to step `k`.
+    It represents the set of states reachable in `k` steps. -/
+def controllabilityColumnSpace (a : σ →L[ℂ] σ) (B : ι →L[ℂ] σ) (k : ℕ) : Submodule ℂ σ :=
+  Submodule.span ℂ (⋃ i : Fin k, Set.range (fun v => (a ^ i.val) (B v)))
+
+/-- Helper theorem: Expands the state evolution from zero initial state as a sum. -/
+theorem evolve_from_zero_eq_sum
+(s : DiscreteLinearSystemState σ ι) (u : ℕ → ι)  (k : ℕ) :
+   s.evolve_from_zero u k =
+     ∑ j ∈ Finset.range k, (s.a ^ (k - 1 - j)) (s.B (u j)) := by
+    induction k with
+    | zero =>
+      simp [DiscreteLinearSystemState.evolve_from_zero]
+    | succ k ih =>
+      simp only [evolve_from_zero, add_tsub_cancel_right]
+      rw [ih]
+      simp only [map_sum]
+      rw [Finset.sum_range_succ]
+      simp only [tsub_self, pow_zero, ContinuousLinearMap.one_apply, add_left_inj]
+      apply Finset.sum_congr rfl
+      intro x hx
+      have : s.a.comp (s.a ^ (k - 1 - x)) = s.a ^ (k -  x ) := by
+        rw [← system_power_multiplication_flopped]
+        congr
+        have : x < k := Finset.mem_range.mp hx
+        omega
+      rw [<-ContinuousLinearMap.comp_apply]
+      rw [this]
+
+/-- Matrix form of evolution equation -/
+theorem evolution_eq_matrix_form (s : DiscreteLinearSystemState σ ι) (kf : ℕ) (u : ℕ → ι)  :
+    s.evolve_from_zero u kf =
+    ∑ i : Fin kf, (controllabilityMatrix s.a s.B kf i) (u (kf - 1 - (i : ℕ))) := by
+    rw [evolve_from_zero_eq_sum]
+    simp only [controllabilityMatrix]
+    rw [← Finset.sum_attach]
+    simp only [ContinuousLinearMap.coe_comp', Function.comp_apply]
+    refine Finset.sum_bij'
+      (fun j _ => ⟨kf - 1 - j.val, ?_⟩)
+      (fun i _ => ⟨kf - 1 - i.val, ?_⟩)
+      ?_ ?_ ?_ ?_ ?_
+    · obtain ⟨j, hj⟩ := j; simp at hj ⊢; omega
+    · simp; omega
+    · intros; simp
+    · intros; simp
+    · intro ⟨j, hj⟩ _; ext; simp at hj ⊢; omega
+    · intro i _; ext; simp; omega
+    · intro ⟨j, hj⟩ _; simp only [mem_range] at hj ⊢; (congr; omega)
+
+
+/-- The set of reachable states in `k` steps exactly equals
+the controllability column space of order `k`. -/
+theorem reachable_set_eq_controllability_range
+(s : DiscreteLinearSystemState σ ι) (k : ℕ) (hk : k > 0) :
+    reachableSetInKSteps s k = controllabilityColumnSpace s.a s.B k := by
+  ext x
+  simp only [SetLike.mem_coe]
+  constructor
+  · intro h_reach
+    obtain ⟨u, h_reach⟩ := h_reach
+    rw [evolution_eq_matrix_form s k] at h_reach
+    rw [← h_reach]
+    apply Submodule.sum_mem
+    intro i _
+    apply Submodule.subset_span
+    simp only [Set.mem_iUnion, Set.mem_range]
+    use i
+    use u (k - 1 - ↑i)
+    simp [controllabilityMatrix]
+  · intro h_span
+    induction h_span using Submodule.span_induction with
+    | mem x hx =>
+      simp only [Set.mem_iUnion, Set.mem_range] at hx
+      obtain ⟨i, v, rfl⟩ := hx
+      use fun j => if j = k - 1 - i.val then v else 0
+      rw [evolve_from_zero_eq_sum]
+      have hi_lt : k - 1 - i.val < k := by omega
+      rw [Finset.sum_eq_single_of_mem (k - 1 - i.val) (Finset.mem_range.mpr hi_lt)]
+      · simp only [ite_true]
+        congr 1
+        have hi_bound : i.val < k := i.isLt
+        rw [Nat.sub_sub_self (Nat.lt_iff_le_pred hk |>.mp hi_bound)]
+      · intro b bh b_diff
+        rw [if_neg b_diff]
+        simp only [ContinuousLinearMap.map_zero]
+    | zero =>
+      use fun _ => 0
+      rw [evolve_from_zero_eq_sum]
+      simp only [ContinuousLinearMap.map_zero, Finset.sum_const_zero]
+    | add x y _ _ ihx ihy =>
+      obtain ⟨ux, hux⟩ := ihx
+      obtain ⟨uy, huy⟩ := ihy
+      use fun j => ux j + uy j
+      rw [evolve_from_zero_eq_sum]
+      simp only [ContinuousLinearMap.map_add, Finset.sum_add_distrib]
+      rw [← evolve_from_zero_eq_sum, ← evolve_from_zero_eq_sum]
+      rw [hux, huy]
+    | smul c x _ ihx =>
+      obtain ⟨ux, hux⟩ := ihx
+      use fun j => c • ux j
+      rw [evolve_from_zero_eq_sum]
+      simp only [ContinuousLinearMap.map_smul]
+      rw [← Finset.smul_sum]
+      congr 1
+      rw [evolve_from_zero_eq_sum] at hux
+      exact hux
+
+
+/-- The controllability column space is non-decreasing with respect to `k`. -/
+theorem controllabilityColumnSpace_mono (s : DiscreteLinearSystemState σ ι)
+{k₁ k₂ : ℕ} (h : k₁ ≤ k₂) :
+    controllabilityColumnSpace s.a s.B k₁ ≤ controllabilityColumnSpace s.a s.B k₂ := by
+  apply Submodule.span_mono
+  intro x hx
+  simp only [Set.mem_iUnion, Set.mem_range] at hx ⊢
+
+  obtain ⟨i, v, rfl⟩ := hx
+  exact ⟨⟨i.val, Nat.lt_of_lt_of_le i.isLt h⟩, v, rfl⟩
+
+
+/-- Stabilizatlon theorem: The controllability column space stabilizes after `n` steps,
+    where `n` is the dimension of the state space.
+    This is a consequence of the Cayley-Hamilton theorem. -/
+theorem controllabilityColumnSpace_stabilizes
+    [FiniteDimensional ℂ σ]
+    (s : DiscreteLinearSystemState σ ι) (n : ℕ) (hn : n > 0)
+    (h_dim : Module.finrank ℂ σ = n) :
+    ∀ k ≥ n, controllabilityColumnSpace s.a s.B k = controllabilityColumnSpace s.a s.B n := by
+  intro k hk
+  apply le_antisymm
+  · apply Submodule.span_le.mpr
+    intro x hx
+    simp only [Set.mem_iUnion, Set.mem_range] at hx
+    obtain ⟨i, v, rfl⟩ := hx
+    by_cases hi : i.val < n
+    · apply Submodule.subset_span
+      simp only [Set.mem_iUnion, Set.mem_range]
+      exact ⟨⟨i.val, hi⟩, v, rfl⟩
+    · push_neg at hi
+      obtain ⟨d, hd⟩ := Nat.exists_eq_add_of_le hi
+      induction d with
+      | zero =>
+        simp only [Nat.add_zero] at hd
+        rw [hd]
+        apply cayley_hamilton_controllability' s.a s.B n hn h_dim
+        trivial
+      | succ m ih =>
+        have h_ge : n + (m + 1) ≥ n := Nat.le_add_right n (m + 1)
+        apply cayley_hamilton_controllability' s.a s.B n hn h_dim
+        rw [hd]
+        trivial
+  · apply controllabilityColumnSpace_mono
+    trivial
+
+
+/-- The total reachable submodule is the limit (supremum) of the controllability column spaces. -/
+def totalReachableSubmodule  (sys : DiscreteLinearSystemState σ ι) : Submodule ℂ σ :=
+  ⨆ k : ℕ, controllabilityColumnSpace sys.a sys.B k
+
+
+/-- If the system is reachable, then the total reachable set is the entire space. -/
+theorem reachable_implies_total_reachable_eq_univ
+    (sys : DiscreteLinearSystemState σ ι)
+    (h_reach : sys.IsReachable) : totalReachableSet sys = Set.univ := by
+  ext x
+  simp only [totalReachableSet, Set.mem_iUnion]
+  simp only [Set.mem_univ, iff_true]
+  -- By definition of IsReachable, for any x there exists k and u such that we reach x
+  obtain ⟨k, u, hk_pos, hx⟩ := h_reach x
+  exact ⟨k, u, hx⟩
+
+
+/-- For a finite dimensional system of dimension `n`, the total reachable submodule is equal to
+    the controllability column space at step `n`. -/
+theorem totalReachableSubmodule_eq_controllabilityColumnSpace
+    [FiniteDimensional ℂ σ]
+    (sys : DiscreteLinearSystemState σ ι) (n : ℕ) (hn : n > 0)
+    (h_dim : Module.finrank ℂ σ = n) :
+    totalReachableSubmodule sys = controllabilityColumnSpace sys.a sys.B n := by
+  apply le_antisymm
+  · apply iSup_le
+    intro k
+    by_cases hk : k ≤ n
+    · exact controllabilityColumnSpace_mono sys hk
+    · push_neg at hk
+      rw [controllabilityColumnSpace_stabilizes sys n hn h_dim k (le_of_lt hk)]
+  · exact le_iSup (controllabilityColumnSpace sys.a sys.B) n
+
+/-- If a system is reachable, then its controllability column space
+at step `n` is the whole space (Top). -/
+theorem reachable_implies_controllabilityColumnSpace_eq_top
+    [FiniteDimensional ℂ σ]
+    (sys : DiscreteLinearSystemState σ ι) (n : ℕ) (hn : n > 0)
+    (h_dim : Module.finrank ℂ σ = n)
+    (h_reach : sys.IsReachable) :
+    controllabilityColumnSpace sys.a sys.B n = ⊤ := by
+  rw [← totalReachableSubmodule_eq_controllabilityColumnSpace sys n hn h_dim]
+  rw [eq_top_iff]
+  intro x _
+  obtain ⟨k, u, hk_pos, hx⟩ := h_reach x
+  have h_in_k : x ∈ reachableSetInKSteps sys k := ⟨u, hx⟩
+  rw [reachable_set_eq_controllability_range sys k hk_pos] at h_in_k
+  exact Submodule.mem_iSup_of_mem k h_in_k
+
+/-- Alias for `reachable_implies_controllabilityColumnSpace_eq_top`. -/
+theorem reachability_implies_full_rank
+    [FiniteDimensional ℂ σ]
+    (sys : DiscreteLinearSystemState σ ι) (n : ℕ) (hn : n > 0)
+    (h_dim : Module.finrank ℂ σ = n)
+    (h_reach : sys.IsReachable) :
+    controllabilityColumnSpace sys.a sys.B n = ⊤ :=
+  reachable_implies_controllabilityColumnSpace_eq_top sys n hn h_dim h_reach
+
+
+
+/-- Reachability implies that the dimension of the controllability subspace is equal
+to the dimension of the state space. -/
+theorem reachability_implies_rank_eq_dim
+    [FiniteDimensional ℂ σ]
+    (sys : DiscreteLinearSystemState σ ι)
+    (h_reach : sys.IsReachable)
+    (hn : Module.finrank ℂ σ > 0) :
+    Module.finrank ℂ (controllabilityColumnSpace sys.a sys.B (Module.finrank ℂ σ)) =
+    Module.finrank ℂ σ := by
+  have h_top := reachability_implies_full_rank sys (Module.finrank ℂ σ) hn rfl h_reach
+  rw [h_top]
+  exact finrank_top ℂ σ
+
+/-- Alternative formulation: The controllability submodule has full rank -/
+theorem reachability_implies_full_finrank
+    [FiniteDimensional ℂ σ]
+    (sys : DiscreteLinearSystemState σ ι) (n : ℕ)
+    (h_dim : Module.finrank ℂ σ = n)
+    (hn : n > 0)
+    (h_reach : sys.IsReachable) :
+    Module.finrank ℂ (controllabilityColumnSpace sys.a sys.B n) = n := by
+  have h_top := reachability_implies_full_rank sys n hn h_dim h_reach
+  rw [h_top]
+  rw [finrank_top]
+  exact h_dim
+
+
+
+/-- The Kalman Controllability Rank Condition: If the controllability matrix
+ has full rank, the system is reachable. -/
+theorem full_finrank_implies_reachability
+    [FiniteDimensional ℂ σ]
+    (sys : DiscreteLinearSystemState σ ι) (n : ℕ)
+    (h_dim : Module.finrank ℂ σ = n)
+    (hn : n > 0)
+    (h_top : Module.finrank ℂ (controllabilityColumnSpace sys.a sys.B n) = n) :
+    sys.IsReachable := by
+    intro x
+    have h_sub_top : controllabilityColumnSpace sys.a sys.B n = ⊤ := by
+      apply Submodule.eq_top_of_finrank_eq
+      rw [h_top, h_dim]
+    have h_reach_sys : reachableSetInKSteps sys n = controllabilityColumnSpace sys.a sys.B n := by
+      rw [reachable_set_eq_controllability_range]
+      exact hn
+    have hx_in : x ∈ reachableSetInKSteps sys n := by
+      rw [h_reach_sys]
+      rw [h_sub_top]
+      simp
+    obtain ⟨u, hu⟩ := hx_in
+    use n, u
+
+
+/-- Full equivalence: A system is reachable
+if and only if the controllability matrix has full rank. -/
+theorem full_finrank_equivalent_to_reachability
+    [FiniteDimensional ℂ σ]
+    (sys : DiscreteLinearSystemState σ ι) (n : ℕ)
+    (h_dim : Module.finrank ℂ σ = n)
+    (hn : n > 0) :
+    sys.IsReachable ↔ Module.finrank ℂ (controllabilityColumnSpace sys.a sys.B n) = n := by
+    constructor
+    · intro h_reach
+      exact reachability_implies_full_finrank sys n h_dim hn h_reach
+    · intro h_top
+      exact full_finrank_implies_reachability sys n h_dim hn h_top
diff --git a/Cslib/CPS/DiscreteLinearTime/Reachability.lean b/Cslib/CPS/DiscreteLinearTime/Reachability.lean
new file mode 100644
index 00000000..2e332179
--- /dev/null
+++ b/Cslib/CPS/DiscreteLinearTime/Reachability.lean
@@ -0,0 +1,54 @@
+/-
+Copyright (c) 2026 Bashar Hamade. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Bashar Hamade
+-/
+
+module
+
+public import Cslib.Init
+public import Cslib.CPS.DiscreteLinearTime.Basic
+
+@[expose] public section
+
+universe u v
+
+variable {σ : Type u} {ι : Type v}
+variable [TopologicalSpace σ] [NormedAddCommGroup σ] [NormedSpace ℂ σ]
+variable [TopologicalSpace ι] [NormedAddCommGroup ι] [NormedSpace ℂ ι]
+variable {sys : DiscreteLinearSystemState σ ι}
+
+
+/-!
+# Reachability Analysis
+
+This module defines the concept of reachability for discrete-time linear systems.
+It defines the set of states reachable in a specific number of steps (`reachableSetInKSteps`)
+and the total reachable set (`totalReachableSet`).
+
+## Main Definitions
+* `IsReachable`: Proposition that the system can reach any target state from the zero state.
+* `reachableSetInKSteps`: The set of states reachable in exactly `k` steps from zero.
+* `totalReachableSet`: The union of all reachable states over all time steps.
+
+## References
+http://www.dii.unimo.it/~zanasi/didattica/Teoria_dei_Sistemi/Luc_TDS_ING_2016_Reachability_and_Controllability.pdf
+https://www.cds.caltech.edu/~murray/books/AM08/pdf/fbs-public_24Jul2020.pdf
+-/
+
+/-- Reachability: For any target state x_f ∈ σ, there exists a positive integer k_f
+    and an input sequence U = (u[0], u[1], ..., u[k_f-1]) such that starting from
+    x[0] = 0, the system reaches x[k_f] = x_f -/
+def DiscreteLinearSystemState.IsReachable : Prop :=
+  ∀ x_f : σ, ∃ (k_f : ℕ) (u : ℕ → ι) , k_f > 0 ∧
+  DiscreteLinearSystemState.evolve_from_zero u sys k_f = x_f
+
+/-- The set of states reachable in exactly k steps -/
+def DiscreteLinearSystemState.reachableSetInKSteps
+(sys : DiscreteLinearSystemState σ ι) (k : ℕ) : Set σ :=
+  {x : σ | ∃ u : ℕ → ι, DiscreteLinearSystemState.evolve_from_zero u sys k = x}
+
+
+/-- The total reachable set is the union of reachable sets over all possible time steps `k`. -/
+def DiscreteLinearSystemState.totalReachableSet (sys : DiscreteLinearSystemState σ ι) : Set σ :=
+  ⋃ k : ℕ, reachableSetInKSteps sys k
diff --git a/Cslib/Foundations/Semantics/GameSemantics/Basic.lean b/Cslib/Foundations/Semantics/GameSemantics/Basic.lean
new file mode 100644
index 00000000..965c8019
--- /dev/null
+++ b/Cslib/Foundations/Semantics/GameSemantics/Basic.lean
@@ -0,0 +1,397 @@
+/-
+Copyright (c) 2026 Bashar Hamade. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Bashar Hamade
+-/
+module
+public import Cslib.Foundations.Semantics.LTS.Basic
+
+@[expose] public section
+
+/-!
+# Reachability Games
+
+This module formalizes reachability games, a fundamental concept in game theory with
+applications in formal verification, concurrency theory, and logic.
+
+A reachability game is played on a graph (represented as an LTS) between two players:
+- **Defender** (also called Player 0), who wins if the play reaches a designated target set
+- **Attacker** (also called Player 1), who tries to avoid the target set forever
+
+These games are particularly useful for specifying and verifying properties of concurrent
+systems, such as safety and liveness properties.
+
+## Main definitions
+
+- `ReachabilityGame` is a structure defining the game graph, consisting of positions,
+  transition labels, a target set of positions for the Defender, an LTS representing
+  possible moves, and an initial position.
+
+- `Attacker` is the set of positions belonging to the Attacker (complement of the Defender's
+  target set).
+
+- `step` formalizes the existence of a valid move between two positions.
+
+- `Play` is a structure representing a (potentially infinite) play in the game, with
+  fields ensuring the play starts at the initial position and follows valid transitions.
+
+- `DefenderWins` and `AttackerWins` define winning conditions for the respective players.
+
+- `AttackerStrategy` and `DefenderStrategy` represent deterministic strategies for each player.
+
+- `AttackerWinsFromPos` and `AttackerWinRegion` define the set of positions from which
+  the Attacker can force a win.
+
+## Main statements
+
+- `disjoint_Gd_Attacker`: The Defender's target set and the Attacker positions are disjoint.
+
+- `union_Gd_Attacker`: Every position is either in the Defender's target set or belongs to
+  the Attacker.
+
+- `finite_or_infinite`: Every play is either finite or infinite.
+
+- `not_both_finite_infinite`: A play cannot be both finite and infinite.
+
+- `Defender_or_Attacker_win`: Exactly one player wins each play.
+
+- `AttackerStrategy.isDeterministic` and `DefenderStrategy.isDeterministic`:
+  Characterizations of deterministic strategies.
+
+## References
+
+  * [B. Bisping, D. N. Jansen, and U. Nestmann,
+  *Deciding All Behavioral Equivalences at Once:
+  A Game for Linear-Time–Branching-Time Spectroscopy*][BispingEtAl2022]
+
+-/
+
+namespace Cslib
+
+universe u v
+
+/--
+A ReachabilityGame consists of:
+- `Label`: The type of transition labels
+- `Pos`: The type of positions in the game graph
+- `G_d`: The set of positions designated as the Defender's target set
+- `GameMoves`: An LTS defining the valid moves in the game
+- `g_0`: The initial position of the game
+
+The game is played between a Defender (who tries to reach `G_d`) and an Attacker
+(who tries to avoid `G_d` forever).
+-/
+structure Game where
+  /-- The type of transition labels. -/
+  Label : Type v
+  /-- The type of positions in the game graph. -/
+  Pos : Type u
+  /-- The set of positions constituting the Defender's target set. -/
+  G_d : Set Pos
+  /-- The labelled transition system defining valid moves. -/
+  GameMoves : LTS Pos Label
+  /-- The initial position of the game. -/
+  g_0 : Pos
+
+namespace Game
+
+variable {G : Game}
+
+/--
+The set of positions belonging to the Attacker.
+A position belongs to the Attacker if it is not in the Defender's target set `G_d`.
+-/
+def Attacker : Set G.Pos := { p : G.Pos | ¬G.G_d p }
+
+/-- The Defender's target set and the Attacker positions are disjoint. -/
+theorem disjoint_Gd_Attacker : G.G_d ∩ G.Attacker = ∅ := by
+  ext x; simp only [Attacker, Set.mem_inter_iff, Set.mem_setOf_eq, Set.mem_empty_iff_false,
+    iff_false, not_and, not_not]
+  intro x_in_Gd
+  trivial
+
+/-- Every position is either in the Defender's target set or belongs to the Attacker. -/
+theorem union_Gd_Attacker : G.G_d ∪ G.Attacker = Set.univ := by
+  ext x;
+  simp only [Attacker, Set.mem_union, Set.mem_setOf_eq, Set.mem_univ, iff_true]; apply Classical.em
+
+/--
+A step exists from position `p` to position `p'` if there is a transition
+in the game LTS from `p` to `p'` with some label.
+-/
+def step (p p' : G.Pos) : Prop := ∃ (l : G.Label), G.GameMoves.Tr p l p'
+
+/-!
+## Plays
+
+This section defines plays, which are (potentially infinite) sequences of positions
+in a reachability game.
+-/
+
+/--
+A Play is a (potentially infinite) sequence of positions in the game.
+
+A play consists of:
+- `π`: A partial function from natural numbers to positions, representing the
+  position at each step of the play
+- `start_eq_initial`: The play must start at the initial position `g_0`
+- `continuous`: If a position exists at step n+1, then a position must exist at step n
+- `coherent`: If a position exists at step n, all earlier positions must exist
+- `is_path`: The sequence must follow valid transitions in the game
+-/
+structure Play where
+  /-- Position at index i (partial function). -/
+  π : ℕ → Option G.Pos
+  /-- First position exists and equals `g_0`. -/
+  start_eq_initial : π 0 = some G.g_0
+  /-- Continuity: if position n+1 exists, then position n must exist. -/
+  continuous : ∀ n, π (n + 1) ≠ none → π n ≠ none
+  /-- Coherent: if position n exists, all earlier positions exist. -/
+  coherent : ∀ n m, m < n → π n ≠ none → π m ≠ none
+  /-- Path property: valid transitions between consecutive positions. -/
+  is_path : ∀ n p p', π n = some p → π (n + 1) = some p' → step p p'
+
+namespace Play
+
+
+open Classical in
+/--
+The length of a finite play, or `none` if the play is infinite.
+If the play terminates (some `π n = none`), returns `some (Nat.find h - 1)`,
+which is the last index where a position exists.
+-/
+noncomputable def len (play : @Play G) : Option ℕ :=
+  if h : ∃ n, play.π n = none
+  then some (Nat.find h - 1)
+  else none
+
+/-!
+### Finite and Infinite Plays
+
+This section defines properties for classifying plays as finite or infinite.
+-/
+
+/-- A play is finite if there exists a position at step n but no position at step n+1. -/
+def isFinite (play : @Play G) : Prop := ∃ n, play.π n ≠ none ∧ play.π (n + 1) = none
+
+/-- A play is infinite if for every n, whenever a position exists at step n,
+a position also exists at step n+1. -/
+def isInfinite (play : @Play G) : Prop := ∀ n, play.π n ≠ none → play.π (n + 1) ≠ none
+
+/-- The initial position π 0 always exists (never none). -/
+lemma π_zero_not_none (play : @Play G) : play.π 0 ≠ none := by
+  rw [play.start_eq_initial]
+  simp
+
+/-- Every play is either finite or infinite. -/
+theorem finite_or_infinite (play : @Play G) : play.isFinite ∨ play.isInfinite := by
+  by_cases h : ∃ n, play.π n = none
+  · -- Case: There exists some n where π n = none
+    left
+    unfold isFinite
+    -- Use Nat.find to get the minimal such n
+    let n := Nat.find h
+    have hn : play.π n = none := Nat.find_spec h
+     -- Show that n > 0 (because π 0 ≠ none)
+    have n_pos : n > 0 := by
+      by_contra h_neg
+      simp only [gt_iff_lt, not_lt, nonpos_iff_eq_zero] at h_neg
+      have : play.π 0 = none := by rwa [h_neg] at hn
+      exact π_zero_not_none play this
+    -- So n - 1 exists and is a natural number
+    use (n - 1)
+    constructor
+    · -- Show π (n-1) ≠ none
+      have : n - 1 < n := by
+        simp only [tsub_lt_self_iff, Nat.lt_add_one, and_true]
+        trivial
+      exact Nat.find_min h this
+    · -- Show π n = none (which equals π ((n-1)+1) when n > 0)
+      have : n = (n - 1) + 1 := by omega
+      rwa [← this]
+  · -- Case: For all n, π n ≠ none
+    right
+    unfold isInfinite
+    push_neg at h
+    exact fun n _ => h (n + 1)
+
+/-- A play cannot be both finite and infinite (these properties are mutually exclusive). -/
+theorem not_both_finite_infinite (play : @Play G) : ¬(play.isFinite ∧ play.isInfinite) := by
+  intro ⟨hfin, hinf⟩
+  obtain ⟨n, hn_not_none, hn1_none⟩ := hfin
+  have : play.π (n + 1) ≠ none := hinf n hn_not_none
+  contradiction
+
+
+open Classical in
+/-- For a finite play, returns the index of the first `none` value in π. -/
+def final_pos_index (play : @Play G) (h : play.isFinite) : ℕ :=
+  -- Convert isFinite to the existential Nat.find expects
+  let h' : ∃ n, play.π n = none := by
+    obtain ⟨n, _, hn1⟩ := h
+    exact ⟨n + 1, hn1⟩
+  Nat.find h'
+
+/-!
+### Winning Conditions
+
+This section defines winning conditions for the Defender and Attacker.
+-/
+
+/--
+The Defender wins a play if either:
+1. The play is infinite (Attacker cannot force termination), or
+2. The play is finite and ends in a position belonging to `G_d` (the Defender's target set)
+-/
+def DefenderWins (play : @Play G) : Prop :=
+  play.isInfinite ∨ (∃ h : play.isFinite,
+    ∃ p, play.π (final_pos_index play h - 1) = some p ∧ p ∈ Attacker)
+
+
+/-- The Attacker wins a play if the Defender does not win. -/
+def AttackerWins (play : @Play G) : Prop :=
+  ¬ (DefenderWins play)
+
+/-- Exactly one player wins each play (excluded middle for winning conditions). -/
+theorem Defender_or_Attacker_win (play : @Play G) : DefenderWins play ∨ AttackerWins play := by
+  have fin_or_infin := finite_or_infinite play
+  cases fin_or_infin with
+  | inl hfin =>
+    -- Finite case
+    by_cases h : ∃ p, play.π (final_pos_index play hfin - 1) = some p ∧ p ∈ Attacker
+    · -- Defender wins
+      left
+      right
+      exact ⟨hfin, h⟩
+    · -- Attacker wins
+      right
+      unfold AttackerWins DefenderWins
+      push_neg
+      constructor
+      · unfold isInfinite
+        push_neg
+        exact hfin
+      · intro h_fin p exists_some_p
+        push_neg at h
+        specialize h p
+        apply h
+        trivial
+  | inr hinf =>
+    -- Infinite case: Defender wins
+    left
+    left
+    exact hinf
+
+
+/-!
+### Strategies
+
+This section defines strategies for the Attacker and Defender players.
+-/
+
+/--
+An AttackerStrategy is a function that, given an Attacker position `g` and proof that
+`g` belongs to the Attacker, returns a next position `g'` along with a proof that
+`step g g'` holds (i.e., there is a valid move from `g` to `g'`).
+-/
+def AttackerStrategy :=
+  (g : G.Pos) → g ∈ Attacker → { g' : G.Pos // G.step g g' }
+
+/--
+A DefenderStrategy is a function that, given a Defender position `g` (i.e., `g ∈ G_d`)
+and proof that `g` belongs to the Defender's target set, returns a next position `g'`
+along with a proof that `step g g'` holds.
+-/
+def DefenderStrategy :=
+  (g : G.Pos) → g ∈ G.G_d → { g' : G.Pos // G.step g g' }
+
+/-!
+#### Deterministic Strategies
+
+This section characterizes when strategies are deterministic.
+-/
+
+/-- An AttackerStrategy is deterministic if it always selects the same next position
+for the same input position (i.e., it is a function, not a relation). -/
+def AttackerStrategy.isDeterministic (strategy : @AttackerStrategy G) : Prop :=
+  ∀ (g : G.Pos) (hg : g ∈ Attacker),
+    ∀ (g1 g2 : { g' : G.Pos // step g g' }),
+    strategy g hg = g1 → strategy g hg = g2 → g1.val = g2.val
+
+
+/-- A DefenderStrategy is deterministic if it always selects the same next position
+for the same input position (i.e., it is a function, not a relation). -/
+def DefenderStrategy.isDeterministic (strategy : @DefenderStrategy G) : Prop :=
+  ∀ (g : G.Pos) (hg : g ∈ G.G_d),
+    ∀ (g1 g2 : { g' : G.Pos // step g g' }),
+    strategy g hg = g1 → strategy g hg = g2 → g1.val = g2.val
+
+/--
+A play follows an AttackerStrategy if, whenever the play is at an Attacker position,
+the next position is the one prescribed by the strategy.
+-/
+def followsAttackerStrategy (play : @Play G) (σ : @AttackerStrategy G) : Prop :=
+  ∀ n p, play.π n = some p →
+    (hg : p ∈ Attacker) →
+    play.π (n + 1) = some (σ p hg).val
+
+/-!
+### Winning Regions
+
+This section defines the set of positions from which the Attacker can force a win.
+-/
+
+/--
+AttackerWinsFromPos pos means the Attacker has a strategy that guarantees winning
+from position `pos`, no matter how the Defender responds.
+-/
+def AttackerWinsFromPos (pos : G.Pos) : Prop :=
+  ∃ (strategy : @AttackerStrategy G),
+    ∀ (play : @Play G),
+      play.π 0 = some pos →
+      play.followsAttackerStrategy strategy →
+      play.AttackerWins
+
+/--
+The AttackerWinRegion is the set of Attacker positions from which the Attacker
+can force a win (i.e., positions in the winning region for the Attacker).
+-/
+def AttackerWinRegion : Set G.Pos :=
+  { g : G.Pos | AttackerWinsFromPos g }
+
+
+/--
+A play follows a DefenderStrategy if, whenever the play is at a Defender position,
+the next position is the one prescribed by the strategy.
+-/
+def followsDefenderStrategy (play : @Play G) (σ : @DefenderStrategy G) : Prop :=
+  ∀ n p, play.π n = some p →
+    (hg : p ∈ G.G_d) →
+    play.π (n + 1) = some (σ p hg).val
+
+/--
+DefenderWinsFromPos pos means the Defender has a strategy that guarantees winning
+from position pos, no matter how the Attacker responds.
+-/
+def DefenderWinsFromPos (pos : G.Pos) : Prop :=
+  ∃ (strategy : @DefenderStrategy G),
+  ∀ (play : @Play G),
+    play.π 0 = some pos →
+    play.followsDefenderStrategy strategy →
+    play.DefenderWins
+
+/--
+The DefenderWinRegion is the set of positions from which the Defender
+can force a win.
+-/
+def DefenderWinRegion : Set G.Pos :=
+  { g : G.Pos | DefenderWinsFromPos g }
+
+
+
+
+end Play
+
+end Game
+
+end Cslib
diff --git a/Cslib/Foundations/Semantics/GameSemantics/HMLGame.lean b/Cslib/Foundations/Semantics/GameSemantics/HMLGame.lean
new file mode 100644
index 00000000..8acc0228
--- /dev/null
+++ b/Cslib/Foundations/Semantics/GameSemantics/HMLGame.lean
@@ -0,0 +1,326 @@
+/-
+Copyright (c) 2025 Bashar Hamade. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Bashar Hamade
+-/
+module
+public import Cslib.Init
+public import Cslib.Foundations.Semantics.LTS.Basic
+public import Cslib.Logics.HennessyMilnerLogic.Basic
+public import Cslib.Foundations.Semantics.GameSemantics.Basic
+
+@[expose] public section
+
+
+/-!
+# HML Games
+
+This module formalizes HML games, a game-theoretic characterization of
+Hennessy–Milner Logic (HML) on Labelled Transition Systems (LTS).
+
+HML games provide an alternative, operational characterization of modal
+logic satisfaction. They are played between two players:
+- **Defender** (∃ player): tries to prove that a state satisfies a formula
+- **Attacker** (∀ player): tries to refute the formula
+
+A state satisfies an HML formula if and only if the Defender has a winning
+strategy in the corresponding HML game starting from that state.
+
+## Game Structure
+
+For a transition system S = (P, Σ, →), the HML game G^S_HML[g₀] = (G, G_d, ↝, g₀)
+is played on G = P × HML[Σ], where positions are pairs (p, φ) consisting of
+a process p and a formula φ.
+
+## Main definitions
+
+- `HMLGame` is a structure defining the game components: labels, processes,
+  transition system, and initial position.
+
+- `Gd` defines the set of defender-controlled positions, corresponding to
+  modal observations ⟨a⟩φ and negated conjunctions ¬⋀ᵢφᵢ.
+
+- `Ga` defines the set of attacker-controlled positions (complement of Gd).
+
+- `move_defender` formalizes valid moves from defender positions.
+
+- `move_attacker` formalizes valid moves from attacker positions.
+
+- `move` is the union of defender and attacker moves.
+
+- `toGame` converts an HMLGame to a generic Game structure.
+
+- `defenderWinsHML` gives a direct recursive characterization of when
+  the defender wins.
+
+## Main statements
+
+- `game_semantics_soundness`: Proves that the game-theoretic characterization
+  of HML satisfaction is equivalent to the standard Tarski-style semantics.
+  This establishes the fundamental connection between games and modal logic.
+
+## References
+
+* [M. Hennessy and R. Milner,
+*Algebraic Laws for Non-Determinism and Concurrency*][HennessyMilner1985]
+
+* [B. Bisping, D. N. Jansen, and U. Nestmann,
+*Deciding All Behavioral Equivalences at Once:
+A Game for Linear-Time–Branching-Time Spectroscopy*][BispingEtAl2022]
+
+-/
+
+namespace Cslib
+
+open HennessyMilner
+
+universe u v
+
+variable {State : Type u} {Label : Type v}
+
+/--
+HML game structure: For a transition system S = (P, Σ, →), the HML game
+G^S_HML[g₀] = (G, G_d, ↝, g₀) is played on G = P × HML[Σ].
+
+A position is a pair (p, φ) where p is a process and φ is an HML formula.
+The game is played between a Defender (who tries to prove φ holds at p) and
+an Attacker (who tries to refute it).
+-/
+structure HMLGame where
+  /-- The type of transition labels (alphabet Σ). -/
+  Label : Type v
+  /-- The type of processes/states in the transition system. -/
+  Process : Type u
+  /-- The labelled transition system S = (P, Σ, →). -/
+  S : LTS Process Label
+  /-- The initial position g₀ = (p, φ). -/
+  g0 : Process × Formula Label
+
+namespace HMLGame
+
+variable {G : HMLGame}
+
+/-!
+## Game Positions and Moves
+
+This section defines the key components of the HML game: the sets of
+positions controlled by each player and the valid moves they can make.
+-/
+
+/--
+The set of defender-controlled positions (G_d).
+
+The Defender controls positions where the formula requires an existential
+interpretation:
+- Modal observation ⟨a⟩φ: Defender must find an a-transition
+- Negated conjunction ¬⋀ᵢφᵢ: Defender must pick one disjunct to defend
+
+Positions in G_d are winning for the Defender if they reach a true formula.
+-/
+def Gd : Set (G.Process × Formula G.Label) :=
+  { pos | match pos.2 with
+    | .modal _ _ => True
+    | .neg (.conj _) => True
+    | _ => False }
+
+/--
+Alternative definition of defender-controlled positions using explicit
+set-builder notation with existential quantification.
+-/
+def Gd' : Set (G.Process × Formula G.Label) :=
+  { pos | ∃ (a : G.Label) (φ : Formula G.Label), pos.2 = .modal a φ } ∪
+  { pos | ∃ (φs : List (Formula G.Label)), pos.2 = .neg (.conj φs) }
+
+/--
+The set of attacker-controlled positions (G_a), defined as the complement
+of the defender's positions.
+
+The Attacker controls positions where the formula requires a universal
+interpretation:
+- Conjunction ⋀ᵢφᵢ: Attacker can choose any conjunct to challenge
+- Negated modal ¬⟨a⟩φ: Attacker must find a transition to refute
+- Double negation ¬¬φ: Attacker can eliminate one negation
+-/
+def Ga : Set (G.Process × Formula G.Label) :=
+  { pos | ¬ Gd pos }
+
+/-!
+### Defender's Moves
+
+This section formalizes the valid moves available to the Defender player.
+-/
+
+/--
+The Defender's move relation from position x to position y.
+
+At a modal observation ⟨a⟩φ, the Defender chooses a transition p -a→ p'.
+The next position is (p', φ) at the same formula.
+
+At a negated conjunction ¬⋀ᵢφᵢ, the Defender chooses one disjunct φ'ᵢ
+(from the negated list of conjuncts) to defend. The process stays the same.
+-/
+def move_defender : (G.Process × Formula G.Label) → (G.Process × Formula G.Label) → Prop
+  | (p, .modal a φ), (p', φ') =>
+      G.S.Tr p a p' ∧ φ' = φ
+  | (p, .neg (.conj φs)), (p', φ') =>
+      p' = p ∧ φ' ∈ φs.map Formula.neg
+  | _, _ => False
+
+/-!
+### Attacker's Moves
+
+This section formalizes the valid moves available to the Attacker player.
+-/
+
+/--
+The Attacker's move relation from position x to position y.
+
+At a negated modal ¬⟨a⟩φ, the Attacker chooses a transition p -a→ p'.
+The next position is (p', ¬φ).
+
+At a conjunction ⋀ᵢφᵢ, the Attacker chooses one conjunct φᵢ to challenge.
+The process stays the same.
+
+At a double negation ¬¬φ, the Attacker can eliminate one negation,
+moving to (p, φ).
+-/
+def move_attacker : (G.Process × Formula G.Label) → (G.Process × Formula G.Label) → Prop
+  | (p, .neg (.modal a φ)), (p', φ') =>
+      G.S.Tr p a p' ∧ φ' = .neg φ
+  | (p, .conj φs), (p', φ') =>
+      p' = p ∧ φ' ∈ φs
+  | (p, .neg (.neg φ)), (p', φ') =>
+      p' = p ∧ φ' = φ
+  | _, _ => False
+
+/--
+The complete move relation for the HML game, defined as the union of
+defender and attacker moves.
+
+This captures all valid transitions between positions in the game graph.
+-/
+def move (x y : G.Process × Formula G.Label) : Prop :=
+  move_defender x y ∨ move_attacker x y
+
+/--
+Converts an HMLGame to the generic Game structure defined in Basic.lean.
+
+This allows HML games to be analyzed using the general game theory
+framework (plays, strategies, winning regions, etc.).
+-/
+def toGame (HG : HMLGame) : Game where
+  Label := HG.Label
+  Pos := HG.Process × Formula HG.Label
+  G_d := HG.Gd
+  GameMoves := { Tr := fun p _ p' => HG.move p p' }
+  g_0 := HG.g0
+
+/-!
+## Direct Characterization of Winning
+
+This section provides a direct (non-game-theoretic) characterization of
+when the Defender wins the HML game, by recursion on the formula structure.
+-/
+
+/--
+Direct characterization of Defender's winning positions in an HML game.
+
+This function gives a recursive characterization equivalent to "the Defender
+has a winning strategy", but defined directly by recursion on the formula:
+- ⊤: Defender trivially wins (true formula)
+- ⊥: Defender cannot win (false formula)
+- ⟨a⟩φ: Defender wins if there exists an a-transition to a state where they win
+- ⋀ᵢφᵢ: Defender wins if they win on all conjuncts
+- ¬φ: Defender wins if they do NOT win on φ (Attacker loses)
+-/
+def defenderWinsHML (lts : LTS State Label) : State → Formula Label → Prop
+  | _, .true => True
+  | _, .false => False
+  | p, .modal a φ => ∃ p', lts.Tr p a p' ∧ defenderWinsHML lts p' φ
+  | p, .conj φs => ∀ φ ∈ φs, defenderWinsHML lts p φ
+  | p, .neg φ => ¬defenderWinsHML lts p φ
+
+/-!
+## Correctness Theorem
+
+This section proves the fundamental correctness result: the game-theoretic
+characterization is equivalent to the standard Tarski semantics.
+-/
+
+/--
+**Main Correctness Theorem**: The game-theoretic characterization of HML
+satisfaction is equivalent to the standard Tarski-style semantics.
+
+That is, for any LTS lts, state p, and HML formula φ:
+  defenderWinsHML lts p φ  ↔  lts ⊨ p φ
+
+This theorem establishes that the Defender has a winning strategy in the
+HML game starting from (p, φ) if and only if p satisfies φ according to
+the standard modal logic semantics.
+
+The proof proceeds by structural induction on the formula φ, showing that
+each case of the recursive definition of `defenderWinsHML` corresponds
+exactly to the semantic clause for that formula constructor.
+-/
+theorem game_semantics_soundness (lts : LTS State Label) (p : State) (φ : Formula Label) :
+    defenderWinsHML lts p φ ↔ satisfies lts p φ := by
+  induction φ using Formula.ind_on generalizing p with
+  | h_true =>
+    -- Base Case 1: φ = ⊤
+    simp [defenderWinsHML, satisfies]
+  | h_false =>
+    -- Base Case 2: φ = ⊥
+    simp [defenderWinsHML, satisfies]
+  | h_modal a ψ ih =>
+      constructor
+      · intro h_def_win
+        simp only [defenderWinsHML] at h_def_win
+        obtain ⟨p', h_tr, h_def_win_ψ⟩ := h_def_win
+        simp only [satisfies]
+        use p'
+        constructor
+        · exact h_tr
+        · specialize ih p'
+          exact ih.mp h_def_win_ψ
+      · simp only [satisfies]
+        intro ⟨p', h_tr, h_sat_ψ⟩
+        simp only [defenderWinsHML]
+        use p'
+        constructor
+        · exact h_tr
+        · specialize ih p'
+          exact ih.mpr h_sat_ψ
+  | h_conj φs ih_list =>
+      constructor
+      · intro h_def_win
+        simp only [defenderWinsHML] at h_def_win
+        simp only [satisfies]
+        intro φ h_mem
+        specialize ih_list φ h_mem p
+        exact ih_list.mp (h_def_win φ h_mem)
+      · simp only [satisfies]
+        intro h_sat_all
+        simp only [defenderWinsHML]
+        intro φ h_mem
+        specialize ih_list φ h_mem p
+        exact ih_list.mpr (h_sat_all φ h_mem)
+  | h_neg ψ ih =>
+      constructor
+      · intro h_def_win
+        simp only [defenderWinsHML] at h_def_win
+        simp only [satisfies]
+        intro h_sat_ψ
+        specialize ih p
+        have h_def_ψ := ih.mpr h_sat_ψ
+        contradiction
+      · simp only [satisfies]
+        intro h_not_sat_ψ
+        simp only [defenderWinsHML]
+        intro h_def_ψ
+        specialize ih p
+        have h_sat_ψ := ih.mp h_def_ψ
+        contradiction
+
+end HMLGame
+
+end Cslib
diff --git a/Cslib/Logics/HennessyMilnerLogic/Basic.lean b/Cslib/Logics/HennessyMilnerLogic/Basic.lean
new file mode 100644
index 00000000..15161b1b
--- /dev/null
+++ b/Cslib/Logics/HennessyMilnerLogic/Basic.lean
@@ -0,0 +1,195 @@
+module
+
+
+public import Cslib.Init
+public import Cslib.Foundations.Semantics.LTS.Basic
+@[expose] public section
+
+
+/-! # Hennessy–Milner Logic
+
+Hennessy–Milner logic (HML) is a modal logic used for reasoning about the behaviour of
+Labelled Transition Systems (LTS). It was introduced by Matthew Hennessy and Robin Milner
+to capture notion of observational equivalence in concurrent systems.
+
+## Main definitions
+
+- `Formula` defines the inductive syntax of HML formulas over an alphabet Σ
+- `satisfies` provides the semantics for evaluating formulas on LTS states
+- Scoped notation for modal operators (⟨a⟩ and ¬)
+
+## References
+
+* [M. Hennessy and R. Milner,
+*Algebraic Laws for Non-Determinism and Concurrency*][HennessyMilner1985]
+* [B. Bisping, D. N. Jansen, and U. Nestmann,
+*Deciding All Behavioral Equivalences at Once:
+A Game for Linear-Time–Branching-Time Spectroscopy*][BispingEtAl2022]
+
+-/
+
+namespace Cslib
+
+universe u v
+
+variable {State : Type u} {Label : Type v}
+
+namespace HennessyMilner
+
+/--
+Hennessy–Milner logic formulas over an alphabet Σ (Label type).
+
+(Hennessy–Milner logic): Given an alphabet Σ, the syntax of
+Hennessy-Milner logic formulas, HML[Σ], is inductively defined as follows:
+1. **Observations:** If φ ∈ HML[Σ] and a ∈ Σ, then ⟨a⟩φ ∈ HML[Σ]
+2. **Conjunctions:** If φ_i ∈ HML[Σ] for all i from an index set I, then ⋀_{i∈I} φ_i ∈ HML[Σ]
+3. **Negations:** If φ ∈ HML[Σ], then ¬φ ∈ HML[Σ]
+-/
+inductive Formula (Label : Type v) : Type v where
+  /-- The constant truth formula. -/
+  | true : Formula Label
+  /-- The constant false formula. -/
+  | false : Formula Label
+  /-- Modal observation operator: ⟨a⟩φ means
+  "there exists an a-transition to a state satisfying φ". -/
+  | modal (a : Label) (φ : Formula Label) : Formula Label
+  /-- Conjunction of a finite list of formulas. -/
+  | conj (φs : List (Formula Label)) : Formula Label
+  /-- Negation operator: ¬φ means "not φ". -/
+  | neg (φ : Formula Label) : Formula Label
+
+/-- Disjunction of a finite list of formulas (derived connective). -/
+def Formula.disj (φs : List (Formula Label)) : Formula Label :=
+  Formula.neg (Formula.conj (φs.map Formula.neg))
+
+/-- Implication φ ⊃ ψ (derived connective). -/
+def Formula.impl (φ ψ : Formula Label) : Formula Label :=
+  Formula.disj [Formula.neg φ, ψ]
+
+/-- Semantic satisfaction relation for HML formulas on LTS states. -/
+def satisfies (lts : LTS State Label) (s : State) : Formula Label → Prop
+  | .true => True
+  | .false => False
+  | .modal a φ => ∃ s', lts.Tr s a s' ∧ satisfies lts s' φ
+  | .conj φs => ∀ φ ∈ φs, satisfies lts s φ
+  | .neg φ => ¬satisfies lts s φ
+
+/-- Shorthand for the satisfaction relation. -/
+scoped infix:50 " ⊨ " => satisfies
+
+/-- A state satisfies a formula if there exists some transition leading
+to a state that satisfies it. -/
+theorem satisfies_modal_intro (lts : LTS State Label) (s : State) (a : Label) (φ : Formula Label) :
+  ∀ s', lts.Tr s a s' → satisfies lts s' φ → satisfies lts s (Formula.modal a φ) := by
+  intro s' htr hφ
+  simp only [satisfies]
+  exists s'
+
+
+
+/-- A state satisfies a conjunction if and only if it satisfies all conjuncts. -/
+theorem satisfies_conj (lts : LTS State Label) (s : State) (φs : List (Formula Label)) :
+  satisfies lts s (Formula.conj φs) ↔ ∀ φ ∈ φs, satisfies lts s φ := by
+  simp only [satisfies]
+
+
+/-- Double negation elimination for HML. -/
+theorem satisfies_double_neg (lts : LTS State Label) (s : State) (φ : Formula Label) :
+  satisfies lts s (Formula.neg (Formula.neg φ)) ↔ satisfies lts s φ := by
+  simp only [satisfies, not_not]
+
+
+/-- Size of a formula (for well-founded recursion) -/
+def Formula.size : Formula Label → Nat
+  | .true => 1
+  | .false => 1
+  | .modal _ φ => 1 + φ.size
+  | .conj φs => 1 + φs.foldr (fun φ acc => φ.size + acc) 0
+  | .neg φ => 1 + φ.size
+
+
+
+/-- induction on formula structure -/
+theorem Formula.ind_on {Label : Type v} {P : Formula Label → Prop}
+  (φ : Formula Label)
+  (h_true : P .true)
+  (h_false : P .false)
+  (h_modal : ∀ a ψ, P ψ → P (.modal a ψ))
+  (h_conj : ∀ φs, (∀ φ ∈ φs, P φ) → P (.conj φs))
+  (h_neg : ∀ ψ, P ψ → P (.neg ψ))
+  : P φ := by
+  apply Formula.rec
+    (motive_1 := P)
+    (motive_2 := fun φs => ∀ φ ∈ φs, P φ)
+  · -- true case
+    exact h_true
+  · -- false case
+    exact h_false
+  · -- modal case
+    exact h_modal
+  · -- conj case
+    intro φs h_all
+    exact h_conj φs h_all
+  · -- neg case
+    exact h_neg
+  · -- nil case (empty list)
+    intros φ h_mem
+    cases h_mem
+  · -- cons case
+    intros head tail h_head h_tail φ h_mem
+    cases h_mem with
+    | head h_eq => exact h_head
+    | tail h_in =>
+      apply h_tail
+      trivial
+
+/-- If two LTS have the same transition relation, then they satisfy the same formulas -/
+theorem satisfies_independent_of_lts_structure
+  {State : Type u} {Label : Type v}
+  (lts1 lts2 : LTS State Label)
+  (h_same_tr : ∀ s a s', lts1.Tr s a s' ↔ lts2.Tr s a s')
+  (s : State) (φ : Formula Label) :
+  satisfies lts1 s φ ↔ satisfies lts2 s φ := by
+  induction φ using Formula.ind_on generalizing s with
+  | h_true =>
+    -- True case: both are True
+    simp [satisfies]
+  | h_false =>
+    -- False case: both are False
+    simp [satisfies]
+  | h_modal a ψ ih =>
+    -- Modal case: ⟨a⟩ψ
+    simp only [satisfies]
+    constructor
+    · intro ⟨s', ⟨h_tr1, h_sat1⟩⟩
+      use s'
+      constructor
+      · rw [← h_same_tr]
+        exact h_tr1
+      · rw [← ih s']
+        exact h_sat1
+    · intro ⟨s', ⟨h_tr2, h_sat2⟩⟩
+      use s'
+      constructor
+      · rw [h_same_tr]
+        exact h_tr2
+      · rw [ih s']
+        exact h_sat2
+  | h_conj φs ih_list =>
+    -- Conjunction case: ⋀ᵢ φᵢ
+    simp only [satisfies]
+    constructor
+    · intro h_all φ h_mem
+      rw [← ih_list φ h_mem s]
+      exact h_all φ h_mem
+    · intro h_all φ h_mem
+      rw [ih_list φ h_mem s]
+      exact h_all φ h_mem
+  | h_neg ψ ih =>
+    -- Negation case: ¬ψ
+    simp only [satisfies]
+    rw [ih s]
+
+end HennessyMilner
+
+end Cslib
diff --git a/CslibTests.lean b/CslibTests.lean
index 2e443461..e4ce07eb 100644
--- a/CslibTests.lean
+++ b/CslibTests.lean
@@ -5,6 +5,7 @@ public import CslibTests.CCS
 public import CslibTests.DFA
 public import CslibTests.FreeMonad
 public import CslibTests.GrindLint
+public import CslibTests.HMLGame
 public import CslibTests.HasFresh
 public import CslibTests.ImportWithMathlib
 public import CslibTests.LTS
diff --git a/CslibTests/HMLGame.lean b/CslibTests/HMLGame.lean
new file mode 100644
index 00000000..34329ef4
--- /dev/null
+++ b/CslibTests/HMLGame.lean
@@ -0,0 +1,157 @@
+import Cslib.Foundations.Semantics.GameSemantics.Basic
+import Cslib.Foundations.Semantics.GameSemantics.HMLGame
+
+namespace CslibTests
+
+open Cslib
+open Cslib.HennessyMilner
+open Cslib.HMLGame
+
+/-!
+## Example 2.9 from Bisping et al.
+
+Testing that ⟦⟨a⟩¬⟨d⟩⊤⟧^CCS_{P₂} is false.
+-/
+
+/-- Actions for the CCS processes -/
+inductive Action where
+  | a : Action
+  | b : Action
+  | c : Action
+  | d : Action
+  deriving DecidableEq, Repr
+
+/-- States for process P₁ -/
+inductive StateP1 where
+  | P1 : StateP1       -- Initial state
+  | bc : StateP1       -- b + c
+  | d : StateP1        -- d
+  | zero : StateP1     -- 0 (terminated)
+  deriving DecidableEq, Repr
+
+/-- States for process P₂ -/
+inductive StateP2 where
+  | P2 : StateP2       -- Initial state
+  | bd : StateP2       -- b + d
+  | cd : StateP2       -- c + d
+  | zero : StateP2     -- 0 (terminated)
+  deriving DecidableEq, Repr
+
+/-- Transition relation for P₁ -/
+inductive TrP1 : StateP1 → Action → StateP1 → Prop where
+  | p1_a_bc : TrP1 .P1 .a .bc
+  | p1_a_d  : TrP1 .P1 .a .d
+  | bc_b_0  : TrP1 .bc .b .zero
+  | bc_c_0  : TrP1 .bc .c .zero
+  | d_d_0   : TrP1 .d  .d .zero
+
+/-- Transition relation for P₂ -/
+inductive TrP2 : StateP2 → Action → StateP2 → Prop where
+  | p2_a_bd : TrP2 .P2 .a .bd
+  | p2_a_cd : TrP2 .P2 .a .cd
+  | bd_b_0  : TrP2 .bd .b .zero
+  | bd_d_0  : TrP2 .bd .d .zero
+  | cd_c_0  : TrP2 .cd .c .zero
+  | cd_d_0  : TrP2 .cd .d .zero
+
+/-- LTS for process P₁ -/
+def ltsP1 : LTS StateP1 Action := ⟨TrP1⟩
+
+/-- LTS for process P₂ -/
+def ltsP2 : LTS StateP2 Action := ⟨TrP2⟩
+
+/-!
+### Formula: ⟨a⟩¬⟨d⟩⊤
+
+This formula says: "there exists an a-transition to a state where
+there is no d-transition to any state satisfying ⊤"
+
+Since ⊤ is always true, ¬⟨d⟩⊤ means "there is no d-transition".
+-/
+
+/-- The formula ⊤ (empty conjunction) -/
+def formulaTrue : Formula Action := .conj []
+
+/-- The formula ⟨d⟩⊤ -/
+def formulaDiamondD : Formula Action := .modal .d formulaTrue
+
+/-- The formula ¬⟨d⟩⊤ -/
+def formulaNegDiamondD : Formula Action := .neg formulaDiamondD
+
+/-- The formula ⟨a⟩¬⟨d⟩⊤ -/
+def formulaExample : Formula Action := .modal .a formulaNegDiamondD
+
+/-!
+### Example 2.9: ⟦⟨a⟩¬⟨d⟩⊤⟧^CCS_{P₂} is false
+
+From P₂, the defender can move via 'a' to either (b+d) or (c+d).
+In both cases, there IS a d-transition to 0, so ¬⟨d⟩⊤ is false.
+Therefore ⟨a⟩¬⟨d⟩⊤ is false at P₂.
+-/
+
+/-- Helper: 0 satisfies ⊤ -/
+theorem zero_satisfies_true : satisfies ltsP2 .zero formulaTrue := by
+  simp [formulaTrue, satisfies]
+
+/-- Helper: b+d has a d-transition -/
+theorem bd_has_d_transition : ltsP2.Tr .bd .d .zero := TrP2.bd_d_0
+
+/-- Helper: c+d has a d-transition -/
+theorem cd_has_d_transition : ltsP2.Tr .cd .d .zero := TrP2.cd_d_0
+
+/-- Helper: b+d satisfies ⟨d⟩⊤ -/
+theorem bd_satisfies_diamondD : satisfies ltsP2 .bd formulaDiamondD := by
+  simp only [formulaDiamondD, formulaTrue, satisfies]
+  simp only [List.not_mem_nil, IsEmpty.forall_iff, implies_true, and_true]
+  exists .zero
+  constructor
+
+
+/-- Helper: c+d satisfies ⟨d⟩⊤ -/
+theorem cd_satisfies_diamondD : satisfies ltsP2 .cd formulaDiamondD := by
+  simp only [formulaDiamondD, formulaTrue, satisfies, List.not_mem_nil, IsEmpty.forall_iff,
+    implies_true, and_true]
+  exact ⟨.zero, by constructor⟩
+
+/-- Helper: b+d does NOT satisfy ¬⟨d⟩⊤ -/
+theorem bd_not_satisfies_negDiamondD : ¬ satisfies ltsP2 .bd formulaNegDiamondD := by
+  simp only [formulaNegDiamondD, satisfies, not_not]
+  exact bd_satisfies_diamondD
+
+/-- Helper: c+d does NOT satisfy ¬⟨d⟩⊤ -/
+theorem cd_not_satisfies_negDiamondD : ¬ satisfies ltsP2 .cd formulaNegDiamondD := by
+  simp only [formulaNegDiamondD, satisfies, not_not]
+  exact cd_satisfies_diamondD
+
+/--
+Example 2.9: P₂ does NOT satisfy ⟨a⟩¬⟨d⟩⊤
+
+The defender cannot win because both a-successors (b+d and c+d)
+have d-transitions, so ¬⟨d⟩⊤ fails at both.
+-/
+theorem example_2_9 : ¬ satisfies ltsP2 .P2 formulaExample := by
+  simp only [formulaExample, formulaNegDiamondD, formulaDiamondD, formulaTrue, satisfies]
+  intro ⟨p', h_tr, h_neg⟩
+  -- p' must be either bd or cd (the only a-successors of P2)
+  cases h_tr with
+  | p2_a_bd =>
+    -- p' = bd, but bd has a d-transition to zero
+    apply h_neg
+    simp only [List.not_mem_nil, IsEmpty.forall_iff, implies_true, and_true]
+    exact ⟨.zero, by constructor⟩
+  | p2_a_cd =>
+    -- p' = cd, but cd has a d-transition to zero
+    apply h_neg
+    simp only [List.not_mem_nil, IsEmpty.forall_iff, implies_true, and_true]
+    exact ⟨.zero,  by constructor⟩
+
+/--
+Using the game semantics theorem to show the same result:
+The defender does NOT win G^S_HML[(P₂, ⟨a⟩¬⟨d⟩⊤)]
+-/
+theorem example_2_9_game : ¬ defenderWinsHML ltsP2 .P2 formulaExample := by
+  rw [game_semantics_soundness]
+  exact example_2_9
+
+
+end CslibTests
diff --git a/references.bib b/references.bib
index 551d6819..ae0f32aa 100644
--- a/references.bib
+++ b/references.bib
@@ -195,3 +195,13 @@ @incollection{ Thomas1990
   publisher    = {Elsevier and {MIT} Press},
   year         = {1990}
 }
+
+@book{AstromMurray2020,
+  author    = {Åström, Karl J. and Murray, Richard M.},
+  title     = {Feedback Systems: An Introduction for Scientists and Engineers},
+  edition   = {2nd},
+  publisher = {Princeton University Press},
+  year      = {2020},
+  url       = {https://www.cds.caltech.edu/~murray/books/AM08/pdf/fbs-public_24Jul2020.pdf},
+  note      = {Electronic public draft, Version v3.1.5 (24 Jul 2020)}
+}
diff --git a/test_error.lean b/test_error.lean
new file mode 100644
index 00000000..d1a1f6c7
--- /dev/null
+++ b/test_error.lean
@@ -0,0 +1,15 @@
+import Cslib.Foundations.Semantics.LTS.Basic
+import Cslib.Logics.HennessyMilnerLogic.Basic
+
+namespace Cslib
+open HennessyMilner
+
+universe u v
+
+structure HMLGame where
+  Label : Type v
+  Process : Type u
+  S : LTS Process Label
+  G : Process × Formula Label
+  -- This is the problematic line adapted
+  G_d : {(p, φ) : Process × Formula Label | φ = Formula.modal _ _ }