diff --git a/docs/src/index.md b/docs/src/index.md
index 1dd2420..e536830 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -10,5 +10,5 @@ Documentation for [TensorKitTensors](https://github.com/QuantumKitHub/TensorKitT
 ```
 
 ```@autodocs
-Modules = [TensorKitTensors, TensorKitTensors.SpinOperators, TensorKitTensors.BosonOperators, TensorKitTensors.HubbardOperators, TensorKitTensors.TJOperators]
+Modules = [TensorKitTensors, TensorKitTensors.SpinOperators, TensorKitTensors.BosonOperators, TensorKitTensors.HubbardOperators, TensorKitTensors.TJOperators, TensorKitTensors.FermionOperators]
 ```
diff --git a/src/TensorKitTensors.jl b/src/TensorKitTensors.jl
index a0ca6d1..60c1744 100644
--- a/src/TensorKitTensors.jl
+++ b/src/TensorKitTensors.jl
@@ -4,10 +4,12 @@ export SpinOperators
 export BosonOperators
 export HubbardOperators
 export TJOperators
+export FermionOperators
 
 include("spinoperators.jl")
 include("bosonoperators.jl")
 include("hubbardoperators.jl")
 include("tjoperators.jl")
+include("fermionoperators.jl")
 
 end
diff --git a/src/fermionoperators.jl b/src/fermionoperators.jl
new file mode 100644
index 0000000..707fef6
--- /dev/null
+++ b/src/fermionoperators.jl
@@ -0,0 +1,107 @@
+#===========================================================================================
+    Spinless fermions
+===========================================================================================#
+
+module FermionOperators
+
+using TensorKit
+
+export fermion_space
+export c_num
+export c_plus_c_min, c_min_c_plus, c_plus_c_plus, c_min_c_min
+export n
+export c⁺c⁻, c⁻c⁺, c⁺c⁺, c⁻c⁻
+
+"""
+    fermion_space()
+
+Return the local hilbert space for a model of spinless fermions.
+"""
+function fermion_space()
+    return Vect[fℤ₂](0 => 1, 1 => 1)
+end
+
+# Single-site operators
+# ---------------------
+function single_site_operator(T)
+    V = fermion_space()
+    return zeros(T, V ← V)
+end
+
+@doc """
+    c_num([elt::Type{<:Number}=ComplexF64])
+    n([elt::Type{<:Number}=ComplexF64])
+
+Return the one-body operator that counts the nunber of particles.
+""" c_num
+function c_num(T::Type{<:Number}=ComplexF64)
+    t = single_site_operator(T)
+    block(t, fℤ₂(1)) .= one(T)
+    return t
+end
+const n = c_num
+
+# Two site operators
+# ------------------
+function two_site_operator(T::Type{<:Number}=ComplexF64)
+    V = fermion_space()
+    return zeros(T, V ⊗ V ← V ⊗ V)
+end
+
+@doc """
+    c_plus_c_min([elt::Type{<:Number}=ComplexF64])
+    c⁺c⁻([elt::Type{<:Number}=ComplexF64])
+
+Return the two-body operator that creates a particle at the first site and annihilates a particle at the second.
+""" c_plus_c_min
+function c_plus_c_min(T::Type{<:Number}=ComplexF64)
+    t = two_site_operator(T)
+    I = sectortype(t)
+    t[(I(1), I(0), dual(I(0)), dual(I(1)))] .= 1
+    return t
+end
+const c⁺c⁻ = c_plus_c_min
+
+@doc """
+    c_min_c_plus([elt::Type{<:Number}=ComplexF64])
+    c⁻c⁺([elt::Type{<:Number}=ComplexF64])
+
+Return the two-body operator that annihilates a particle at the first site and creates a particle at the second.
+""" c_min_c_plus
+function c_min_c_plus(T::Type{<:Number}=ComplexF64)
+    t = two_site_operator(T)
+    I = sectortype(t)
+    t[(I(0), I(1), dual(I(1)), dual(I(0)))] .= 1
+    return t
+end
+const c⁻c⁺ = c_min_c_plus
+
+@doc """
+    c_plus_c_plus([elt::Type{<:Number}=ComplexF64])
+    c⁺c⁺([elt::Type{<:Number}=ComplexF64])
+
+Return the two-body operator that creates a particle at the first and at the second site.
+""" c_plus_c_plus
+function c_plus_c_plus(T::Type{<:Number}=ComplexF64)
+    t = two_site_operator(T)
+    I = sectortype(t)
+    t[(I(1), I(1), dual(I(0)), dual(I(0)))] .= 1
+    return t
+end
+const c⁺c⁺ = c_plus_c_plus
+
+@doc """
+    c_min_c_min([elt::Type{<:Number}=ComplexF64])
+    c⁻c⁻([elt::Type{<:Number}=ComplexF64])
+
+Return the two-body operator that annihilates a particle at the first and at the second site.
+""" c_min_c_min
+function c_min_c_min(T::Type{<:Number}=ComplexF64)
+    t = two_site_operator(T)
+    I = sectortype(t)
+    t[(I(0), I(0), dual(I(1)), dual(I(1)))] .= 1
+    return t
+end
+const c⁻c⁻ = c_min_c_min
+
+end
diff --git a/test/fermionoperators.jl b/test/fermionoperators.jl
new file mode 100644
index 0000000..840762b
--- /dev/null
+++ b/test/fermionoperators.jl
@@ -0,0 +1,45 @@
+using TensorKit
+using Test
+include("testsetup.jl")
+using TensorKitTensors
+using .TensorKitTensorsTestSetup
+using TensorKitTensors.FermionOperators
+using StableRNGs
+
+# anticommutation relations
+# {cᵢ†, cⱼ†} = 0 = {cᵢ, cⱼ}
+# {cᵢ, cⱼ†} = δᵢⱼ
+
+@testset "simple fermions" begin
+    @test c⁻c⁻() ≈ -permute(c⁻c⁻(), ((2, 1), (4, 3)))
+    @test c⁺c⁺() ≈ -permute(c⁺c⁺(), ((2, 1), (4, 3)))
+
+    # the following doesn't hold
+    # I don't think I can get all of these to hold simultaneously?
+    # @test cc⁺ ≈ -permute(c⁺c, (2, 1), (4, 3))
+
+    @test c⁻c⁺()' ≈ c⁺c⁻()
+    @test c⁻c⁻()' ≈ c⁺c⁺()
+    @test (c⁺c⁻() + c⁻c⁺())' ≈ c⁻c⁺() + c⁺c⁻()
+    @test (c⁺c⁻() - c⁻c⁺())' ≈ c⁻c⁺() - c⁺c⁻()
+
+    @plansor c_number[-1; -2] := c⁺c⁻()[-1 1; 3 2] * τ[3 2; -2 1]
+    @test c_number ≈ c_num()
+end
+
+@testset "Exact Diagonalization" begin
+    rng = StableRNG(123)
+
+    L = 2
+    t, V, mu = rand(rng, 3)
+    pspace = Vect[fℤ₂](0 => 1, 1 => 1)
+
+    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])
+
+    eigenvals = get_lowest_eigenvalues(H, -1; L)
+    @test eigenvals ≈ true_eigenvalues
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 0fdbff3..06d4729 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,5 +1,13 @@
 using SafeTestsets
 
+@time @safetestset "spin operators" begin
+    include("spinoperators.jl")
+end
+
+@time @safetestset "boson operators" begin
+    include("bosonoperators.jl")
+end
+
 @time @safetestset "Hubbard operators" begin
     include("hubbardoperators.jl")
 end
@@ -8,12 +16,8 @@ end
     include("tjoperators.jl")
 end
 
-@time @safetestset "spin operators" begin
-    include("spinoperators.jl")
-end
-
-@time @safetestset "boson operators" begin
-    include("bosonoperators.jl")
+@time @safetestset "spinless fermion operators" begin
+    include("fermionoperators.jl")
 end
 
 @time @safetestset "Aqua" begin
diff --git a/test/testsetup.jl b/test/testsetup.jl
index 3f93623..b1799ce 100644
--- a/test/testsetup.jl
+++ b/test/testsetup.jl
@@ -1,6 +1,6 @@
 module TensorKitTensorsTestSetup
 
-export test_operator, operator_sum
+export test_operator, operator_sum, get_lowest_eigenvalues
 
 using Test
 using TensorKit
@@ -29,4 +29,16 @@ function test_operator(O1::AbstractTensorMap, O2::AbstractTensorMap; L::Int=4,
                    isapproxkwargs...)
 end
 
+function get_lowest_eigenvalues(O1::AbstractTensorMap, n::Int; L::Int=4)
+    H = operator_sum(O1; L)
+    eigenvals = mapreduce(vcat, eigvals(H)) do (c, vals)
+        return repeat(vals, dim(c))
+    end
+    sort!(eigenvals; by=real)
+    if n == -1
+        return eigenvals
+    end
+    return eigenvals[1:n]
+end
+
 end