diff --git a/ext/QuantumCliffordQOpticsExt/QuantumCliffordQOpticsExt.jl b/ext/QuantumCliffordQOpticsExt/QuantumCliffordQOpticsExt.jl
index cca3011df..da7fe3c5d 100644
--- a/ext/QuantumCliffordQOpticsExt/QuantumCliffordQOpticsExt.jl
+++ b/ext/QuantumCliffordQOpticsExt/QuantumCliffordQOpticsExt.jl
@@ -6,6 +6,7 @@ using QuantumOpticsBase
 using Graphs
 using DocStringExtensions
 import QuantumOpticsBase: Ket, Operator
+import QuantumClifford: cauchy_schwarz_check
 
 const _b2 = SpinBasis(1//2)
 const _l0 = spinup(_b2)
@@ -234,6 +235,24 @@ function Operator(c::CliffordOperator)
     cliff_to_unitary(c)
 end
 
+"""
+$TYPEDSIGNATURES
+
+Determines whether two [`GeneralizedStabilizer`](@ref) states, `sm₁` and `sm₂`, satisfy the
+*Cauchy-Schwarz inequality* within the the Hilbert space for `n`-qubit density matrices.
+
+The *Cauchy-Schwarz inequality* for `n`-qubit density matrices is given by:
+
+```math
+\\left| \\text{Tr}[sm_1' sm_2] \\right| \\leq \\sqrt{\\text{Tr}[(sm_1')^2] \\, \\text{Tr}[sm_2^2]}
+```
+
+Equality holds if and only if: `sm₁ ==  sm₂`
+"""
+function cauchy_schwarz_check(sm₁::GeneralizedStabilizer, sm₂::GeneralizedStabilizer)
+    return abs(real(tr(Operator(sm₁)' * Operator(sm₂)))) <= sqrt(real(tr((Operator(sm₁)')^2) * tr(Operator(sm₂)^2)))
+end
+
 end
 
 ##
diff --git a/src/QuantumClifford.jl b/src/QuantumClifford.jl
index 6732cc6f8..8ef574ec3 100644
--- a/src/QuantumClifford.jl
+++ b/src/QuantumClifford.jl
@@ -81,7 +81,7 @@ export
     # petrajectories
     petrajectories, applybranches,
     # nonclifford
-    GeneralizedStabilizer, UnitaryPauliChannel, PauliChannel, pcT,
+    GeneralizedStabilizer, UnitaryPauliChannel, PauliChannel, pcT, cauchy_schwarz_check,
     # makie plotting -- defined only when extension is loaded
     stabilizerplot, stabilizerplot_axis,
     # sum types
diff --git a/src/nonclifford.jl b/src/nonclifford.jl
index 524085025..4f9ebf4c4 100644
--- a/src/nonclifford.jl
+++ b/src/nonclifford.jl
@@ -462,3 +462,15 @@ const pcT = UnitaryPauliChannel(
     (I, Z),
     ((1+exp(im*π/4))/2, (1-exp(im*π/4))/2)
 )
+
+##
+# QuantumOpticsBaseExt methods
+##
+
+function cauchy_schwarz_check(args...)
+    ext = Base.get_extension(QuantumClifford, :QuantumCliffordOpticsBase)
+    if isnothing(ext)
+        throw("The `cauchy_schwarz_check` depends on the package `QuantumOpticsBase` but you have not installed or imported it yet. Immediately after you import `QuantumOpticsBase`, the `cauchy_schwarz_check` will be available.")
+    end
+    return ext.cauchy_schwarz_check(args...)
+end
diff --git a/test/test_nonclifford_quantumoptics.jl b/test/test_nonclifford_quantumoptics.jl
index 3b4ba307f..1c5f2f394 100644
--- a/test/test_nonclifford_quantumoptics.jl
+++ b/test/test_nonclifford_quantumoptics.jl
@@ -1,5 +1,5 @@
 using QuantumClifford
-using QuantumClifford: GeneralizedStabilizer, rowdecompose, PauliChannel, mul_left!, mul_right!
+using QuantumClifford: GeneralizedStabilizer, rowdecompose, PauliChannel, mul_left!, mul_right!, cauchy_schwarz_check
 using QuantumClifford: @S_str, random_stabilizer
 using QuantumOpticsBase
 using LinearAlgebra
@@ -82,6 +82,7 @@ end
         apply!(sm, embed(n, i, pcT))
         smcopy = copy(sm)
         @test smcopy == sm
+        @test cauchy_schwarz_check(copy(smcopy), sm)
         nc = embed(n, rand(1:n), pcT)
         @test copy(nc) == nc
     end