diff --git a/src/hubbardoperators.jl b/src/hubbardoperators.jl
index 806b179..343a8f8 100644
--- a/src/hubbardoperators.jl
+++ b/src/hubbardoperators.jl
@@ -245,7 +245,7 @@ function u_plus_u_min(T, ::Type{Trivial}, ::Type{U1Irrep})
     return t
 end
 function u_plus_u_min(T, ::Type{Trivial}, ::Type{SU2Irrep})
-    return error("Not implemented")
+    throw(ArgumentError("`u_plus_u_min` is not symmetric under `SU2Irrep` spin symmetry"))
 end
 function u_plus_u_min(T, ::Type{U1Irrep}, ::Type{Trivial})
     t = two_site_operator(T, U1Irrep, Trivial)
@@ -266,7 +266,7 @@ function u_plus_u_min(T, ::Type{U1Irrep}, ::Type{U1Irrep})
     return t
 end
 function u_plus_u_min(T, ::Type{U1Irrep}, ::Type{SU2Irrep})
-    return error("Not implemented")
+    throw(ArgumentError("`u_plus_u_min` is not symmetric under `SU2Irrep` spin symmetry"))
 end
 function u_plus_u_min(T, ::Type{SU2Irrep}, ::Type{Trivial})
     return error("Not implemented")
@@ -275,7 +275,7 @@ function u_plus_u_min(T, ::Type{SU2Irrep}, ::Type{U1Irrep})
     return error("Not implemented")
 end
 function u_plus_u_min(T, ::Type{SU2Irrep}, ::Type{SU2Irrep})
-    return error("Not implemented")
+    throw(ArgumentError("`u_plus_u_min` is not symmetric under `SU2Irrep` spin symmetry"))
 end
 const u⁺u⁻ = u_plus_u_min
 
@@ -305,7 +305,7 @@ function d_plus_d_min(T, ::Type{Trivial}, ::Type{U1Irrep})
     return t
 end
 function d_plus_d_min(T, ::Type{Trivial}, ::Type{SU2Irrep})
-    return error("Not implemented")
+    throw(ArgumentError("`d_plus_d_min` is not symmetric under `SU2Irrep` spin symmetry"))
 end
 function d_plus_d_min(T, ::Type{U1Irrep}, ::Type{Trivial})
     t = two_site_operator(T, U1Irrep, Trivial)
@@ -326,7 +326,7 @@ function d_plus_d_min(T, ::Type{U1Irrep}, ::Type{U1Irrep})
     return t
 end
 function d_plus_d_min(T, ::Type{U1Irrep}, ::Type{SU2Irrep})
-    return error("Not implemented")
+    throw(ArgumentError("`d_plus_d_min` is not symmetric under `SU2Irrep` spin symmetry"))
 end
 function d_plus_d_min(T, ::Type{SU2Irrep}, ::Type{Trivial})
     return error("Not implemented")
@@ -335,7 +335,7 @@ function d_plus_d_min(T, ::Type{SU2Irrep}, ::Type{U1Irrep})
     return error("Not implemented")
 end
 function d_plus_d_min(T, ::Type{SU2Irrep}, ::Type{SU2Irrep})
-    return error("Not implemented")
+    throw(ArgumentError("`d_plus_d_min` is not symmetric under `SU2Irrep` spin symmetry"))
 end
 const d⁺d⁻ = d_plus_d_min
 
diff --git a/test/bosonoperators.jl b/test/bosonoperators.jl
index 7ab75d1..cf8c92c 100644
--- a/test/bosonoperators.jl
+++ b/test/bosonoperators.jl
@@ -6,9 +6,9 @@ using .TensorKitTensorsTestSetup
 using TensorKitTensors.BosonOperators
 using StableRNGs
 
-const cutoff = 4
-
 @testset "Non-symmetric bosonic operators" begin
+    cutoff = 4
+
     # inferrability
     A = @inferred a(; cutoff)
     A⁺ = @inferred a⁺(; cutoff)
@@ -42,6 +42,8 @@ const cutoff = 4
 end
 
 @testset "U1-symmetric bosonic operators" begin
+    cutoff = 4
+
     rng = StableRNG(123)
     # inferrability
     N = @inferred n(U1Irrep; cutoff)
@@ -66,3 +68,24 @@ end
 
     test_operator(O_u1, O_triv; L)
 end
