diff --git a/.gitignore b/.gitignore
index 06ba5eb..deec166 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,2 +1,4 @@
 docs/build
 Manifest.toml
+.vscode/
+.DS_Store
diff --git a/src/MPSKitModels.jl b/src/MPSKitModels.jl
index bf61ede..e46d026 100644
--- a/src/MPSKitModels.jl
+++ b/src/MPSKitModels.jl
@@ -35,6 +35,7 @@ export c⁺, c⁻, c⁺⁺, c⁻⁻, c⁺⁻, c⁻⁺
 export e_plus, e_min, e_plusplus, e_minmin, e_plusmin, e_minplus
 export e_number, e_number_up, e_number_down, e_number_updown
 export e⁺⁺, e⁻⁻, e⁺⁻, e⁻⁺
+export TJOperators
 
 export transverse_field_ising
 export kitaev_model
@@ -42,6 +43,7 @@ export quantum_potts
 export heisenberg_XXX, heisenberg_XXZ, heisenberg_XYZ
 export bilinear_biquadratic_model
 export hubbard_model, bose_hubbard_model
+export tj_model
 export quantum_chemistry_hamiltonian
 
 export classical_ising
@@ -65,6 +67,10 @@ include("operators/spinoperators.jl")
 include("operators/fermionoperators.jl")
 include("operators/hubbardoperators.jl")
 using .HubbardOperators
+# TJOperators share operator names with HubbardOperators
+# and is only imported to avoid name conflicts
+include("operators/tjoperators.jl")
+import .TJOperators
 include("operators/bosonoperators.jl")
 
 include("models/hamiltonians.jl")
diff --git a/src/models/hamiltonians.jl b/src/models/hamiltonians.jl
index c500588..1d441f2 100644
--- a/src/models/hamiltonians.jl
+++ b/src/models/hamiltonians.jl
@@ -370,3 +370,54 @@ function bose_hubbard_model(elt::Type{<:Number}=ComplexF64,
 
     return H
 end
+
+#===========================================================================================
+    t-J models
+===========================================================================================#
+
+"""
+    tj_model([elt::Type{<:Number}], [particle_symmetry::Type{<:Sector}],
+                  [spin_symmetry::Type{<:Sector}], [lattice::AbstractLattice];
+                  t, J, mu, slave_fermion::Bool=false)
+
+MPO for the hamiltonian of the t-J model, as defined by
+```math
+H = -t \\sum_{\\langle i,j \\rangle, \\sigma}
+    (\\tilde{e}^\\dagger_{i,\\sigma} \\tilde{e}_{j,\\sigma} + h.c.)
+    + J \\sum_{\\langle i,j \\rangle}(\\mathbf{S}_i \\cdot \\mathbf{S}_j - \\frac{1}{4} n_i n_j)
+    - \\mu \\sum_i n_i
+```
+where ``\\tilde{e}_{i,\\sigma}`` is the electron operator with spin ``\\sigma`` restrict to the no-double-occupancy subspace. 
+"""
+function tj_model end
+function tj_model(lattice::AbstractLattice; kwargs...)
+    return tj_model(ComplexF64, Trivial, Trivial, lattice; kwargs...)
+end
+function tj_model(particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector};
+                  kwargs...)
+    return tj_model(ComplexF64, particle_symmetry, spin_symmetry; kwargs...)
+end
+function tj_model(elt::Type{<:Number}, lattice::AbstractLattice; kwargs...)
+    return tj_model(elt, Trivial, Trivial, lattice; kwargs...)
+end
+function tj_model(T::Type{<:Number}=ComplexF64,
+                  particle_symmetry::Type{<:Sector}=Trivial,
+                  spin_symmetry::Type{<:Sector}=Trivial,
+                  lattice::AbstractLattice=InfiniteChain(1);
+                  t=2.5, J=1.0, mu=0.0, slave_fermion::Bool=false)
+    hopping = TJOperators.e_plusmin(T, particle_symmetry, spin_symmetry; slave_fermion) +
+              TJOperators.e_minplus(T, particle_symmetry, spin_symmetry; slave_fermion)
+    num = TJOperators.e_number(T, particle_symmetry, spin_symmetry; slave_fermion)
+    heisenberg = TJOperators.S_exchange(T, particle_symmetry, spin_symmetry;
+                                        slave_fermion) -
+                 (1 / 4) * (num ⊗ num)
+    return @mpoham begin
+        sum(nearest_neighbours(lattice)) do (i, j)
+            return (-t) * hopping{i,j} + J * heisenberg{i,j}
+        end + sum(vertices(lattice)) do i
+            return (-mu) * num{i}
+        end
+    end
+end
+
+# TODO: add (hardcore) bosonic t-J model (https://arxiv.org/abs/2409.15424)
diff --git a/src/operators/hubbardoperators.jl b/src/operators/hubbardoperators.jl
index 6b60c9f..8c4e259 100644
--- a/src/operators/hubbardoperators.jl
+++ b/src/operators/hubbardoperators.jl
@@ -64,7 +64,7 @@ end
 """
     e_plusmin_up(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector})
 
-Return the two-body operator that creates a spin-up electron at the first site and annihilates a spin-up electron at the second.
+Return the two-body operator ``e†_{1,↑}, e_{2,↑}`` that creates a spin-up electron at the first site and annihilates a spin-up electron at the second.
 """
 e_plusmin_up(P::Type{<:Sector}, S::Type{<:Sector}) = e_plusmin_up(ComplexF64, P, S)
 function e_plusmin_up(T, ::Type{Trivial}, ::Type{Trivial})
@@ -111,7 +111,7 @@ const e⁺e⁻ꜛ = e_plusmin_up
 """
     e_plusmin_down(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector})
 
-Return the two-body operator that creates a spin-down electron at the first site and annihilates a spin-down electron at the second.
+Return the two-body operator ``e†_{1,↓}, e_{2,↓}`` that creates a spin-down electron at the first site and annihilates a spin-down electron at the second.
 """
 e_plusmin_down(P::Type{<:Sector}, S::Type{<:Sector}) = e_plusmin_down(ComplexF64, P, S)
 function e_plusmin_down(T, ::Type{Trivial}, ::Type{Trivial})
