diff --git a/src/SUNRepresentations.jl b/src/SUNRepresentations.jl
index 0f2bcec..c659d74 100644
--- a/src/SUNRepresentations.jl
+++ b/src/SUNRepresentations.jl
@@ -13,6 +13,7 @@ export SUNIrrep, basis, weight, Zweight, creation, annihilation, highest_weight,
 export directproduct, CGC
 export SU, SU₃, SU₄, SU₅, SU3Irrep, SU4Irrep, SU5Irrep
 export dynkin_label, congruency
+export casimir
 
 """
     struct SUNIrrep{N} <: AbstractIrrep{SU{N}}
@@ -101,6 +102,34 @@ end
 
 include("gtpatterns.jl")
 
+"""
+    casimir(k::Int, irrep::SUNIrrep)
+
+Return the eigenvalue of the `k`-th order Casimir operator in the representation `irrep`.
+
+The formula is:
+```math
+C_k(\\lambda) = \\frac{1}{2} \\left[ \\sum_{i=1}^{N} L_i^k - \\sum_{i=1}^{N} \\rho_i^k \\right]
+```
+where ``\\rho_i = (N+1-2i)/2`` are the Weyl vector components and
+``L_i = (\\lambda_i - \\bar\\lambda) + \\rho_i`` are the shifted traceless weights
+(``\\bar\\lambda = \\sum_j \\lambda_j / N``).
+
+The independent primitive Casimir operators have orders ``k = 2, 3, \\ldots, N``.
+Other values of `k` are valid but give dependent (or zero) results.
+"""
+function casimir(k::Int, irrep::SUNIrrep{N}) where {N}
+    λ = weight(irrep)
+    λ̄ = sum(λ) // N
+    c = zero(Rational{Int})
+    for i in 1:N
+        ρᵢ = (N + 1 - 2i) // 2
+        Lᵢ = (λ[i] - λ̄) + ρᵢ
+        c += Lᵢ^k - ρᵢ^k
+    end
+    return c / 2
+end
+
 basis(s::SUNIrrep) = GTPatternIterator(s)
 
 # direct product: return dictionary with new irreps as keys, outer multiplicities as value
diff --git a/test/casimir.jl b/test/casimir.jl
new file mode 100644
index 0000000..47ceb4c
--- /dev/null
+++ b/test/casimir.jl
@@ -0,0 +1,29 @@
+@testset "casimir" begin
+    # --- SU(2): C₂(m) = m(m+2)/4, all odd C_k = 0 ---
+    for m in 0:10
+        @test casimir(2, SUNIrrep{2}([m])) == m * (m + 2) // 4
+        @test casimir(3, SUNIrrep{2}([m])) == 0
+    end
+
+    # --- SU(3): C₂(p,q) = (p²+q²+pq+3p+3q)/3, C₃(p,q) = (p-q)(2p+q+3)(p+2q+3)/18 ---
+    for p in 0:10, q in 0:10
+        irrep = SUNIrrep{3}([p, q])
+        @test casimir(2, irrep) == (p^2 + q^2 + p * q + 3p + 3q) // 3
+        @test casimir(3, irrep) == (p - q) * (2p + q + 3) * (p + 2q + 3) // 18
+    end
+
+    # --- k > 3: structural properties only ---
+    # Singlet always 0
+    for k in 2:5, N in 2:5
+        @test casimir(k, one(SUNIrrep{N})) == 0 // 1
+    end
+
+    # Conjugation: C_k(dual(λ)) = (-1)^k * C_k(λ)
+    for N in 2:5
+        for irrep in Iterators.take(values(SUNIrrep{N}), 20)
+            for k in 2:5
+                @test casimir(k, dual(irrep)) == (-1)^k * casimir(k, irrep)
+            end
+        end
+    end
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index d46695a..f178d8e 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -46,6 +46,7 @@ module GenericTests
 end
 
 include("caching.jl")
+include("casimir.jl")
 sectorlist = (SUNIrrep{3}, SUNIrrep{4}, SUNIrrep{5}, SUNIrrep{3} ⊠ SUNIrrep{3})
 include("sectors.jl")
 sectorlist = (SUNIrrep{3}, SUNIrrep{4}, SUNIrrep{5})