+
+@testset "Exact Diagonalization" begin
+    cutoff = 1
+    L = 2
+    for symmetry in [Trivial U1Irrep]
+        rng = StableRNG(123)
+        # inferrability
+        N = @inferred n(U1Irrep; cutoff)
+        A⁺A = @inferred a⁺a(U1Irrep; cutoff)
+        AA⁺ = @inferred aa⁺(U1Irrep; cutoff)
+        V = @inferred boson_space(U1Irrep; cutoff)
+
+        a_pm, a_mp, a_n = rand(rng, 3)
+        O = (N ⊗ id(V) + id(V) ⊗ N) * a_n + AA⁺ * a_mp + A⁺A * a_pm
+
+        true_eigenvals = sort([0, 2 * a_n, a_n + sqrt(a_mp * a_pm),
+                               a_n - sqrt(a_mp * a_pm)])
+        eigenvals = expanded_eigenvalues(O; L)
+        @test eigenvals ≈ true_eigenvals
+    end
+end
diff --git a/test/fermionoperators.jl b/test/fermionoperators.jl
index 13d1af7..1c63a38 100644
--- a/test/fermionoperators.jl
+++ b/test/fermionoperators.jl
@@ -37,9 +37,9 @@ end
     H = -t * (c⁻c⁺() + c⁺c⁻()) +
         V * ((n() - 0.5 * id(pspace)) ⊗ (n() - 0.5 * id(pspace))) -
         0.5 * mu * (n() ⊗ id(pspace) + id(pspace) ⊗ n())
-    # Values based on https://arxiv.org/abs/1610.05003v1. Half-Chain Entanglement Entropy in the One-Dimensional Spinless Fermion Model
-    true_eigenvalues = sort([V / 4, V / 4 - mu, -V / 4 - mu / 2 + t, -V / 4 - mu / 2 - t])
 
+    # Values based on https://arxiv.org/abs/1610.05003v1. Half-Chain Entanglement Entropy in the One-Dimensional Spinless Fermion Model
+    true_eigenvals = sort([V / 4, V / 4 - mu, -V / 4 - mu / 2 + t, -V / 4 - mu / 2 - t])
     eigenvals = expanded_eigenvalues(H; L)
-    @test eigenvals ≈ true_eigenvalues
+    @test eigenvals ≈ true_eigenvals
 end
diff --git a/test/hubbardoperators.jl b/test/hubbardoperators.jl
index 4905663..f8944c7 100644
--- a/test/hubbardoperators.jl
+++ b/test/hubbardoperators.jl
@@ -64,10 +64,15 @@ end
                       d_num(particle_symmetry, spin_symmetry) ≈
                       d_num(particle_symmetry, spin_symmetry) *
                       u_num(particle_symmetry, spin_symmetry)
+            else
+                @test_throws ArgumentError u_plus_u_min(particle_symmetry, spin_symmetry)
+                @test_throws ArgumentError d_plus_d_min(particle_symmetry, spin_symmetry)
             end
         else
             @test_broken c_plus_c_min(particle_symmetry, spin_symmetry)
             @test_broken c_min_c_plus(particle_symmetry, spin_symmetry)
+            @test_broken d_plus_d_min(particle_symmetry, spin_symmetry)
+            @test_broken u_plus_u_min(particle_symmetry, spin_symmetry)
         end
     end
 end
@@ -114,3 +119,28 @@ end
         @test vals_triv ≈ vals_symm
     end
 end
+
+@testset "Exact diagonalisation" begin
+    for particle_symmetry in [Trivial, U1Irrep, SU2Irrep],
+        spin_symmetry in [Trivial, U1Irrep, SU2Irrep]
+
+        if (particle_symmetry, spin_symmetry) in implemented_symmetries
+            rng = StableRNG(123)
+
+            L = 2
+            t, U = rand(rng, 5)
+            mu = 0.0
+            E⁻ = U / 2 - sqrt((U / 2)^2 + 4 * t^2)
+            E⁺ = U / 2 + sqrt((U / 2)^2 + 4 * t^2)
+            H_triv = hubbard_hamiltonian(particle_symmetry, spin_symmetry; t, U, mu, L)
+
+            # Values based on https://arxiv.org/pdf/0807.4878. Introduction to Hubbard Model and Exact Diagonalization
+            true_eigenvals = sort(vcat(repeat([-t], 2), [E⁻], repeat([0], 4),
+                                       repeat([t], 2),
+                                       repeat([U - t], 2), [U], [E⁺], repeat([U + t], 2),
+                                       [2 * U]))
+            eigenvals = expanded_eigenvalues(H_triv; L)
+            @test eigenvals ≈ true_eigenvals
+        end
+    end
+end
diff --git a/test/spinoperators.jl b/test/spinoperators.jl
index 8feb6fd..05b9e02 100644
--- a/test/spinoperators.jl
+++ b/test/spinoperators.jl
@@ -125,3 +125,20 @@ end
              S_exchange(; spin) * a_exchange
     test_operator(O_u1, O_triv; L)
 end
+
+@testset "Exact diagonalisation for $sector symmetry" for sector in [Trivial U1Irrep]
+    spin = 1
+
+    ZZ = @inferred S_zz(sector; spin)
+    plusmin = @inferred S_plusmin(sector; spin)
+    minplus = @inferred S_minplus(sector; spin)
+    O = ZZ + 0.5 * (plusmin + minplus)
+
+    true_eigenvals = vcat([-2.0], repeat([-1.0], 3), repeat([1.0], 5))
+    eigenvals = expanded_eigenvalues(O; L=2)
+    @test eigenvals ≈ true_eigenvals
+
+    # Value based on https://doi.org/10.1088/0953-8984/2/26/010. Exact diagonalisations of open spin-1 chains
+    eigenvals = expanded_eigenvalues(O; L=4)
+    @test eigenvals[2] - eigenvals[1] ≈ 0.509170 atol = 1e-6
+end
diff --git a/test/tjoperators.jl b/test/tjoperators.jl
index 32e8c97..378a9ba 100644
--- a/test/tjoperators.jl
+++ b/test/tjoperators.jl
@@ -53,14 +53,18 @@ end
                     @test u_plus_u_min(particle_symmetry, spin_symmetry; slave_fermion)' ≈
                           u_min_u_plus(particle_symmetry, spin_symmetry; slave_fermion)
                 else