@@ -158,7 +158,9 @@ const e⁺e⁻ꜜ = e_plusmin_down
 """
     e_minplus_up(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector})
 
-Return the two-body operator that annihilates a spin-up electron at the first site and creates a spin-up electron at the second.
+Return the Hermitian conjugate of `e_plusmin_up`, i.e.
+``(e†_{1,↑}, e_{2,↑})† = -e_{1,↑}, e†_{2,↑}`` (note the extra minus sign). 
+It annihilates a spin-up electron at the first site and creates a spin-up electron at the second.
 """
 e_minplus_up(P::Type{<:Sector}, S::Type{<:Sector}) = e_minplus_up(ComplexF64, P, S)
 function e_minplus_up(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector})
@@ -169,7 +171,9 @@ const e⁻⁺ꜛ = e_minplus_up
 """
     e_minplus_down(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector})
 
-Return the two-body operator that annihilates a spin-down electron at the first site and creates a spin-down electron at the second.
+Return the Hermitian conjugate of `e_plusmin_down`, i.e.
+``(e†_{1,↓}, e_{2,↓})† = -e_{1,↓}, e†_{2,↓}`` (note the extra minus sign). 
+It annihilates a spin-down electron at the first site and creates a spin-down electron at the second.
 """
 e_minplus_down(P::Type{<:Sector}, S::Type{<:Sector}) = e_minplus_down(ComplexF64, P, S)
 function e_minplus_down(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector})
