diff --git a/src/hubbardoperators.jl b/src/hubbardoperators.jl
index ad20f1e..3954b9d 100644
--- a/src/hubbardoperators.jl
+++ b/src/hubbardoperators.jl
@@ -3,7 +3,7 @@ module HubbardOperators
 using TensorKit
 
 export hubbard_space
-export e_num, u_num, d_num, ud_num
+export e_num, u_num, d_num, ud_num, half_ud_num
 export S_x, S_y, S_z, S_plus, S_min
 export u_plus_u_min, d_plus_d_min
 export u_min_u_plus, d_min_d_plus
@@ -56,15 +56,17 @@ function hubbard_space(::Type{U1Irrep}, ::Type{SU2Irrep})
     )
 end
 function hubbard_space(::Type{SU2Irrep}, ::Type{Trivial})
-    return Vect[FermionParity ⊠ SU2Irrep]((0, 0) => 2, (1, 1 // 2) => 1)
+    return Vect[FermionParity ⊠ SU2Irrep]((0, 1 // 2) => 1, (1, 0) => 2)
 end
 function hubbard_space(::Type{SU2Irrep}, ::Type{U1Irrep})
     return Vect[FermionParity ⊠ SU2Irrep ⊠ U1Irrep](
-        (0, 0, 0) => 1, (1, 1 // 2, 1) => 1, (0, 0, 2) => 1
+        (0, 1 // 2, 0) => 1, (1, 0, -1 // 2) => 1, (1, 0, 1 // 2) => 1
     )
 end
 function hubbard_space(::Type{SU2Irrep}, ::Type{SU2Irrep})
-    return Vect[FermionParity ⊠ SU2Irrep ⊠ SU2Irrep]((1, 1 // 2, 1 // 2) => 1)
+    return Vect[FermionParity ⊠ SU2Irrep ⊠ SU2Irrep](
+        (0, 1 // 2, 0) => 1, (1, 0, 1 // 2) => 1
+    )
 end
 
 # Single-site operators
@@ -239,6 +241,27 @@ function ud_num(elt::Type{<:Number}, ::Type{U1Irrep}, ::Type{SU2Irrep})
 end
 const nꜛꜜ = ud_num
 
+@doc """
+    half_ud_num([elt::Type{<:Number}], [particle_symmetry::Type{<:Sector}], [spin_symmetry::Type{<:Sector}])
+
+Return the the one-body operator that is equivalent to `(nꜛ - 1/2)(nꜜ - 1/2)`, which respects the particle-hole symmetry.
+"""
+half_ud_num(P::Type{<:Sector}, S::Type{<:Sector}) = half_ud_num(ComplexF64, P, S)
+function half_ud_num(
+        elt::Type{<:Number}, particle_symmetry::Type{<:Sector},
+        spin_symmetry::Type{<:Sector}
+    )
+    I = id(hubbard_space(particle_symmetry, spin_symmetry))
+    return (u_num(elt, particle_symmetry, spin_symmetry) - I / 2) *
+        (d_num(elt, particle_symmetry, spin_symmetry) - I / 2)
+end
+function half_ud_num(elt::Type{<:Number}, ::Type{SU2Irrep}, ::Type{SU2Irrep})
+    t = single_site_operator(elt, SU2Irrep, SU2Irrep)
+    block(t, sectortype(t)(0, 1 // 2, 0)) .= 1 // 4
+    block(t, sectortype(t)(1, 0, 1 // 2)) .= -1 // 4
+    return t
+end
+
 @doc """
     S_plus(elt::Type{<:Number}, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector})
     S⁺(elt::Type{<:Number}, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector})
@@ -586,6 +609,22 @@ function e_hopping(
     return e_plus_e_min(elt, particle_symmetry, spin_symmetry) -
         e_min_e_plus(elt, particle_symmetry, spin_symmetry)
 end
+function e_hopping(elt::Type{<:Number}, ::Type{SU2Irrep}, ::Type{SU2Irrep})
+    elt <: Complex || throw(DomainError(elt, "SU₂ × SU₂ symmetry requires complex entries"))
+    t = two_site_operator(elt, SU2Irrep, SU2Irrep)
+    I = sectortype(t)
+    even = I(0, 1 // 2, 0)
+    odd = I(1, 0, 1 // 2)
+    f1 = only(fusiontrees((odd, odd), one(I)))
+    f2 = only(fusiontrees((even, even), one(I)))
+    t[f1, f2] .= 2im
+    t[f2, f1] .= -2im
+    f3 = only(fusiontrees((even, odd), I((1, 1 // 2, 1 // 2))))
+    f4 = only(fusiontrees((odd, even), I((1, 1 // 2, 1 // 2))))
+    t[f3, f4] .= im
+    t[f4, f3] .= -im
+    return t
+end
 const e_hop = e_hopping
 
 @doc """
diff --git a/test/hubbardoperators.jl b/test/hubbardoperators.jl
index 8d214a5..fec0d55 100644
--- a/test/hubbardoperators.jl
+++ b/test/hubbardoperators.jl
@@ -9,6 +9,7 @@ using StableRNGs
 implemented_symmetries = [
     (Trivial, Trivial), (Trivial, U1Irrep), (Trivial, SU2Irrep),
     (U1Irrep, Trivial), (U1Irrep, U1Irrep), (U1Irrep, SU2Irrep),
+    (SU2Irrep, SU2Irrep),
 ]
 
 @testset "Compare symmetric with trivial tensors" begin
@@ -16,7 +17,18 @@ implemented_symmetries = [
             spin_symmetry in [Trivial, U1Irrep, SU2Irrep]
 
         if (particle_symmetry, spin_symmetry) in implemented_symmetries
-            space = hubbard_space(particle_symmetry, spin_symmetry)
+            space = @inferred hubbard_space(particle_symmetry, spin_symmetry)
+
+            if particle_symmetry == spin_symmetry == SU2Irrep
+                O = e_hopping(ComplexF64, SU2Irrep, SU2Irrep)
+                O_triv = e_hopping(ComplexF64, Trivial, Trivial)
+                test_operator(O, O_triv)
+
+                O = half_ud_num(ComplexF64, SU2Irrep, SU2Irrep)
+                O_triv = half_ud_num(ComplexF64, Trivial, Trivial)
+                test_operator(O, O_triv)
+                continue
+            end
 
             O = e_plus_e_min(ComplexF64, particle_symmetry, spin_symmetry)
             O_triv = e_plus_e_min(ComplexF64, Trivial, Trivial)
@@ -50,6 +62,7 @@ end
         @test dim(space) == 4
 
         if (particle_symmetry, spin_symmetry) in implemented_symmetries
+            particle_symmetry == spin_symmetry == SU2Irrep && continue
             # test hopping operator
             epem = e_plus_e_min(particle_symmetry, spin_symmetry)
             emep = e_min_e_plus(particle_symmetry, spin_symmetry)
@@ -168,46 +181,77 @@ end
     end
 end
 
-function hubbard_hamiltonian(particle_symmetry, spin_symmetry; t, U, mu, L)
-    hopping = -t * e_hopping(particle_symmetry, spin_symmetry)
-    interaction = U * ud_num(particle_symmetry, spin_symmetry)
-    chemical_potential = -mu * e_num(particle_symmetry, spin_symmetry)
+function hubbard_hamiltonian(particle_symmetry, spin_symmetry; t, U, mu)
+    L = length(t) + 1
+    @assert length(t) + 1 == length(U) == length(mu)
+    hopping = e_hopping(particle_symmetry, spin_symmetry)
+    interaction = ud_num(particle_symmetry, spin_symmetry)
+    chemical_potential = e_num(particle_symmetry, spin_symmetry)
     I = id(hubbard_space(particle_symmetry, spin_symmetry))
     H = sum(1:(L - 1)) do i
-        return reduce(⊗, insert!(collect(Any, fill(I, L - 2)), i, hopping))
+        return reduce(⊗, insert!(collect(Any, fill(I, L - 2)), i, hopping * -t[i]))
     end +
         sum(1:L) do i
-        return reduce(⊗, insert!(collect(Any, fill(I, L - 1)), i, interaction))
+        return reduce(⊗, insert!(collect(Any, fill(I, L - 1)), i, interaction * U[i]))
     end +
         sum(1:L) do i
-        return reduce(⊗, insert!(collect(Any, fill(I, L - 1)), i, chemical_potential))
+        return reduce(⊗, insert!(collect(Any, fill(I, L - 1)), i, chemical_potential * -mu[i]))
+    end
+    return H
+end
+function hubbard_hamiltonian(::Type{SU2Irrep}, ::Type{SU2Irrep}; t, U, mu = U ./ 2)
+    L = length(t) + 1
+    @assert length(t) + 1 == length(U) == length(mu)
+    @assert mu ≈ U / 2
+    hopping = e_hopping(SU2Irrep, SU2Irrep)
+    interaction = half_ud_num(SU2Irrep, SU2Irrep)
+    I = id(hubbard_space(SU2Irrep, SU2Irrep))
+    H = sum(1:(L - 1)) do i
+        return reduce(⊗, insert!(collect(Any, fill(I, L - 2)), i, hopping * -t[i]))
+    end + sum(1:L) do i
+        return reduce(⊗, insert!(collect(Any, fill(I, L - 1)), i, interaction * U[i]))
     end
     return H
 end
 
 @testset "spectrum" begin
     L = 4
-    t = randn()
-    U = randn()
-    mu = randn()
+    t = randn(L - 1)
+    U = randn(L)
+    mu = randn(L)
 
-    H_triv = hubbard_hamiltonian(Trivial, Trivial; t, U, mu, L)
+    H_triv = hubbard_hamiltonian(Trivial, Trivial; t, U, mu)
     vals_triv = mapreduce(vcat, eigvals(H_triv)) do (c, v)
         return repeat(real.(v), dim(c))
     end
     sort!(vals_triv)
 
     for (particle_symmetry, spin_symmetry) in implemented_symmetries
-        if (particle_symmetry, spin_symmetry) == (Trivial, Trivial)
+        if particle_symmetry == spin_symmetry == Trivial ||
+                particle_symmetry == spin_symmetry == SU2Irrep
             continue
         end
-        H_symm = hubbard_hamiltonian(particle_symmetry, spin_symmetry; t, U, mu, L)
+        H_symm = hubbard_hamiltonian(particle_symmetry, spin_symmetry; t, U, mu)
         vals_symm = mapreduce(vcat, eigvals(H_symm)) do (c, v)
             return repeat(real.(v), dim(c))
         end
         sort!(vals_symm)
         @test vals_triv ≈ vals_symm
     end
+
+    mu = U ./ 2
+    H_triv = hubbard_hamiltonian(Trivial, Trivial; t, U, mu)
+    vals_triv = mapreduce(vcat, eigvals(H_triv)) do (c, v)
+        return repeat(real.(v), dim(c)) .+ sum(U) / 4
+    end
+    sort!(vals_triv)
+
+    H_symm = hubbard_hamiltonian(SU2Irrep, SU2Irrep; t, U, mu)
+    vals_symm = mapreduce(vcat, eigvals(H_symm)) do (c, v)
+        return repeat(real.(v), dim(c))
+    end
+    sort!(vals_symm)
+    @test vals_triv ≈ vals_symm
 end
 
 @testset "Exact diagonalisation" begin
@@ -215,6 +259,7 @@ end
             spin_symmetry in [Trivial, U1Irrep, SU2Irrep]
 
         if (particle_symmetry, spin_symmetry) in implemented_symmetries
+            particle_symmetry == spin_symmetry == SU2Irrep && continue
             rng = StableRNG(123)
 
             L = 2
@@ -222,7 +267,10 @@ end
             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)
+            H_triv = hubbard_hamiltonian(
+                particle_symmetry, spin_symmetry;
+                t = fill(t, L - 1), U = fill(U, L), mu = fill(mu, L)
+            )
 
             # Values based on https://arxiv.org/pdf/0807.4878. Introduction to Hubbard Model and Exact Diagonalization
             true_eigenvals = sort(