-                    @test_broken d_plus_d_min(particle_symmetry, spin_symmetry;
-                                              slave_fermion)
-                    @test_broken d_min_d_plus(particle_symmetry, spin_symmetry;
-                                              slave_fermion)
-                    @test_broken u_plus_u_min(particle_symmetry, spin_symmetry;
-                                              slave_fermion)
-                    @test_broken u_min_u_plus(particle_symmetry, spin_symmetry;
-                                              slave_fermion)
+                    @test_throws ArgumentError d_plus_d_min(particle_symmetry,
+                                                            spin_symmetry;
+                                                            slave_fermion)
+                    @test_throws ArgumentError d_min_d_plus(particle_symmetry,
+                                                            spin_symmetry;
+                                                            slave_fermion)
+                    @test_throws ArgumentError u_plus_u_min(particle_symmetry,
+                                                            spin_symmetry;
+                                                            slave_fermion)
+                    @test_throws ArgumentError u_min_u_plus(particle_symmetry,
+                                                            spin_symmetry;
+                                                            slave_fermion)
                 end
 
                 # test number operator
@@ -77,8 +81,10 @@ end
                           c_num_hole(particle_symmetry, spin_symmetry; slave_fermion) +
                           c_num(particle_symmetry, spin_symmetry; slave_fermion)
                 else
-                    @test_broken u_num(particle_symmetry, spin_symmetry; slave_fermion)
-                    @test_broken d_num(particle_symmetry, spin_symmetry; slave_fermion)
+                    @test_throws ArgumentError u_num(particle_symmetry, spin_symmetry;
+                                                     slave_fermion)
+                    @test_throws ArgumentError d_num(particle_symmetry, spin_symmetry;
+                                                     slave_fermion)
                 end
 
                 # test spin operator
@@ -88,11 +94,12 @@ end
                            d_min_u_min(particle_symmetry, spin_symmetry; slave_fermion)) /
                           sqrt(2)
                 else
-                    @test_broken c_singlet(particle_symmetry, spin_symmetry; slave_fermion)
-                    @test_broken u_min_d_min(particle_symmetry, spin_symmetry;
-                                             slave_fermion)
-                    @test_broken d_min_u_min(particle_symmetry, spin_symmetry;
-                                             slave_fermion)
+                    @test_throws ArgumentError c_singlet(particle_symmetry, spin_symmetry;
+                                                         slave_fermion)
+                    @test_throws ArgumentError u_min_d_min(particle_symmetry, spin_symmetry;
+                                                           slave_fermion)
+                    @test_throws ArgumentError d_min_u_min(particle_symmetry, spin_symmetry;
+                                                           slave_fermion)
                 end
 
                 if spin_symmetry == Trivial
@@ -166,12 +173,11 @@ function tjhamiltonian(particle_symmetry, spin_symmetry; t, J, mu, L, slave_ferm
 end
 
 @testset "spectrum" begin
+    rng = StableRNG(123)
     L = 4
-    t = randn()
-    J = randn()
-    mu = randn()
 
     for slave_fermion in (false, true)
+        t, J, mu = rand(rng, 3)
         H_triv = tjhamiltonian(Trivial, Trivial; t, J, mu, L, slave_fermion)
         vals_triv = mapreduce(vcat, eigvals(H_triv)) do (c, v)
             return repeat(real.(v), dim(c))
@@ -199,3 +205,30 @@ end
         end
     end
 end
+
+@testset "Exact Diagonalisation" begin
+    rng = StableRNG(123)
+    L = 2
+
+    for particle_symmetry in [Trivial, U1Irrep],
+        spin_symmetry in [Trivial, U1Irrep, SU2Irrep]
+
+        for slave_fermion in (false, true)
+            if (particle_symmetry, spin_symmetry) in implemented_symmetries
+                t, J = rand(rng, 2)
+                num = c_num(particle_symmetry, spin_symmetry; slave_fermion)
+                H = (-t) *
+                    (c_plus_c_min(particle_symmetry, spin_symmetry; slave_fermion) +
+                     c_min_c_plus(particle_symmetry, spin_symmetry; slave_fermion)) +
+                    J *
+                    (S_exchange(particle_symmetry, spin_symmetry; slave_fermion) -
+                     (1 / 4) * (num ⊗ num))
+
+                true_eigenvals = sort(vcat([-J], repeat([-t], 2), repeat([t], 2),
+                                           repeat([0.0], 4)))
+                eigenvals = expanded_eigenvalues(H; L)
+                @test eigenvals ≈ true_eigenvals
+            end
+        end
+    end
+end