diff --git a/src/operators/tjoperators.jl b/src/operators/tjoperators.jl
new file mode 100644
index 0000000..a93536f
--- /dev/null
+++ b/src/operators/tjoperators.jl
@@ -0,0 +1,629 @@
+#= Operators that act on t-J-type models
+i.e. the local hilbert space consists of 
+
+- usual basis states: 
+    |∅⟩, |↑⟩, |↓⟩
+- slave-fermion basis states (c_σ = h† b_σ; holon h is fermionic, spinon b_σ is bosonic): 
+    |h⟩ = h†|∅⟩, |↑'⟩ = (b↑)†|∅⟩, |↓'⟩ = (b↓)†|∅⟩
+=#
+module TJOperators
+
+using TensorKit
+
+export tj_space
+export e_plusmin, e_plusmin_up, e_plusmin_down
+export e_minplus, e_minplus_up, e_minplus_down
+export e_minmin_updown, e_minmin_downup, e_singlet
+export e_number, e_number_up, e_number_down, e_number_hole
+export S_x, S_y, S_z, S_plus, S_min
+export S_plusmin, S_minplus, S_exchange
+
+export e⁺e⁻, e⁺e⁻ꜛ, e⁺e⁻ꜜ, e⁻e⁺, e⁻e⁺ꜛ, e⁻e⁺ꜜ
+export nꜛ, nꜜ, nʰ
+export Sˣ, Sʸ, Sᶻ, S⁺, S⁻
+# not exported because namespace: export n
+
+"""
+    tj_space(particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; slave_fermion::Bool = false)
+
+Return the local hilbert space for a t-J-type model with the given particle and spin symmetries.
+The possible symmetries are 
+- Particle number: `Trivial`, `U1Irrep`;
+- Spin: `Trivial`, `U1Irrep`, `SU2Irrep`.
+
+Setting `slave_fermion = true` switches to the slave-fermion basis. 
+"""
+function tj_space(::Type{Trivial}=Trivial, ::Type{Trivial}=Trivial;
+                  slave_fermion::Bool=false)
+    return slave_fermion ? Vect[FermionParity](0 => 2, 1 => 1) :
+           Vect[FermionParity](0 => 1, 1 => 2)
+end
+function tj_space(::Type{Trivial}, ::Type{U1Irrep}; slave_fermion::Bool=false)
+    return if slave_fermion
+        Vect[FermionParity ⊠ U1Irrep]((1, 0) => 1, (0, 1 // 2) => 1, (0, -1 // 2) => 1)
+    else
+        Vect[FermionParity ⊠ U1Irrep]((0, 0) => 1, (1, 1 // 2) => 1, (1, -1 // 2) => 1)
+    end
+end
+function tj_space(::Type{Trivial}, ::Type{SU2Irrep}; slave_fermion::Bool=false)
+    return error("Not implemented")
+end
+function tj_space(::Type{U1Irrep}, ::Type{Trivial}; slave_fermion::Bool=false)
+    return if slave_fermion
+        Vect[FermionParity ⊠ U1Irrep]((1, 0) => 1, (0, 1) => 2)
+    else
+        Vect[FermionParity ⊠ U1Irrep]((0, 0) => 1, (1, 1) => 2)
+    end
+end
+function tj_space(::Type{U1Irrep}, ::Type{U1Irrep}; slave_fermion::Bool=false)
+    return if slave_fermion
+        Vect[FermionParity ⊠ U1Irrep ⊠ U1Irrep]((1, 0, 0) => 1, (0, 1, 1 // 2) => 1,
+                                                (0, 1, -1 // 2) => 1)
+    else
+        Vect[FermionParity ⊠ U1Irrep ⊠ U1Irrep]((0, 0, 0) => 1, (1, 1, 1 // 2) => 1,
+                                                (1, 1, -1 // 2) => 1)
+    end
+end
+function tj_space(::Type{U1Irrep}, ::Type{SU2Irrep}; slave_fermion::Bool=false)
+    return error("Not implemented")
+end
+
+function single_site_operator(T, particle_symmetry::Type{<:Sector},
+                              spin_symmetry::Type{<:Sector}; slave_fermion::Bool=false)
+    V = tj_space(particle_symmetry, spin_symmetry; slave_fermion)
+    return TensorMap(zeros, T, V ← V)
+end
+
+function two_site_operator(T, particle_symmetry::Type{<:Sector},
+                           spin_symmetry::Type{<:Sector}; slave_fermion::Bool=false)
+    V = tj_space(particle_symmetry, spin_symmetry; slave_fermion)
+    return TensorMap(zeros, T, V ⊗ V ← V ⊗ V)
+end
+
+"""
+    e_plusmin_up(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; slave_fermion::Bool = false)
+
+Return the two-body operator ``e†_{1,↑}, e_{2,↑}`` that creates a spin-up electron at the first site and annihilates a spin-up electron at the second.
+The only nonzero matrix element corresponds to `|↑0⟩ <-- |0↑⟩`.
+"""
+function e_plusmin_up(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool=false)
+    return e_plusmin_up(ComplexF64, P, S; slave_fermion)
+end
+function e_plusmin_up(T, ::Type{Trivial}, ::Type{Trivial}; slave_fermion::Bool=false)
+    t = two_site_operator(T, Trivial, Trivial; slave_fermion)
+    I = sectortype(t)
+    (h, b, sgn) = slave_fermion ? (1, 0, -1) : (0, 1, 1)
+    #= The extra minus sign in slave-fermion basis:
+    c†_{1,↑} c_{2,↑} |0↑⟩
+    = h_1 b†_{1,↑} h†_2 b_{2,↑} h†_1 b†_{2,↑}|vac⟩
+    = -b†_{1,↑} h†_2 h_1 h†_1 b_{2,↑} b†_{2,↑}|vac⟩
+    = -b†_{1,↑} h†_2 |vac⟩
+    = -|↑0⟩
+    =#
+    t[(I(b), I(h), dual(I(h)), dual(I(b)))][1, 1, 1, 1] = sgn * 1
+    return t
+end
+function e_plusmin_up(T, ::Type{Trivial}, ::Type{U1Irrep}; slave_fermion::Bool=false)
+    t = two_site_operator(T, Trivial, U1Irrep; slave_fermion)
+    I = sectortype(t)
+    (h, b, sgn) = slave_fermion ? (1, 0, -1) : (0, 1, 1)
+    t[(I(b, 1 // 2), I(h, 0), dual(I(h, 0)), dual(I(b, 1 // 2)))] .= sgn * 1
+    return t
+end
+function e_plusmin_up(T, ::Type{Trivial}, ::Type{SU2Irrep}; slave_fermion::Bool=false)
+    return error("Not implemented")
+end
+function e_plusmin_up(T, ::Type{U1Irrep}, ::Type{Trivial}; slave_fermion::Bool=false)
+    t = two_site_operator(T, U1Irrep, Trivial; slave_fermion)
+    I = sectortype(t)
+    (h, b, sgn) = slave_fermion ? (1, 0, -1) : (0, 1, 1)
+    t[(I(b, 1), I(h, 0), dual(I(h, 0)), dual(I(b, 1)))][1, 1, 1, 1] = sgn * 1
+    return t
+end
+function e_plusmin_up(T, ::Type{U1Irrep}, ::Type{U1Irrep}; slave_fermion::Bool=false)
+    t = two_site_operator(T, U1Irrep, U1Irrep; slave_fermion)
+    I = sectortype(t)
+    (h, b, sgn) = slave_fermion ? (1, 0, -1) : (0, 1, 1)
+    t[(I(b, 1, 1 // 2), I(h, 0, 0), dual(I(h, 0, 0)), dual(I(b, 1, 1 // 2)))] .= sgn * 1
+    return t
+end
+function e_plusmin_up(T, ::Type{U1Irrep}, ::Type{SU2Irrep}; slave_fermion::Bool=false)
+    return error("Not implemented")
+end
+const e⁺e⁻ꜛ = e_plusmin_up
+
+"""
+    e_plusmin_down(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; slave_fermion::Bool = false)
+
+Return the two-body operator ``e†_{1,↓}, e_{2,↓}`` that creates a spin-down electron at the first site and annihilates a spin-down electron at the second.
+The only nonzero matrix element corresponds to `|↓0⟩ <-- |0↓⟩`.
+"""
+function e_plusmin_down(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool=false)
+    return e_plusmin_down(ComplexF64, P, S; slave_fermion)
+end
+function e_plusmin_down(T, ::Type{Trivial}, ::Type{Trivial}; slave_fermion::Bool=false)
+    t = two_site_operator(T, Trivial, Trivial; slave_fermion)
+    I = sectortype(t)
+    (h, b, sgn) = slave_fermion ? (1, 0, -1) : (0, 1, 1)
+    t[(I(b), I(h), dual(I(h)), dual(I(b)))][2, 1, 1, 2] = sgn * 1
+    return t
+end
+function e_plusmin_down(T, ::Type{Trivial}, ::Type{U1Irrep}; slave_fermion::Bool=false)
+    t = two_site_operator(T, Trivial, U1Irrep; slave_fermion)
+    I = sectortype(t)
+    (h, b, sgn) = slave_fermion ? (1, 0, -1) : (0, 1, 1)
+    t[(I(b, -1 // 2), I(h, 0), dual(I(h, 0)), dual(I(b, -1 // 2)))] .= sgn * 1
+    return t
+end
+function e_plusmin_down(T, ::Type{Trivial}, ::Type{SU2Irrep}; slave_fermion::Bool=false)
+    return error("Not implemented")
+end
+function e_plusmin_down(T, ::Type{U1Irrep}, ::Type{Trivial}; slave_fermion::Bool=false)
+    t = two_site_operator(T, U1Irrep, Trivial; slave_fermion)
+    I = sectortype(t)
+    (h, b, sgn) = slave_fermion ? (1, 0, -1) : (0, 1, 1)
+    t[(I(b, 1), I(h, 0), dual(I(h, 0)), dual(I(b, 1)))][2, 1, 1, 2] = sgn * 1
+    return t
+end
+function e_plusmin_down(T, ::Type{U1Irrep}, ::Type{U1Irrep}; slave_fermion::Bool=false)
+    t = two_site_operator(T, U1Irrep, U1Irrep; slave_fermion)
+    I = sectortype(t)
+    (h, b, sgn) = slave_fermion ? (1, 0, -1) : (0, 1, 1)
+    t[(I(b, 1, -1 // 2), I(h, 0, 0), dual(I(h, 0, 0)), dual(I(b, 1, -1 // 2)))] .= sgn * 1
+    return t
+end
+function e_plusmin_down(T, ::Type{U1Irrep}, ::Type{SU2Irrep}; slave_fermion::Bool=false)
+    return error("Not implemented")
+end
+const e⁺e⁻ꜜ = e_plusmin_down
+
+"""
+    e_minplus_up(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; slave_fermion::Bool = false)
+
+Return the Hermitian conjugate of `e_plusmin_up`, i.e.
+``(e†_{1,↑}, e_{2,↑})† = -e_{1,↑}, e†_{2,↑}`` (note the extra minus sign). 
+It annihilates a spin-up electron at the first site and creates a spin-up electron at the second.
+The only nonzero matrix element corresponds to `|0↑⟩ <-- |↑0⟩`.
+"""
+function e_minplus_up(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool=false)
+    return e_minplus_up(ComplexF64, P, S; slave_fermion)
+end
+function e_minplus_up(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector};
+                      slave_fermion::Bool=false)
+    return copy(adjoint(e_plusmin_up(T, particle_symmetry, spin_symmetry; slave_fermion)))
+end
+const e⁻⁺ꜛ = e_minplus_up
+
+"""
+    e_minplus_down(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; slave_fermion::Bool = false)
+
+Return the Hermitian conjugate of `e_plusmin_down`, i.e.
+``(e†_{1,↓}, e_{2,↓})† = -e_{1,↓}, e†_{2,↓}`` (note the extra minus sign). 
+It annihilates a spin-down electron at the first site and creates a spin-down electron at the second.
+The only nonzero matrix element corresponds to `|0↓⟩ <-- |↓0⟩`.
+"""
+function e_minplus_down(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool=false)
+    return e_minplus_down(ComplexF64, P, S; slave_fermion)
+end
+function e_minplus_down(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector};
+                        slave_fermion::Bool=false)
+    return copy(adjoint(e_plusmin_down(T, particle_symmetry, spin_symmetry; slave_fermion)))
+end
+const e⁻e⁺ꜜ = e_minplus_down
+
+"""
+    e_plusmin(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; slave_fermion::Bool = false)
+
+Return the two-body operator that creates a particle at the first site and annihilates a particle at the second.
+This is the sum of `e_plusmin_up` and `e_plusmin_down`.
+"""
+function e_plusmin(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool=false)
+    return e_plusmin(ComplexF64, P, S; slave_fermion)
+end
+function e_plusmin(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector};
+                   slave_fermion::Bool=false)
+    return e_plusmin_up(T, particle_symmetry, spin_symmetry; slave_fermion) +
+           e_plusmin_down(T, particle_symmetry, spin_symmetry; slave_fermion)
+end
+const e⁺e⁻ = e_plusmin
+
+"""
+    e_minplus(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; slave_fermion::Bool = false)
+
+Return the two-body operator that annihilates a particle at the first site and creates a particle at the second.
+This is the sum of `e_minplus_up` and `e_minplus_down`.
+"""
+function e_minplus(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool=false)
+    return e_minplus(ComplexF64, P, S; slave_fermion)
+end
+function e_minplus(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector};
+                   slave_fermion::Bool=false)
+    return copy(adjoint(e_plusmin(T, particle_symmetry, spin_symmetry; slave_fermion)))
+end
+const e⁻e⁺ = e_minplus
+
+"""
+    e_minmin_updown(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; slave_fermion::Bool = false)
+
+Return the two-body operator ``e_{1,↑} e_{2,↓}`` that annihilates a spin-up particle at the first site and a spin-down particle at the second site.
+The only nonzero matrix element corresponds to `|00⟩ <-- |↑↓⟩`.
+"""
+function e_minmin_updown(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool=false)
+    return e_minmin_updown(ComplexF64, P, S; slave_fermion)
+end
+function e_minmin_updown(T, ::Type{Trivial}, ::Type{Trivial}; slave_fermion::Bool=false)
+    t = two_site_operator(T, Trivial, Trivial; slave_fermion)
+    I = sectortype(t)
+    (h, b, sgn) = slave_fermion ? (1, 0, -1) : (0, 1, 1)
+    t[(I(h), I(h), dual(I(b)), dual(I(b)))][1, 1, 1, 2] = -sgn * 1
+    return t
+end
+function e_minmin_updown(T, ::Type{Trivial}, ::Type{U1Irrep}; slave_fermion::Bool=false)
+    t = two_site_operator(T, Trivial, U1Irrep; slave_fermion)
+    I = sectortype(t)
+    (h, b, sgn) = slave_fermion ? (1, 0, -1) : (0, 1, 1)
+    t[(I(h, 0), I(h, 0), dual(I(b, 1 // 2)), dual(I(b, -1 // 2)))] .= -sgn * 1
+    return t
+end
+function e_minmin_updown(T, ::Type{U1Irrep}, ::Type{<:Sector}; slave_fermion::Bool=false)
+    throw(ArgumentError("`e_minmin_updown` is not symmetric under `U1Irrep` particle symmetry"))
+end
+
+"""
+    e_minmin_downup(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; slave_fermion::Bool = false)
+
+Return the two-body operator ``e_{1,↓} e_{2,↑}`` that annihilates a spin-down particle at the first site and a spin-up particle at the second site.
+The only nonzero matrix element corresponds to `|00⟩ <-- |↓↑⟩`.
+"""
+function e_minmin_downup(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool=false)
+    return e_minmin_downup(ComplexF64, P, S; slave_fermion)
+end
+function e_minmin_downup(T, ::Type{Trivial}, ::Type{Trivial}; slave_fermion::Bool=false)
+    t = two_site_operator(T, Trivial, Trivial; slave_fermion)
+    I = sectortype(t)
+    (h, b, sgn) = slave_fermion ? (1, 0, -1) : (0, 1, 1)
+    t[(I(h), I(h), dual(I(b)), dual(I(b)))][1, 1, 2, 1] = -sgn * 1
+    return t
+end
+function e_minmin_downup(T, ::Type{Trivial}, ::Type{U1Irrep}; slave_fermion::Bool=false)
+    t = two_site_operator(T, Trivial, U1Irrep; slave_fermion)
+    I = sectortype(t)
+    (h, b, sgn) = slave_fermion ? (1, 0, -1) : (0, 1, 1)
+    t[(I(h, 0), I(h, 0), dual(I(b, -1 // 2)), dual(I(b, 1 // 2)))] .= -sgn * 1
+    return t
+end
+function e_minmin_downup(T, ::Type{U1Irrep}, ::Type{<:Sector}; slave_fermion::Bool=false)
+    throw(ArgumentError("`e_minmin_downup` is not symmetric under `U1Irrep` particle symmetry"))
+end
+
+"""
+    e_singlet(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; slave_fermion::Bool = false)
+
+Return the two-body singlet operator ``(e_{1,↓} e_{2,↑} - e_{1,↓} e_{2,↑}) / sqrt(2)``.
+"""
+function e_singlet(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool=false)
+    return e_singlet(ComplexF64, P, S; slave_fermion)
+end
+function e_singlet(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector};
+                   slave_fermion::Bool=false)
+    return (e_minmin_updown(T, particle_symmetry, spin_symmetry; slave_fermion) -
+            e_minmin_downup(T, particle_symmetry, spin_symmetry; slave_fermion)) / sqrt(2)
+end
+
+"""
+    e_number_up(particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; slave_fermion::Bool = false)
+
+Return the one-body operator that counts the number of spin-up electrons.
+"""
+function e_number_up(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool=false)
+    return e_number_up(ComplexF64, P, S; slave_fermion)
+end
+function e_number_up(T::Type{<:Number}, ::Type{Trivial}, ::Type{Trivial};
+                     slave_fermion::Bool=false)
+    t = single_site_operator(T, Trivial, Trivial; slave_fermion)
+    I = sectortype(t)
+    b = slave_fermion ? 0 : 1
+    t[(I(b), dual(I(b)))][1, 1] = 1
+    return t
+end
+function e_number_up(T, ::Type{Trivial}, ::Type{U1Irrep}; slave_fermion::Bool=false)
+    t = single_site_operator(T, Trivial, U1Irrep; slave_fermion)
+    I = sectortype(t)
+    b = slave_fermion ? 0 : 1
+    t[(I(b, 1 // 2), dual(I(b, 1 // 2)))][1, 1] = 1
+    return t
+end
+function e_number_up(T, ::Type{Trivial}, ::Type{SU2Irrep}; slave_fermion::Bool=false)
+    throw(ArgumentError("`e_number_up` is not symmetric under `SU2Irrep` spin symmetry"))
+end
+function e_number_up(T, ::Type{U1Irrep}, ::Type{Trivial}; slave_fermion::Bool=false)
+    t = single_site_operator(T, U1Irrep, Trivial; slave_fermion)
+    I = sectortype(t)
+    b = slave_fermion ? 0 : 1
+    t[(I(b, 1), dual(I(b, 1)))][1, 1] = 1
+    return t
+end
+function e_number_up(T, ::Type{U1Irrep}, ::Type{U1Irrep}; slave_fermion::Bool=false)
+    t = single_site_operator(T, U1Irrep, U1Irrep; slave_fermion)
+    I = sectortype(t)
+    b = slave_fermion ? 0 : 1
+    t[(I(b, 1, 1 // 2), dual(I(b, 1, 1 // 2)))] .= 1
+    return t
+end
+function e_number_up(T, ::Type{U1Irrep}, ::Type{SU2Irrep}; slave_fermion::Bool=false)
+    throw(ArgumentError("`e_number_up` is not symmetric under `SU2Irrep` spin symmetry"))
+end
+const nꜛ = e_number_up
+
+"""
+    e_number_down(particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; slave_fermion::Bool=false)
+
+Return the one-body operator that counts the number of spin-down electrons.
+"""
+function e_number_down(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool=false)
+    return e_number_down(ComplexF64, P, S; slave_fermion)
+end
+function e_number_down(T::Type{<:Number}, ::Type{Trivial}, ::Type{Trivial};
+                       slave_fermion::Bool=false)
+    t = single_site_operator(T, Trivial, Trivial; slave_fermion)
+    I = sectortype(t)
+    b = slave_fermion ? 0 : 1
+    t[(I(b), dual(I(b)))][2, 2] = 1
+    return t
+end
+function e_number_down(T, ::Type{Trivial}, ::Type{U1Irrep}; slave_fermion::Bool=false)
+    t = single_site_operator(T, Trivial, U1Irrep; slave_fermion)
+    I = sectortype(t)
+    b = slave_fermion ? 0 : 1
+    t[(I(b, -1 // 2), dual(I(b, -1 // 2)))][1, 1] = 1
+    return t
+end
+function e_number_down(T, ::Type{Trivial}, ::Type{SU2Irrep}; slave_fermion::Bool=false)
+    throw(ArgumentError("`e_number_down` is not symmetric under `SU2Irrep` spin symmetry"))
+end
+function e_number_down(T, ::Type{U1Irrep}, ::Type{Trivial}; slave_fermion::Bool=false)
+    t = single_site_operator(T, U1Irrep, Trivial; slave_fermion)
+    I = sectortype(t)
+    b = slave_fermion ? 0 : 1
+    t[(I(b, 1), dual(I(b, 1)))][2, 2] = 1
+    return t
+end
+function e_number_down(T, ::Type{U1Irrep}, ::Type{U1Irrep}; slave_fermion::Bool=false)
+    t = single_site_operator(T, U1Irrep, U1Irrep; slave_fermion)
+    I = sectortype(t)
+    b = slave_fermion ? 0 : 1
+    t[(I(b, 1, -1 // 2), dual(I(b, 1, -1 // 2)))] .= 1
+    return t
+end
+function e_number_down(T, ::Type{U1Irrep}, ::Type{SU2Irrep}; slave_fermion::Bool=false)
+    throw(ArgumentError("`e_number_down` is not symmetric under `SU2Irrep` spin symmetry"))
+end
+const nꜜ = e_number_down
+
+"""
+    e_number(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; slave_fermion::Bool=false)
+
+Return the one-body operator that counts the number of particles.
+"""
+function e_number(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool=false)
+    return e_number(ComplexF64, P, S; slave_fermion)
+end
+function e_number(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector};
+                  slave_fermion::Bool=false)
+    return e_number_up(T, particle_symmetry, spin_symmetry; slave_fermion) +
+           e_number_down(T, particle_symmetry, spin_symmetry; slave_fermion)
+end
+const n = e_number
+
+"""
+    e_number_hole(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; slave_fermion::Bool=false)
+
+Return the one-body operator that counts the number of holes.
+"""
+function e_number_hole(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool=false)
+    return e_number_hole(ComplexF64, P, S; slave_fermion)
+end
+function e_number_hole(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector};
+                       slave_fermion::Bool=false)
+    iden = TensorKit.id(tj_space(particle_symmetry, spin_symmetry; slave_fermion))
+    return iden - e_number(T, particle_symmetry, spin_symmetry; slave_fermion)
+end
+const nʰ = e_number_hole
+
+"""
+    S_x(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; slave_fermion::Bool=false)
+
+Return the one-body spin-1/2 x-operator on the electrons.
+"""
+function S_x(P::Type{<:Sector}=Trivial, S::Type{<:Sector}=Trivial;
+             slave_fermion::Bool=false)
+    return S_x(ComplexF64, P, S; slave_fermion)
+end
+function S_x(T::Type{<:Number}, ::Type{Trivial}, ::Type{Trivial}; slave_fermion::Bool=false)
+    t = single_site_operator(T, Trivial, Trivial; slave_fermion)
+    I = sectortype(t)
+    b = slave_fermion ? 0 : 1
+    t[(I(b), dual(I(b)))][1, 2] = 0.5
+    t[(I(b), dual(I(b)))][2, 1] = 0.5
+    return t
+end
+function S_x(T::Type{<:Number}, ::Type{U1Irrep}, ::Type{Trivial}; slave_fermion::Bool=false)
+    t = single_site_operator(T, U1Irrep, Trivial; slave_fermion)
+    I = sectortype(t)
+    b = slave_fermion ? 0 : 1
+    t[(I(b, 1), dual(I(b, 1)))][1, 2] = 0.5
+    t[(I(b, 1), dual(I(b, 1)))][2, 1] = 0.5
+    return t
+end
+const Sˣ = S_x
+
+"""
+    S_y(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; slave_fermion::Bool=false)
+
+Return the one-body spin-1/2 x-operator on the electrons (only defined for `Trivial` symmetry). 
+"""
+function S_y(P::Type{<:Sector}=Trivial, S::Type{<:Sector}=Trivial;
+             slave_fermion::Bool=false)
+    return S_y(ComplexF64, P, S; slave_fermion)
+end
+function S_y(T::Type{<:Number}, ::Type{Trivial}, ::Type{Trivial}; slave_fermion::Bool=false)
+    t = single_site_operator(T, Trivial, Trivial; slave_fermion)
+    I = sectortype(t)
+    b = slave_fermion ? 0 : 1
+    t[(I(b), dual(I(b)))][1, 2] = -0.5im
+    t[(I(b), dual(I(b)))][2, 1] = 0.5im
+    return t
+end
+function S_y(T::Type{<:Number}, ::Type{U1Irrep}, ::Type{Trivial}; slave_fermion::Bool=false)
+    t = single_site_operator(T, U1Irrep, Trivial; slave_fermion)
+    I = sectortype(t)
+    b = slave_fermion ? 0 : 1
+    t[(I(b, 1), dual(I(b, 1)))][1, 2] = -0.5im
+    t[(I(b, 1), dual(I(b, 1)))][2, 1] = 0.5im
+    return t
+end
+const Sʸ = S_y
+
+"""
+    S_z(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; slave_fermion::Bool=false)
+
+Return the one-body spin-1/2 z-operator on the electrons. 
+"""
+function S_z(P::Type{<:Sector}=Trivial, S::Type{<:Sector}=Trivial;
+             slave_fermion::Bool=false)
+    return S_z(ComplexF64, P, S; slave_fermion)
+end
+function S_z(T::Type{<:Number}, ::Type{Trivial}, ::Type{Trivial}; slave_fermion::Bool=false)
+    t = single_site_operator(T, Trivial, Trivial; slave_fermion)
+    I = sectortype(t)
+    b = slave_fermion ? 0 : 1
+    t[(I(b), dual(I(b)))][1, 1] = 0.5
+    t[(I(b), dual(I(b)))][2, 2] = -0.5
+    return t
+end
+function S_z(T::Type{<:Number}, ::Type{Trivial}, ::Type{U1Irrep}; slave_fermion::Bool=false)
+    t = single_site_operator(T, Trivial, U1Irrep; slave_fermion)
+    I = sectortype(t)
+    b = slave_fermion ? 0 : 1
+    t[(I(b, 1 // 2), dual(I(b, 1 // 2)))] .= 0.5
+    t[(I(b, -1 // 2), dual(I(b, -1 // 2)))] .= -0.5
+    return t
+end
+function S_z(T::Type{<:Number}, ::Type{U1Irrep}, ::Type{Trivial}; slave_fermion::Bool=false)
+    t = single_site_operator(T, U1Irrep, Trivial; slave_fermion)
+    I = sectortype(t)
+    b = slave_fermion ? 0 : 1
+    t[(I(b, 1), dual(I(b, 1)))][1, 1] = 0.5
+    t[(I(b, 1), dual(I(b, 1)))][2, 2] = -0.5
+    return t
+end
+function S_z(T::Type{<:Number}, ::Type{U1Irrep}, ::Type{U1Irrep}; slave_fermion::Bool=false)
+    t = single_site_operator(T, U1Irrep, U1Irrep; slave_fermion)
+    I = sectortype(t)
+    b = slave_fermion ? 0 : 1
+    t[(I(b, 1, 1 // 2), dual(I(b, 1, 1 // 2)))] .= 0.5
+    t[(I(b, 1, -1 // 2), dual(I(b, 1, -1 // 2)))] .= -0.5
+    return t
+end
+const Sᶻ = S_z
+
+"""
+    S_plus(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; slave_fermion::Bool=false)
+
+Return the spin-plus operator.
+"""
+function S_plus(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool=false)
+    return S_plus(ComplexF64, P, S; slave_fermion)
+end
+function S_plus(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector};
+                slave_fermion::Bool=false)
+    return S_x(T, particle_symmetry, spin_symmetry; slave_fermion) +
+           1im * S_y(T, particle_symmetry, spin_symmetry; slave_fermion)
+end
+const S⁺ = S_plus
+
+"""
+    S_min(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; slave_fermion::Bool=false)
+
+Return the spin-minus operator.
+"""
+function S_min(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool=false)
+    return S_min(ComplexF64, P, S; slave_fermion)
+end
+function S_min(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector};
+               slave_fermion::Bool=false)
+    return S_x(T, particle_symmetry, spin_symmetry; slave_fermion) -
+           1im * S_y(T, particle_symmetry, spin_symmetry; slave_fermion)
+end
+const S⁻ = S_min
+
+"""
+    S_plusmin(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; slave_fermion::Bool = false)
+
+Return the two-body operator S⁺S⁻.
+The only nonzero matrix element corresponds to `|↑↓⟩ <-- |↓↑⟩`.
+"""
+function S_plusmin(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool=false)
+    return S_plusmin(ComplexF64, P, S; slave_fermion)
+end
+function S_plusmin(T, ::Type{Trivial}, ::Type{Trivial}; slave_fermion::Bool=false)
+    t = two_site_operator(T, Trivial, Trivial; slave_fermion)
+    I = sectortype(t)
+    b = slave_fermion ? 0 : 1
+    t[(I(b), I(b), dual(I(b)), dual(I(b)))][1, 2, 2, 1] = 1
+    return t
+end
+function S_plusmin(T, ::Type{Trivial}, ::Type{U1Irrep}; slave_fermion::Bool=false)
+    t = two_site_operator(T, Trivial, U1Irrep; slave_fermion)
+    I = sectortype(t)
+    b = slave_fermion ? 0 : 1
+    t[(I(b, 1 // 2), I(b, -1 // 2), dual(I(b, -1 // 2)), dual(I(b, 1 // 2)))] .= 1
+    return t
+end
+function S_plusmin(T, ::Type{U1Irrep}, ::Type{Trivial}; slave_fermion::Bool=false)
+    t = two_site_operator(T, U1Irrep, Trivial; slave_fermion)
+    I = sectortype(t)
+    b = slave_fermion ? 0 : 1
+    t[(I(b, 1), I(b, 1), dual(I(b, 1)), dual(I(b, 1)))][1, 2, 2, 1] = 1
+    return t
+end
+function S_plusmin(T, ::Type{U1Irrep}, ::Type{U1Irrep}; slave_fermion::Bool=false)
+    t = two_site_operator(T, U1Irrep, U1Irrep; slave_fermion)
+    I = sectortype(t)
+    b = slave_fermion ? 0 : 1
+    t[(I(b, 1, 1 // 2), I(b, 1, -1 // 2), dual(I(b, 1, -1 // 2)), dual(I(b, 1, 1 // 2)))] .= 1
+    return t
+end
+
+"""
+    S_minplus(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; slave_fermion::Bool = false)
+
+Return the two-body operator S⁻S⁺.
+The only nonzero matrix element corresponds to `|↓↑⟩ <-- |↑↓⟩`.
+"""
+function S_minplus(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool=false)
+    return S_minplus(ComplexF64, P, S; slave_fermion)
+end
+function S_minplus(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector};
+                   slave_fermion::Bool=false)
+    return copy(adjoint(S_plusmin(T, particle_symmetry, spin_symmetry; slave_fermion)))
+end
+
+"""
+    S_exchange(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; slave_fermion::Bool = false)
+
+Return the spin exchange operator S⋅S.
+"""
+function S_exchange(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool=false)
+    return S_exchange(ComplexF64, P, S; slave_fermion)
+end
+function S_exchange(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector};
+                    slave_fermion::Bool=false)
+    Sz = S_z(T, particle_symmetry, spin_symmetry; slave_fermion)
+    return (1 / 2) * (S_plusmin(T, particle_symmetry, spin_symmetry; slave_fermion)
+                      +
+                      S_minplus(T, particle_symmetry, spin_symmetry; slave_fermion)) +
+           Sz ⊗ Sz
+end
+
+end
diff --git a/test/hubbardoperators.jl b/test/hubbardoperators.jl
index 105ac22..1a92d87 100644
--- a/test/hubbardoperators.jl
+++ b/test/hubbardoperators.jl
@@ -3,9 +3,7 @@ using TensorKit
 using MPSKitModels.HubbardOperators
 using LinearAlgebra: eigvals
 
-implemented_symmetries = [(Trivial, Trivial),
-                          (U1Irrep, U1Irrep),
-                          (U1Irrep, SU2Irrep)]
+implemented_symmetries = [(Trivial, Trivial), (U1Irrep, U1Irrep), (U1Irrep, SU2Irrep)]
 @testset "basic properties" begin
     for particle_symmetry in (Trivial, U1Irrep, SU2Irrep),
         spin_symmetry in (Trivial, U1Irrep, SU2Irrep)
@@ -51,7 +49,8 @@ function hamiltonian(particle_symmetry, spin_symmetry; t, U, mu, L)
     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))
-        end + sum(1:L) do i
+        end +
+        sum(1:L) do i
             return reduce(⊗, insert!(collect(Any, fill(I, L - 1)), i, interaction))
         end +
         sum(1:L) do i
@@ -66,28 +65,21 @@ end
     U = randn()
     mu = randn()
 
-    hopping = t * (e_plusmin(Trivial, Trivial) + e_minplus(Trivial, Trivial))
-    interaction = U * e_number_updown(Trivial, Trivial)
-    chemical_potential = mu * e_number(Trivial, Trivial)
-    I = id(domain(interaction))
-
     H_triv = hamiltonian(Trivial, Trivial; t, U, mu, L)
     vals_triv = mapreduce(vcat, eigvals(H_triv)) do (c, v)
         return repeat(real.(v), dim(c))
     end
     sort!(vals_triv)
 
-    H_u1_u1 = hamiltonian(U1Irrep, U1Irrep; t, U, mu, L)
-    vals_u1_u1 = mapreduce(vcat, eigvals(H_u1_u1)) do (c, v)
-        return repeat(real.(v), dim(c))
-    end
-    sort!(vals_u1_u1)
-    @test vals_triv ≈ vals_u1_u1
-
-    H_u1_su2 = hamiltonian(U1Irrep, SU2Irrep; t, U, mu, L)
-    vals_u1_su2 = mapreduce(vcat, eigvals(H_u1_su2)) do (c, v)
-        return repeat(real.(v), dim(c))
+    for (particle_symmetry, spin_symmetry) in implemented_symmetries
+        if (particle_symmetry, spin_symmetry) == (Trivial, Trivial)
+            continue
+        end
+        H_symm = hamiltonian(particle_symmetry, spin_symmetry; t, U, mu, L)
+        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
-    sort!(vals_u1_su2)
-    @test vals_triv ≈ vals_u1_su2
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index c897426..08e9b52 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -21,6 +21,10 @@ end
     include("hubbardoperators.jl")
 end
 
+@testset "t-J operators" begin
+    include("tjoperators.jl")
+end
+
 @testset "mpoham" begin
     include("mpoham.jl")
 end
diff --git a/test/tjoperators.jl b/test/tjoperators.jl
new file mode 100644
index 0000000..fa71421
--- /dev/null
+++ b/test/tjoperators.jl
@@ -0,0 +1,120 @@
+using Test
+using TensorKit
+using LinearAlgebra: eigvals
+import MPSKitModels: TJOperators as tJ
+
+implemented_symmetries = [(Trivial, Trivial),
+                          (Trivial, U1Irrep),
+                          (U1Irrep, Trivial),
+                          (U1Irrep, U1Irrep)]
+@testset "basic properties" begin
+    for slave_fermion in (false, true),
+        particle_symmetry in (Trivial, U1Irrep),
+        spin_symmetry in (Trivial, U1Irrep, SU2Irrep)
+
+        if (particle_symmetry, spin_symmetry) in implemented_symmetries
+            # test hermiticity
+            @test tJ.e_plusmin(particle_symmetry, spin_symmetry; slave_fermion)' ≈
+                  tJ.e_minplus(particle_symmetry, spin_symmetry; slave_fermion)
+            if spin_symmetry !== SU2Irrep
+                @test tJ.e_plusmin_down(particle_symmetry, spin_symmetry; slave_fermion)' ≈
+                      tJ.e_minplus_down(particle_symmetry, spin_symmetry; slave_fermion)
+                @test tJ.e_plusmin_up(particle_symmetry, spin_symmetry; slave_fermion)' ≈
+                      tJ.e_minplus_up(particle_symmetry, spin_symmetry; slave_fermion)
+                @test tJ.e_plusmin_down(particle_symmetry, spin_symmetry; slave_fermion)' ≈
+                      tJ.e_minplus_down(particle_symmetry, spin_symmetry; slave_fermion)
+                @test tJ.e_plusmin_up(particle_symmetry, spin_symmetry; slave_fermion)' ≈
+                      tJ.e_minplus_up(particle_symmetry, spin_symmetry; slave_fermion)
+            end
+
+            # test number operator
+            if spin_symmetry !== SU2Irrep
+                pspace = tJ.tj_space(particle_symmetry, spin_symmetry; slave_fermion)
+                @test tJ.e_number(particle_symmetry, spin_symmetry; slave_fermion) ≈
+                      tJ.e_number_up(particle_symmetry, spin_symmetry; slave_fermion) +
+                      tJ.e_number_down(particle_symmetry, spin_symmetry; slave_fermion)
+                @test TensorMap(zeros, pspace, pspace) ≈
+                      tJ.e_number_up(particle_symmetry, spin_symmetry; slave_fermion) *
+                      tJ.e_number_down(particle_symmetry, spin_symmetry; slave_fermion) ≈
+                      tJ.e_number_down(particle_symmetry, spin_symmetry; slave_fermion) *
+                      tJ.e_number_up(particle_symmetry, spin_symmetry; slave_fermion)
+            end
+
+            # test spin operator
+            if spin_symmetry == Trivial
+                ε = zeros(ComplexF64, 3, 3, 3)
+                for i in 1:3
+                    ε[mod1(i, 3), mod1(i + 1, 3), mod1(i + 2, 3)] = 1
+                    ε[mod1(i, 3), mod1(i - 1, 3), mod1(i - 2, 3)] = -1
+                end
+                Svec = [tJ.S_x(particle_symmetry, spin_symmetry; slave_fermion),
+                        tJ.S_y(particle_symmetry, spin_symmetry; slave_fermion),
+                        tJ.S_z(particle_symmetry, spin_symmetry; slave_fermion)]
+                # Hermiticity
+                for s in Svec
+                    @test s' ≈ s
+                end
+                # operators should be normalized
+                S = 1 / 2
+                @test sum(tr(Svec[i]^2) for i in 1:3) / (2S + 1) ≈ S * (S + 1)
+                # test S_plus and S_min
+                @test tJ.S_plusmin(particle_symmetry, spin_symmetry; slave_fermion) ≈
+                      tJ.S_plus(particle_symmetry, spin_symmetry; slave_fermion) ⊗
+                      tJ.S_min(particle_symmetry, spin_symmetry; slave_fermion)
+                # commutation relations
+                for i in 1:3, j in 1:3
+                    @test Svec[i] * Svec[j] - Svec[j] * Svec[i] ≈
+                          sum(im * ε[i, j, k] * Svec[k] for k in 1:3)
+                end
+            end
+        else
+            @test_broken tJ.e_plusmin(particle_symmetry, spin_symmetry; slave_fermion)
+            @test_broken tJ.e_minplus(particle_symmetry, spin_symmetry; slave_fermion)
+        end
+    end
+end
+
+function hamiltonian(particle_symmetry, spin_symmetry; t, J, mu, L, slave_fermion)
+    num = tJ.e_number(particle_symmetry, spin_symmetry; slave_fermion)
+    hop_heis = (-t) * (tJ.e_plusmin(particle_symmetry, spin_symmetry; slave_fermion) +
+                       tJ.e_minplus(particle_symmetry, spin_symmetry; slave_fermion)) +
+               J *
+               (tJ.S_exchange(particle_symmetry, spin_symmetry; slave_fermion) -
+                (1 / 4) * (num ⊗ num))
+    chemical_potential = (-mu) * num
+    I = id(tJ.tj_space(particle_symmetry, spin_symmetry; slave_fermion))
+    H = sum(1:(L - 1)) do i
+        return reduce(⊗, insert!(collect(Any, fill(I, L - 2)), i, hop_heis))
+    end + sum(1:L) do i
+          return reduce(⊗, insert!(collect(Any, fill(I, L - 1)), i, chemical_potential))
+          end
+    return H
+end
+
+@testset "spectrum" begin
+    L = 4
+    t = randn()
+    J = randn()
+    mu = randn()
+
+    for slave_fermion in (false, true)
+        H_triv = hamiltonian(Trivial, Trivial; t, J, mu, L, slave_fermion)
+        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)
+                continue
+            end
+            H_symm = hamiltonian(particle_symmetry, spin_symmetry; t, J, mu, L,
+                                 slave_fermion)
+            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
+    end
+end