diff --git a/tests/python/operators/fermi_hubbard_test.py b/tests/python/operators/fermi_hubbard_test.py
index e4457b6f7..42e3a8d3b 100644
--- a/tests/python/operators/fermi_hubbard_test.py
+++ b/tests/python/operators/fermi_hubbard_test.py
@@ -197,6 +197,193 @@ def test_fermi_hubbard_2d():
 
     # periodic boundary conditions
     op_periodic = fermi_hubbard_2d(
+        norb_x=3,
+        norb_y=3,
+        tunneling=1,
+        interaction=2,
+        chemical_potential=3,
+        nearest_neighbor_interaction=4,
+        periodic=True,
+    )
+    np.testing.assert_equal(
+        dict(op_periodic),
+        {
+            (cre_a(0), des_a(3)): -1,
+            (cre_b(0), des_b(3)): -1,
+            (cre_a(3), des_a(0)): -1,
+            (cre_b(3), des_b(0)): -1,
+            (cre_a(0), des_a(1)): -1,
+            (cre_b(0), des_b(1)): -1,
+            (cre_a(1), des_a(0)): -1,
+            (cre_b(1), des_b(0)): -1,
+            (cre_a(1), des_a(4)): -1,
+            (cre_b(1), des_b(4)): -1,
+            (cre_a(4), des_a(1)): -1,
+            (cre_b(4), des_b(1)): -1,
+            (cre_a(1), des_a(2)): -1,
+            (cre_b(1), des_b(2)): -1,
+            (cre_a(2), des_a(1)): -1,
+            (cre_b(2), des_b(1)): -1,
+            (cre_a(2), des_a(5)): -1,
+            (cre_b(2), des_b(5)): -1,
+            (cre_a(5), des_a(2)): -1,
+            (cre_b(5), des_b(2)): -1,
+            (cre_a(2), des_a(0)): -1,
+            (cre_b(2), des_b(0)): -1,
+            (cre_a(0), des_a(2)): -1,
+            (cre_b(0), des_b(2)): -1,
+            (cre_a(3), des_a(6)): -1,
+            (cre_b(3), des_b(6)): -1,
+            (cre_a(6), des_a(3)): -1,
+            (cre_b(6), des_b(3)): -1,
+            (cre_a(3), des_a(4)): -1,
+            (cre_b(3), des_b(4)): -1,
+            (cre_a(4), des_a(3)): -1,
+            (cre_b(4), des_b(3)): -1,
+            (cre_a(4), des_a(7)): -1,
+            (cre_b(4), des_b(7)): -1,
+            (cre_a(7), des_a(4)): -1,
+            (cre_b(7), des_b(4)): -1,
+            (cre_a(4), des_a(5)): -1,
+            (cre_b(4), des_b(5)): -1,
+            (cre_a(5), des_a(4)): -1,
+            (cre_b(5), des_b(4)): -1,
+            (cre_a(5), des_a(8)): -1,
+            (cre_b(5), des_b(8)): -1,
+            (cre_a(8), des_a(5)): -1,
+            (cre_b(8), des_b(5)): -1,
+            (cre_a(5), des_a(3)): -1,
+            (cre_b(5), des_b(3)): -1,
+            (cre_a(3), des_a(5)): -1,
+            (cre_b(3), des_b(5)): -1,
+            (cre_a(6), des_a(0)): -1,
+            (cre_b(6), des_b(0)): -1,
+            (cre_a(0), des_a(6)): -1,
+            (cre_b(0), des_b(6)): -1,
+            (cre_a(6), des_a(7)): -1,
+            (cre_b(6), des_b(7)): -1,
+            (cre_a(7), des_a(6)): -1,
+            (cre_b(7), des_b(6)): -1,
+            (cre_a(7), des_a(1)): -1,
+            (cre_b(7), des_b(1)): -1,
+            (cre_a(1), des_a(7)): -1,
+            (cre_b(1), des_b(7)): -1,
+            (cre_a(7), des_a(8)): -1,
+            (cre_b(7), des_b(8)): -1,
+            (cre_a(8), des_a(7)): -1,
+            (cre_b(8), des_b(7)): -1,
+            (cre_a(8), des_a(2)): -1,
+            (cre_b(8), des_b(2)): -1,
+            (cre_a(2), des_a(8)): -1,
+            (cre_b(2), des_b(8)): -1,
+            (cre_a(8), des_a(6)): -1,
+            (cre_b(8), des_b(6)): -1,
+            (cre_a(6), des_a(8)): -1,
+            (cre_b(6), des_b(8)): -1,
+            (cre_a(0), des_a(0), cre_b(0), des_b(0)): 2,
+            (cre_a(1), des_a(1), cre_b(1), des_b(1)): 2,
+            (cre_a(2), des_a(2), cre_b(2), des_b(2)): 2,
+            (cre_a(3), des_a(3), cre_b(3), des_b(3)): 2,
+            (cre_a(4), des_a(4), cre_b(4), des_b(4)): 2,
+            (cre_a(5), des_a(5), cre_b(5), des_b(5)): 2,
+            (cre_a(6), des_a(6), cre_b(6), des_b(6)): 2,
+            (cre_a(7), des_a(7), cre_b(7), des_b(7)): 2,
+            (cre_a(8), des_a(8), cre_b(8), des_b(8)): 2,
+            (cre_a(0), des_a(0)): -3,
+            (cre_b(0), des_b(0)): -3,
+            (cre_a(1), des_a(1)): -3,
+            (cre_b(1), des_b(1)): -3,
+            (cre_a(2), des_a(2)): -3,
+            (cre_b(2), des_b(2)): -3,
+            (cre_a(3), des_a(3)): -3,
+            (cre_b(3), des_b(3)): -3,
+            (cre_a(4), des_a(4)): -3,
+            (cre_b(4), des_b(4)): -3,
+            (cre_a(5), des_a(5)): -3,
+            (cre_b(5), des_b(5)): -3,
+            (cre_a(6), des_a(6)): -3,
+            (cre_b(6), des_b(6)): -3,
+            (cre_a(7), des_a(7)): -3,
+            (cre_b(7), des_b(7)): -3,
+            (cre_a(8), des_a(8)): -3,
+            (cre_b(8), des_b(8)): -3,
+            (cre_a(0), des_a(0), cre_a(3), des_a(3)): 4,
+            (cre_a(0), des_a(0), cre_b(3), des_b(3)): 4,
+            (cre_b(0), des_b(0), cre_a(3), des_a(3)): 4,
+            (cre_b(0), des_b(0), cre_b(3), des_b(3)): 4,
+            (cre_a(0), des_a(0), cre_a(1), des_a(1)): 4,
+            (cre_a(0), des_a(0), cre_b(1), des_b(1)): 4,
+            (cre_b(0), des_b(0), cre_a(1), des_a(1)): 4,
+            (cre_b(0), des_b(0), cre_b(1), des_b(1)): 4,
+            (cre_a(1), des_a(1), cre_a(4), des_a(4)): 4,
+            (cre_a(1), des_a(1), cre_b(4), des_b(4)): 4,
+            (cre_b(1), des_b(1), cre_a(4), des_a(4)): 4,
+            (cre_b(1), des_b(1), cre_b(4), des_b(4)): 4,
+            (cre_a(1), des_a(1), cre_a(2), des_a(2)): 4,
+            (cre_a(1), des_a(1), cre_b(2), des_b(2)): 4,
+            (cre_b(1), des_b(1), cre_a(2), des_a(2)): 4,
+            (cre_b(1), des_b(1), cre_b(2), des_b(2)): 4,
+            (cre_a(2), des_a(2), cre_a(5), des_a(5)): 4,
+            (cre_a(2), des_a(2), cre_b(5), des_b(5)): 4,
+            (cre_b(2), des_b(2), cre_a(5), des_a(5)): 4,
+            (cre_b(2), des_b(2), cre_b(5), des_b(5)): 4,
+            (cre_a(2), des_a(2), cre_a(0), des_a(0)): 4,
+            (cre_a(2), des_a(2), cre_b(0), des_b(0)): 4,
+            (cre_b(2), des_b(2), cre_a(0), des_a(0)): 4,
+            (cre_b(2), des_b(2), cre_b(0), des_b(0)): 4,
+            (cre_a(3), des_a(3), cre_a(6), des_a(6)): 4,
+            (cre_a(3), des_a(3), cre_b(6), des_b(6)): 4,
+            (cre_b(3), des_b(3), cre_a(6), des_a(6)): 4,
+            (cre_b(3), des_b(3), cre_b(6), des_b(6)): 4,
+            (cre_a(3), des_a(3), cre_a(4), des_a(4)): 4,
+            (cre_a(3), des_a(3), cre_b(4), des_b(4)): 4,
+            (cre_b(3), des_b(3), cre_a(4), des_a(4)): 4,
+            (cre_b(3), des_b(3), cre_b(4), des_b(4)): 4,
+            (cre_a(4), des_a(4), cre_a(7), des_a(7)): 4,
+            (cre_a(4), des_a(4), cre_b(7), des_b(7)): 4,
+            (cre_b(4), des_b(4), cre_a(7), des_a(7)): 4,
+            (cre_b(4), des_b(4), cre_b(7), des_b(7)): 4,
+            (cre_a(4), des_a(4), cre_a(5), des_a(5)): 4,
+            (cre_a(4), des_a(4), cre_b(5), des_b(5)): 4,
+            (cre_b(4), des_b(4), cre_a(5), des_a(5)): 4,
+            (cre_b(4), des_b(4), cre_b(5), des_b(5)): 4,
+            (cre_a(5), des_a(5), cre_a(8), des_a(8)): 4,
+            (cre_a(5), des_a(5), cre_b(8), des_b(8)): 4,
+            (cre_b(5), des_b(5), cre_a(8), des_a(8)): 4,
+            (cre_b(5), des_b(5), cre_b(8), des_b(8)): 4,
+            (cre_a(5), des_a(5), cre_a(3), des_a(3)): 4,
+            (cre_a(5), des_a(5), cre_b(3), des_b(3)): 4,
+            (cre_b(5), des_b(5), cre_a(3), des_a(3)): 4,
+            (cre_b(5), des_b(5), cre_b(3), des_b(3)): 4,
+            (cre_a(6), des_a(6), cre_a(0), des_a(0)): 4,
+            (cre_a(6), des_a(6), cre_b(0), des_b(0)): 4,
+            (cre_b(6), des_b(6), cre_a(0), des_a(0)): 4,
+            (cre_b(6), des_b(6), cre_b(0), des_b(0)): 4,
+            (cre_a(6), des_a(6), cre_a(7), des_a(7)): 4,
+            (cre_a(6), des_a(6), cre_b(7), des_b(7)): 4,
+            (cre_b(6), des_b(6), cre_a(7), des_a(7)): 4,
+            (cre_b(6), des_b(6), cre_b(7), des_b(7)): 4,
+            (cre_a(7), des_a(7), cre_a(1), des_a(1)): 4,
+            (cre_a(7), des_a(7), cre_b(1), des_b(1)): 4,
+            (cre_b(7), des_b(7), cre_a(1), des_a(1)): 4,
+            (cre_b(7), des_b(7), cre_b(1), des_b(1)): 4,
+            (cre_a(7), des_a(7), cre_a(8), des_a(8)): 4,
+            (cre_a(7), des_a(7), cre_b(8), des_b(8)): 4,
+            (cre_b(7), des_b(7), cre_a(8), des_a(8)): 4,
+            (cre_b(7), des_b(7), cre_b(8), des_b(8)): 4,
+            (cre_a(8), des_a(8), cre_a(2), des_a(2)): 4,
+            (cre_a(8), des_a(8), cre_b(2), des_b(2)): 4,
+            (cre_b(8), des_b(8), cre_a(2), des_a(2)): 4,
+            (cre_b(8), des_b(8), cre_b(2), des_b(2)): 4,
+            (cre_a(8), des_a(8), cre_a(6), des_a(6)): 4,
+            (cre_a(8), des_a(8), cre_b(6), des_b(6)): 4,
+            (cre_b(8), des_b(8), cre_a(6), des_a(6)): 4,
+            (cre_b(8), des_b(8), cre_b(6), des_b(6)): 4,
+        },
+    )
+
+    # periodic boundary conditions (edge case)
+    op_periodic_edge = fermi_hubbard_2d(
         norb_x=2,
         norb_y=2,
         tunneling=1,
@@ -206,7 +393,7 @@ def test_fermi_hubbard_2d():
         periodic=True,
     )
     np.testing.assert_equal(
-        dict(op_periodic),
+        dict(op_periodic_edge),
         {
             (cre_a(0), des_a(2)): -2,
             (cre_b(0), des_b(2)): -2,
@@ -299,11 +486,11 @@ def test_non_interacting_fermi_hubbard_2d_eigenvalue():
     eigs, _ = scipy.sparse.linalg.eigsh(ham, which="SA", k=1)
     np.testing.assert_allclose(eigs[0], -4.000000000000)
 
-    # periodic boundary conditions
-    op_periodic = fermi_hubbard_2d(
+    # periodic boundary conditions (edge case)
+    op_periodic_edge = fermi_hubbard_2d(
         norb_x=2, norb_y=2, tunneling=1, interaction=0, periodic=True
     )
-    ham_periodic = ffsim.linear_operator(op_periodic, norb=4, nelec=(2, 2))
+    ham_periodic = ffsim.linear_operator(op_periodic_edge, norb=4, nelec=(2, 2))
     eigs_periodic, _ = scipy.sparse.linalg.eigsh(ham_periodic, which="SA", k=1)
     print(eigs_periodic)
     np.testing.assert_allclose(eigs_periodic[0], -8.000000000000)
@@ -336,11 +523,11 @@ def test_fermi_hubbard_2d_with_interaction_eigenvalue():
     eigs, _ = scipy.sparse.linalg.eigsh(ham, which="SA", k=1)
     np.testing.assert_allclose(eigs[0], -2.828427124746)
 
-    # periodic boundary conditions
-    op_periodic = fermi_hubbard_2d(
+    # periodic boundary conditions (edge case)
+    op_periodic_edge = fermi_hubbard_2d(
         norb_x=2, norb_y=2, tunneling=1, interaction=2, periodic=True
     )
-    ham_periodic = ffsim.linear_operator(op_periodic, norb=4, nelec=(2, 2))
+    ham_periodic = ffsim.linear_operator(op_periodic_edge, norb=4, nelec=(2, 2))
     eigs_periodic, _ = scipy.sparse.linalg.eigsh(ham_periodic, which="SA", k=1)
     np.testing.assert_allclose(eigs_periodic[0], -6.681695234497)
 
@@ -376,8 +563,8 @@ def test_fermi_hubbard_2d_with_chemical_potential_eigenvalue():
     eigs, _ = scipy.sparse.linalg.eigsh(ham, which="SA", k=1)
     np.testing.assert_allclose(eigs[0], -14.828427124746)
 
-    # periodic boundary conditions
-    op_periodic = fermi_hubbard_2d(
+    # periodic boundary conditions (edge case)
+    op_periodic_edge = fermi_hubbard_2d(
         norb_x=2,
         norb_y=2,
         tunneling=1,
@@ -385,7 +572,7 @@ def test_fermi_hubbard_2d_with_chemical_potential_eigenvalue():
         chemical_potential=3,
         periodic=True,
     )
-    ham_periodic = ffsim.linear_operator(op_periodic, norb=4, nelec=(2, 2))
+    ham_periodic = ffsim.linear_operator(op_periodic_edge, norb=4, nelec=(2, 2))
     eigs_periodic, _ = scipy.sparse.linalg.eigsh(ham_periodic, which="SA", k=1)
     np.testing.assert_allclose(eigs_periodic[0], -18.681695234497)
 
@@ -439,8 +626,8 @@ def test_fermi_hubbard_2d_with_nearest_neighbor_interaction_eigenvalue():
     eigs, _ = scipy.sparse.linalg.eigsh(ham, which="SA", k=1)
     np.testing.assert_allclose(eigs[0], -8.781962448006)
 
-    # periodic boundary conditions
-    op_periodic = fermi_hubbard_2d(
+    # periodic boundary conditions (edge case)
+    op_periodic_edge = fermi_hubbard_2d(
         norb_x=2,
         norb_y=2,
         tunneling=1,
@@ -449,7 +636,7 @@ def test_fermi_hubbard_2d_with_nearest_neighbor_interaction_eigenvalue():
         nearest_neighbor_interaction=4,
         periodic=True,
     )
-    ham_periodic = ffsim.linear_operator(op_periodic, norb=4, nelec=(2, 2))
+    ham_periodic = ffsim.linear_operator(op_periodic_edge, norb=4, nelec=(2, 2))
     eigs_periodic, _ = scipy.sparse.linalg.eigsh(ham_periodic, which="SA", k=1)
     np.testing.assert_allclose(eigs_periodic[0], -9.428197577536)
 
@@ -501,8 +688,8 @@ def test_fermi_hubbard_2d_with_unequal_filling_eigenvalue():
     eigs, _ = scipy.sparse.linalg.eigsh(ham, which="SA", k=1)
     np.testing.assert_allclose(eigs[0], -0.828427124746)
 
-    # periodic boundary conditions
-    op_periodic = fermi_hubbard_2d(
+    # periodic boundary conditions (edge case)
+    op_periodic_edge = fermi_hubbard_2d(
         norb_x=2,
         norb_y=2,
         tunneling=1,
@@ -511,7 +698,7 @@ def test_fermi_hubbard_2d_with_unequal_filling_eigenvalue():
         nearest_neighbor_interaction=4,
         periodic=True,
     )
-    ham_periodic = ffsim.linear_operator(op_periodic, norb=4, nelec=(1, 3))
+    ham_periodic = ffsim.linear_operator(op_periodic_edge, norb=4, nelec=(1, 3))
     eigs_periodic, _ = scipy.sparse.linalg.eigsh(ham_periodic, which="SA", k=1)
     np.testing.assert_allclose(eigs_periodic[0], 8.743352161722)