diff --git a/python/ffsim/__init__.py b/python/ffsim/__init__.py
index cd02318a5..8c916e619 100644
--- a/python/ffsim/__init__.py
+++ b/python/ffsim/__init__.py
@@ -41,6 +41,7 @@
     des_a,
     des_b,
     fermi_hubbard_1d,
+    fermi_hubbard_2d,
     number_operator,
 )
 from ffsim.protocols import (
@@ -126,6 +127,7 @@
     "expectation_one_body_power",
     "expectation_one_body_product",
     "fermi_hubbard_1d",
+    "fermi_hubbard_2d",
     "fermion_operator",
     "hartree_fock_state",
     "indices_to_strings",
diff --git a/python/ffsim/operators/__init__.py b/python/ffsim/operators/__init__.py
index 4f1beddd1..29758742b 100644
--- a/python/ffsim/operators/__init__.py
+++ b/python/ffsim/operators/__init__.py
@@ -11,7 +11,7 @@
 """Operators."""
 
 from ffsim.operators.common_operators import number_operator
-from ffsim.operators.fermi_hubbard import fermi_hubbard_1d
+from ffsim.operators.fermi_hubbard import fermi_hubbard_1d, fermi_hubbard_2d
 from ffsim.operators.fermion_action import (
     FermionAction,
     cre,
@@ -33,5 +33,6 @@
     "des_a",
     "des_b",
     "fermi_hubbard_1d",
+    "fermi_hubbard_2d",
     "number_operator",
 ]
diff --git a/python/ffsim/operators/fermi_hubbard.py b/python/ffsim/operators/fermi_hubbard.py
index d185973e2..02e99b647 100644
--- a/python/ffsim/operators/fermi_hubbard.py
+++ b/python/ffsim/operators/fermi_hubbard.py
@@ -87,3 +87,107 @@ def fermi_hubbard_1d(
             coeffs[(cre_b(p), des_b(p))] = -chemical_potential
 
     return FermionOperator(coeffs)
+
+
+def fermi_hubbard_2d(
+    norb_x: int,
+    norb_y: int,
+    tunneling: float,
+    interaction: float,
+    *,
+    chemical_potential: float = 0,
+    nearest_neighbor_interaction: float = 0,
+    periodic: bool = False,
+) -> FermionOperator:
+    r"""Two-dimensional Fermi-Hubbard model Hamiltonian.
+
+    The Hamiltonian for the two-dimensional Fermi-Hubbard model on a square lattice with
+    :math:`N = N_x \times N_y` spatial orbitals is given by
+
+    .. math::
+
+        H = -t \sum_{\sigma} \sum_{\braket{pq}}
+        (a^\dagger_{\sigma, p} a_{\sigma, q} + \text{h.c.})
+        + U \sum_{p=1}^{N} n_{\alpha, p} n_{\beta, p}
+        - \mu \sum_{p=1}^{N} (n_{\alpha, p} + n_{\beta, p})
+        + V \sum_{\sigma, \sigma'} \sum_{\braket{pq}} n_{\sigma, p} n_{\sigma', q}
+
+    where :math:`\braket{\dots}` denotes nearest-neighbor pairs and
+    :math:`n_{\sigma, p} = a_{\sigma, p}^\dagger a_{\sigma, p}` is the number operator
+    on orbital :math:`p` with spin :math:`\sigma`.
+
+    References:
+        - `The Hubbard Model`_
+
+    Args:
+        norb_x: The number of spatial orbitals in the x-direction :math:`N_x`.
+        norb_y: The number of spatial orbitals in the y-direction :math:`N_y`.
+        tunneling: The tunneling amplitude :math:`t`.
+        interaction: The onsite interaction strength :math:`U`.
+        chemical_potential: The chemical potential :math:`\mu`.
+        nearest_neighbor_interaction: The nearest-neighbor interaction strength
+            :math:`V`.
+        periodic: Whether to use periodic boundary conditions.
+
+    Returns:
+        The two-dimensional Fermi-Hubbard model Hamiltonian.
+
+    .. _The Hubbard Model: https://doi.org/10.1146/annurev-conmatphys-031620-102024
+    """
+    coeffs: dict[tuple[tuple[bool, bool, int], ...], complex] = defaultdict(float)
+
+    for orb in range(norb_x * norb_y):
+        # map from row-major ordering to Cartesian coordinates
+        y, x = divmod(orb, norb_x)
+
+        # get the coordinates/orbitals to the right and up
+        x_right = (x + 1) % norb_x
+        y_up = (y + 1) % norb_y
+        orb_right = x_right + norb_x * y
+        orb_up = x + norb_x * y_up
+
+        if x != norb_x - 1 or periodic:
+            coeffs[(cre_a(orb), des_a(orb_right))] -= tunneling
+            coeffs[(cre_a(orb_right), des_a(orb))] -= tunneling
+            coeffs[(cre_b(orb), des_b(orb_right))] -= tunneling
+            coeffs[(cre_b(orb_right), des_b(orb))] -= tunneling
+            if nearest_neighbor_interaction:
+                coeffs[(cre_a(orb), des_a(orb), cre_a(orb_right), des_a(orb_right))] = (
+                    nearest_neighbor_interaction
+                )
+                coeffs[(cre_a(orb), des_a(orb), cre_b(orb_right), des_b(orb_right))] = (
+                    nearest_neighbor_interaction
+                )
+                coeffs[(cre_b(orb), des_b(orb), cre_a(orb_right), des_a(orb_right))] = (
+                    nearest_neighbor_interaction
+                )
+                coeffs[(cre_b(orb), des_b(orb), cre_b(orb_right), des_b(orb_right))] = (
+                    nearest_neighbor_interaction
+                )
+
+        if y != norb_y - 1 or periodic:
+            coeffs[(cre_a(orb), des_a(orb_up))] -= tunneling
+            coeffs[(cre_a(orb_up), des_a(orb))] -= tunneling
+            coeffs[(cre_b(orb), des_b(orb_up))] -= tunneling
+            coeffs[(cre_b(orb_up), des_b(orb))] -= tunneling
+            if nearest_neighbor_interaction:
+                coeffs[(cre_a(orb), des_a(orb), cre_a(orb_up), des_a(orb_up))] = (
+                    nearest_neighbor_interaction
+                )
+                coeffs[(cre_a(orb), des_a(orb), cre_b(orb_up), des_b(orb_up))] = (
+                    nearest_neighbor_interaction
+                )
+                coeffs[(cre_b(orb), des_b(orb), cre_a(orb_up), des_a(orb_up))] = (
+                    nearest_neighbor_interaction
+                )
+                coeffs[(cre_b(orb), des_b(orb), cre_b(orb_up), des_b(orb_up))] = (
+                    nearest_neighbor_interaction
+                )
+
+        if interaction:
+            coeffs[(cre_a(orb), des_a(orb), cre_b(orb), des_b(orb))] = interaction
+        if chemical_potential:
+            coeffs[(cre_a(orb), des_a(orb))] = -chemical_potential
+            coeffs[(cre_b(orb), des_b(orb))] = -chemical_potential
+
+    return FermionOperator(coeffs)
diff --git a/tests/python/operators/fermi_hubbard_test.py b/tests/python/operators/fermi_hubbard_test.py
index 0b9c39b1e..69a74f6fd 100644
--- a/tests/python/operators/fermi_hubbard_test.py
+++ b/tests/python/operators/fermi_hubbard_test.py
@@ -16,7 +16,7 @@
 import scipy
 
 import ffsim
-from ffsim.operators.fermi_hubbard import fermi_hubbard_1d
+from ffsim.operators.fermi_hubbard import fermi_hubbard_1d, fermi_hubbard_2d
 from ffsim.operators.fermion_action import cre_a, cre_b, des_a, des_b
 
 
@@ -166,6 +166,332 @@ def test_fermi_hubbard_1d():
     )
 
 
+def test_fermi_hubbard_2d():
+    """Test terms of the two-dimensional Fermi-Hubbard model Hamiltonian."""
+
+    # open boundary conditions
+    op = fermi_hubbard_2d(
+        norb_x=2,
+        norb_y=2,
+        tunneling=1,
+        interaction=2,
+        chemical_potential=3,
+        nearest_neighbor_interaction=4,
+    )
+    np.testing.assert_equal(
+        dict(op),
+        {
+            (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(0), des_a(2)): -1,
+            (cre_b(0), des_b(2)): -1,
+            (cre_a(2), des_a(0)): -1,
+            (cre_b(2), des_b(0)): -1,
+            (cre_a(2), des_a(3)): -1,
+            (cre_b(2), des_b(3)): -1,
+            (cre_a(3), des_a(2)): -1,
+            (cre_b(3), des_b(2)): -1,
+            (cre_a(1), des_a(3)): -1,
+            (cre_b(1), des_b(3)): -1,
+            (cre_a(3), des_a(1)): -1,
+            (cre_b(3), des_b(1)): -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(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(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(0), des_a(0), cre_a(2), des_a(2)): 4,
+            (cre_a(0), des_a(0), cre_b(2), des_b(2)): 4,
+            (cre_b(0), des_b(0), cre_a(2), des_a(2)): 4,
+            (cre_b(0), des_b(0), cre_b(2), des_b(2)): 4,
+            (cre_a(2), des_a(2), cre_a(3), des_a(3)): 4,
+            (cre_a(2), des_a(2), cre_b(3), des_b(3)): 4,
+            (cre_b(2), des_b(2), cre_a(3), des_a(3)): 4,
+            (cre_b(2), des_b(2), cre_b(3), des_b(3)): 4,
+            (cre_a(1), des_a(1), cre_a(3), des_a(3)): 4,
+            (cre_a(1), des_a(1), cre_b(3), des_b(3)): 4,
+            (cre_b(1), des_b(1), cre_a(3), des_a(3)): 4,
+            (cre_b(1), des_b(1), cre_b(3), des_b(3)): 4,
+        },
+    )
+
+    # 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(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(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(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(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(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(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(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(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(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(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(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(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(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(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(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(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(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(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(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(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(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(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(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(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(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(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(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(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(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(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(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(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(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(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(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(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,
+            (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,
+        },
+    )
+
+    # periodic boundary conditions (edge case)
+    op_periodic_edge = fermi_hubbard_2d(
+        norb_x=2,
+        norb_y=2,
+        tunneling=1,
+        interaction=2,
+        chemical_potential=3,
+        nearest_neighbor_interaction=4,
+        periodic=True,
+    )
+    np.testing.assert_equal(
+        dict(op_periodic_edge),
+        {
+            (cre_a(0), des_a(1)): -2,
+            (cre_b(0), des_b(1)): -2,
+            (cre_a(1), des_a(0)): -2,
+            (cre_b(1), des_b(0)): -2,
+            (cre_a(0), des_a(2)): -2,
+            (cre_b(0), des_b(2)): -2,
+            (cre_a(2), des_a(0)): -2,
+            (cre_b(2), des_b(0)): -2,
+            (cre_a(2), des_a(3)): -2,
+            (cre_b(2), des_b(3)): -2,
+            (cre_a(3), des_a(2)): -2,
+            (cre_b(3), des_b(2)): -2,
+            (cre_a(1), des_a(3)): -2,
+            (cre_b(1), des_b(3)): -2,
+            (cre_a(3), des_a(1)): -2,
+            (cre_b(3), des_b(1)): -2,
+            (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(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(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(0), des_a(0), cre_a(2), des_a(2)): 4,
+            (cre_a(0), des_a(0), cre_b(2), des_b(2)): 4,
+            (cre_b(0), des_b(0), cre_a(2), des_a(2)): 4,
+            (cre_b(0), des_b(0), cre_b(2), des_b(2)): 4,
+            (cre_a(2), des_a(2), cre_a(3), des_a(3)): 4,
+            (cre_a(2), des_a(2), cre_b(3), des_b(3)): 4,
+            (cre_b(2), des_b(2), cre_a(3), des_a(3)): 4,
+            (cre_b(2), des_b(2), cre_b(3), des_b(3)): 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(1), des_a(1), cre_a(0), des_a(0)): 4,
+            (cre_a(1), des_a(1), cre_b(0), des_b(0)): 4,
+            (cre_b(1), des_b(1), cre_a(0), des_a(0)): 4,
+            (cre_b(1), des_b(1), cre_b(0), des_b(0)): 4,
+            (cre_a(1), des_a(1), cre_a(3), des_a(3)): 4,
+            (cre_a(1), des_a(1), cre_b(3), des_b(3)): 4,
+            (cre_b(1), des_b(1), cre_a(3), des_a(3)): 4,
+            (cre_b(1), des_b(1), cre_b(3), des_b(3)): 4,
+            (cre_a(3), des_a(3), cre_a(2), des_a(2)): 4,
+            (cre_a(3), des_a(3), cre_b(2), des_b(2)): 4,
+            (cre_b(3), des_b(3), cre_a(2), des_a(2)): 4,
+            (cre_b(3), des_b(3), cre_b(2), des_b(2)): 4,
+            (cre_a(3), des_a(3), cre_a(1), des_a(1)): 4,
+            (cre_a(3), des_a(3), cre_b(1), des_b(1)): 4,
+            (cre_b(3), des_b(3), cre_a(1), des_a(1)): 4,
+            (cre_b(3), des_b(3), cre_b(1), des_b(1)): 4,
+        },
+    )
+
+
 def test_non_interacting_fermi_hubbard_1d_eigenvalue():
     """Test ground-state eigenvalue of the non-interacting one-dimensional Fermi-Hubbard
     model Hamiltonian."""
@@ -191,6 +517,26 @@ def test_non_interacting_fermi_hubbard_1d_eigenvalue():
     np.testing.assert_allclose(eigs_periodic[0], -4.000000000000)
 
 
+def test_non_interacting_fermi_hubbard_2d_eigenvalue():
+    """Test ground-state eigenvalue of the non-interacting two-dimensional Fermi-Hubbard
+    model Hamiltonian."""
+
+    # open boundary conditions
+    op = fermi_hubbard_2d(norb_x=2, norb_y=2, tunneling=1, interaction=0)
+    ham = ffsim.linear_operator(op, norb=4, nelec=(2, 2))
+    eigs, _ = scipy.sparse.linalg.eigsh(ham, which="SA", k=1)
+    np.testing.assert_allclose(eigs[0], -4.000000000000)
+
+    # 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_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)
+
+
 def test_fermi_hubbard_1d_with_interaction_eigenvalue():
     """Test ground-state eigenvalue of the one-dimensional Fermi-Hubbard model
     Hamiltonian with onsite interaction."""
@@ -216,6 +562,25 @@ def test_fermi_hubbard_1d_with_interaction_eigenvalue():
     np.testing.assert_allclose(eigs_periodic[0], -3.123105625618)
 
 
+def test_fermi_hubbard_2d_with_interaction_eigenvalue():
+    """Test ground-state eigenvalue of the two-dimensional Fermi-Hubbard model
+    Hamiltonian with onsite interaction."""
+
+    # open boundary conditions
+    op = fermi_hubbard_2d(norb_x=2, norb_y=2, tunneling=1, interaction=2)
+    ham = ffsim.linear_operator(op, norb=4, nelec=(2, 2))
+    eigs, _ = scipy.sparse.linalg.eigsh(ham, which="SA", k=1)
+    np.testing.assert_allclose(eigs[0], -2.828427124746)
+
+    # 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_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)
+
+
 def test_fermi_hubbard_1d_with_chemical_potential_eigenvalue():
     """Test ground-state eigenvalue of the one-dimensional Fermi-Hubbard model
     Hamiltonian with onsite interaction and chemical potential."""
@@ -243,6 +608,32 @@ def test_fermi_hubbard_1d_with_chemical_potential_eigenvalue():
     np.testing.assert_allclose(eigs_periodic[0], -9.123105625618)
 
 
+def test_fermi_hubbard_2d_with_chemical_potential_eigenvalue():
+    """Test ground-state eigenvalue of the two-dimensional Fermi-Hubbard model
+    Hamiltonian with onsite interaction and chemical potential."""
+
+    # open boundary conditions
+    op = fermi_hubbard_2d(
+        norb_x=2, norb_y=2, tunneling=1, interaction=2, chemical_potential=3
+    )
+    ham = ffsim.linear_operator(op, norb=4, nelec=(2, 2))
+    eigs, _ = scipy.sparse.linalg.eigsh(ham, which="SA", k=1)
+    np.testing.assert_allclose(eigs[0], -14.828427124746)
+
+    # periodic boundary conditions (edge case)
+    op_periodic_edge = fermi_hubbard_2d(
+        norb_x=2,
+        norb_y=2,
+        tunneling=1,
+        interaction=2,
+        chemical_potential=3,
+        periodic=True,
+    )
+    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)
+
+
 def test_fermi_hubbard_1d_with_nearest_neighbor_interaction_eigenvalue():
     """Test ground-state eigenvalue of the one-dimensional Fermi-Hubbard model
     Hamiltonian with onsite interaction, chemical potential, and nearest-neighbor
@@ -287,6 +678,39 @@ def test_fermi_hubbard_1d_with_nearest_neighbor_interaction_eigenvalue():
     np.testing.assert_allclose(eigs_periodic[0], -6.000000000000)
 
 
+def test_fermi_hubbard_2d_with_nearest_neighbor_interaction_eigenvalue():
+    """Test ground-state eigenvalue of the two-dimensional Fermi-Hubbard model
+    Hamiltonian with onsite interaction, chemical potential, and nearest-neighbor
+    interaction."""
+
+    # open boundary conditions
+    op = fermi_hubbard_2d(
+        norb_x=2,
+        norb_y=2,
+        tunneling=1,
+        interaction=2,
+        chemical_potential=3,
+        nearest_neighbor_interaction=4,
+    )
+    ham = ffsim.linear_operator(op, norb=4, nelec=(2, 2))
+    eigs, _ = scipy.sparse.linalg.eigsh(ham, which="SA", k=1)
+    np.testing.assert_allclose(eigs[0], -8.781962448006)
+
+    # periodic boundary conditions (edge case)
+    op_periodic_edge = fermi_hubbard_2d(
+        norb_x=2,
+        norb_y=2,
+        tunneling=1,
+        interaction=2,
+        chemical_potential=3,
+        nearest_neighbor_interaction=4,
+        periodic=True,
+    )
+    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)
+
+
 def test_fermi_hubbard_1d_with_unequal_filling_eigenvalue():
     """Test ground-state eigenvalue of the one-dimensional Fermi-Hubbard model
     Hamiltonian with unequal filling."""
@@ -317,6 +741,38 @@ def test_fermi_hubbard_1d_with_unequal_filling_eigenvalue():
     np.testing.assert_allclose(eigs_periodic[0], -0.828427124746)
 
 
+def test_fermi_hubbard_2d_with_unequal_filling_eigenvalue():
+    """Test ground-state eigenvalue of the two-dimensional Fermi-Hubbard model
+    Hamiltonian with unequal filling."""
+
+    # open boundary conditions
+    op = fermi_hubbard_2d(
+        norb_x=2,
+        norb_y=2,
+        tunneling=1,
+        interaction=2,
+        chemical_potential=3,
+        nearest_neighbor_interaction=4,
+    )
+    ham = ffsim.linear_operator(op, norb=4, nelec=(1, 3))
+    eigs, _ = scipy.sparse.linalg.eigsh(ham, which="SA", k=1)
+    np.testing.assert_allclose(eigs[0], -0.828427124746)
+
+    # periodic boundary conditions (edge case)
+    op_periodic_edge = fermi_hubbard_2d(
+        norb_x=2,
+        norb_y=2,
+        tunneling=1,
+        interaction=2,
+        chemical_potential=3,
+        nearest_neighbor_interaction=4,
+        periodic=True,
+    )
+    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)
+
+
 def test_fermi_hubbard_1d_hermiticity():
     """Test hermiticity of the one-dimensional Fermi-Hubbard model Hamiltonian."""
 
@@ -348,3 +804,40 @@ def test_fermi_hubbard_1d_hermiticity():
         op_periodic, norb=n_orbitals, nelec=n_electrons
     )
     np.testing.assert_allclose(ham_periodic @ np.eye(dim), ham_periodic.H @ np.eye(dim))
+
+
+def test_fermi_hubbard_2d_hermiticity():
+    """Test hermiticity of the two-dimensional Fermi-Hubbard model Hamiltonian."""
+
+    n_orbitals_x = 2
+    n_orbitals_y = 2
+    n_orbitals = n_orbitals_x * n_orbitals_y
+    n_electrons = (3, 1)
+    dim = ffsim.dim(n_orbitals, n_electrons)
+
+    # open boundary conditions
+    op = fermi_hubbard_2d(
+        norb_x=n_orbitals_x,
+        norb_y=n_orbitals_y,
+        tunneling=1.1,
+        interaction=1.2,
+        chemical_potential=1.3,
+        nearest_neighbor_interaction=1.4,
+    )
+    ham = ffsim.linear_operator(op, norb=n_orbitals, nelec=n_electrons)
+    np.testing.assert_allclose(ham @ np.eye(dim), ham.H @ np.eye(dim))
+
+    # periodic boundary conditions
+    op_periodic = fermi_hubbard_2d(
+        norb_x=n_orbitals_x,
+        norb_y=n_orbitals_y,
+        tunneling=1.1,
+        interaction=1.2,
+        chemical_potential=1.3,
+        nearest_neighbor_interaction=1.4,
+        periodic=True,
+    )
+    ham_periodic = ffsim.linear_operator(
+        op_periodic, norb=n_orbitals, nelec=n_electrons
+    )
+    np.testing.assert_allclose(ham_periodic @ np.eye(dim), ham_periodic.H @ np.eye(dim))