added Fermionic and Jordan-Wigner MPO encodings of 1D Fermi-Hubbard model

CHANGELOG.md

@@ -11,7 +11,8 @@ This project adheres to [Semantic Versioning], with the exception that minor rel
 
 ### Added
 
-- added deterministic ensemble evolution with optional autocorrelator and two-time correlator outputs, including periodic-wrap two-site observable support on `(L-1, 0)` ([#409]) ([**@Gauthameshwar**])\
+- added Fermionic and Jordan-Wigner MPO encodings of 1D Fermi-Hubbard model ([#220]) ([**@thilomueller**])
+- added deterministic ensemble evolution with optional autocorrelator and two-time correlator outputs, including periodic-wrap two-site observable support on `(L-1, 0)` ([#409]) ([**@Gauthameshwar**])
 
 ### Changed
 
@@ -100,7 +101,7 @@ _📚 Refer to the [GitHub Release Notes](https://github.com/munich-quantum-tool
 
 <!-- Version links -->
 
-[Unreleased]: https://github.com/munich-quantum-toolkit/yaqs/compare/v0.4.0...HEAD
+[Unreleased]: https://github.com/munich-quantum-toolkit/yaqs/compare/v0.5.0...HEAD
 [0.5.0]: https://github.com/munich-quantum-toolkit/yaqs/compare/v0.5.0
 [0.4.0]: https://github.com/munich-quantum-toolkit/yaqs/releases/tag/v0.4.0
 [0.3.3]: https://github.com/munich-quantum-toolkit/yaqs/releases/tag/v0.3.3
@@ -109,6 +110,7 @@ _📚 Refer to the [GitHub Release Notes](https://github.com/munich-quantum-tool
 
 <!-- PR links -->
 
+[#220]: https://github.com/munich-quantum-toolkit/yaqs/pull/220
 [#420]: https://github.com/munich-quantum-toolkit/yaqs/pull/420
 [#409]: https://github.com/munich-quantum-toolkit/yaqs/pull/409
 [#344]: https://github.com/munich-quantum-toolkit/yaqs/pull/344

docs/examples/fermi_hubbard_mpo.md

@@ -0,0 +1,61 @@
+---
+file_format: mystnb
+kernelspec:
+  name: python3
+mystnb:
+  number_source_lines: true
+  execution_timeout: 300
+---
+
+# 1D Fermi-Hubbard MPO
+
+This example shows how to build a 1D Fermi-Hubbard Hamiltonian as an MPO using {class}`~mqt.yaqs.core.data_structures.networks.MPO.fermi_hubbard_1d`.
+
+YAQS supports two representations:
+
+- **Fermionic sites** (default): one site with local dimension 4 per physical lattice site.
+  Ladder operators act on a composite ↑/↓ basis per site; this is not a Jordan–Wigner qubit chain
+  across sites, but matches the standard tensor-product embedding of site Fock spaces.
+- **Jordan-Wigner Pauli chain** (`jordan_wigner=True`): qubits in the order 1↑, 1↓, 2↑, 2↓, … with
+  local dimension 2 and full JW signs between spin orbitals. Use this mode for Pauli-string /
+  qubit simulators.
+
+The Hamiltonian (open boundaries, no chemical potential) is
+
+$$
+H = -t \sum_{i,\sigma} \left(c^\dagger_{i,\sigma} c_{i+1,\sigma} + \mathrm{h.c.}\right)
++ U \sum_i n_{i,\uparrow} n_{i,\downarrow}.
+$$
+
+## Fermionic MPO
+
+```{code-cell} ipython3
+from mqt.yaqs.core.data_structures.networks import MPO
+
+num_sites = 4
+t = 1.0
+u = 0.5
+
+h_mpo = MPO.fermi_hubbard_1d(num_sites, t=t, u=u)
+print(f"sites={h_mpo.length}, local dim={h_mpo.physical_dimension}, matrix shape={h_mpo.to_matrix().shape}")
+```
+
+The single-site basis is $|0\rangle, |\!\downarrow\rangle, |\!\uparrow\rangle, |\!\uparrow\downarrow\rangle$ (NumPy `kron` ordering for $|\!\uparrow\rangle \otimes |\!\downarrow\rangle$).
+
+## Jordan-Wigner MPO
+
+For the same model on $L$ physical sites, pass `length=2 * L` spin orbitals:
+
+```{code-cell} ipython3
+num_orbitals = 2 * num_sites
+
+h_jw = MPO.fermi_hubbard_1d(num_orbitals, t=t, u=u, jordan_wigner=True)
+print(f"orbitals={h_jw.length}, local dim={h_jw.physical_dimension}, matrix shape={h_jw.to_matrix().shape}")
+```
+
+## Relation to the Trotter circuit helper
+
+{func}`~mqt.yaqs.core.libraries.circuit_library.create_1d_fermi_hubbard_circuit` builds a **digital** Trotter circuit on separate ↑ and ↓ registers and can include a chemical potential $\mu$.
+The MPO factories above target the **analog** Hamiltonian without $\mu$ and use either fermionic operators or an interleaved JW layout.
+
+For digital simulation of the circuit model, use the circuit API; for tensor-network evolution of the Hubbard Hamiltonian, use `MPO.fermi_hubbard_1d`.

docs/index.md

@@ -39,6 +39,7 @@ examples/ensemble_evolution
 examples/solver_comparison
 examples/scheduled_jumps
 examples/transmon_emulation
+examples/fermi_hubbard_mpo
 examples/process_tomography
 examples/strong_circuit_simulation
 examples/sample_observable_digital_tjm

src/mqt/yaqs/core/data_structures/networks.py

@@ -1309,6 +1309,7 @@ class MPO:
 
     - ``MPO.ising(...)`` / ``MPO.heisenberg(...)``: qubit Pauli Hamiltonians.
     - ``MPO.hamiltonian(...)``: generic one-/two-body Pauli interactions.
+    - ``MPO.fermi_hubbard_1d(...)``: 1D Fermi-Hubbard (fermionic or Jordan-Wigner Pauli).
     - ``MPO.coupled_transmon(...)``: alternating qubit/resonator chain MPO.
     - ``from_pauli_sum(...)``: in-place build from a sum of Pauli-string terms.
     - ``identity(...)``, ``custom(...)``, ``finite_state_machine(...)``: in-place builders.
@@ -1508,6 +1509,146 @@ def heisenberg(
             n_sweeps=n_sweeps,
         )
 
+    @classmethod
+    def fermi_hubbard_1d(
+        cls,
+        length: int,
+        t: float,
+        u: float,
+        *,
+        jordan_wigner: bool = False,
+    ) -> MPO:
+        r"""Construct a 1D Fermi-Hubbard Hamiltonian MPO.
+
+        Without ``jordan_wigner``, builds the standard fermionic MPO on sites with
+        local dimension 4. The single-site basis is
+        :math:`|0\\rangle, |\\!\\downarrow\\rangle, |\\!\\uparrow\\rangle, |\\!\\uparrow\\downarrow\\rangle`
+        (NumPy ``kron`` ordering for :math:`|\\!\\uparrow\\rangle \\otimes |\\!\\downarrow\\rangle`).
+        The Hamiltonian is
+        :math:`H = -t \\sum_{i,\\sigma} (c^\\dagger_{i,\\sigma} c_{i+1,\\sigma} + \\mathrm{h.c.})
+        + U \\sum_i n_{i,\\uparrow} n_{i,\\downarrow}`.
+
+        With ``jordan_wigner=True``, builds the Jordan-Wigner Pauli-string MPO on an
+        interleaved spin chain 1↑, 1↓, 2↑, 2↓, ... (local dimension 2):
+
+        .. math::
+
+            U n_{i,\\uparrow} n_{i,\\downarrow}
+            = \\frac{U}{4} \\left(I - Z_{i,\\uparrow} - Z_{i,\\downarrow}
+            + Z_{i,\\uparrow} Z_{i,\\downarrow}\\right)
+
+            H = \\sum_i \\frac{U}{4} \\left(I - Z_{i,\\uparrow} - Z_{i,\\downarrow}
+            + Z_{i,\\uparrow} Z_{i,\\downarrow}\\right)
+            - \\frac{t}{2} \\sum_i \\left( X_{\\uparrow,i} Z_{\\downarrow,i} X_{\\uparrow,i+1}
+            + Y_{\\uparrow,i} Z_{\\downarrow,i} Y_{\\uparrow,i+1} \\right)
+            - \\frac{t}{2} \\sum_i \\left( X_{\\downarrow,i} Z_{\\uparrow,i+1} X_{\\downarrow,i+1}
+            + Y_{\\downarrow,i} Z_{\\uparrow,i+1} Y_{\\downarrow,i+1} \\right)
+
+        Without ``jordan_wigner``, the MPO uses fermionic ladder operators on composite
+        dimension-4 sites (hard-core constraint per site). Inter-site algebra matches
+        that embedding; use ``jordan_wigner=True`` for a Pauli-chain representation
+        with full Jordan-Wigner signs between spin orbitals.
+
+        In JW mode ``length`` is the number of **spin orbitals** and must be even and
+        at least 2.
+
+        Args:
+            length: Chain length. Number of fermionic sites if ``jordan_wigner`` is
+                False; number of spin orbitals (even) if True.
+            t: Hopping strength.
+            u: On-site interaction strength.
+            jordan_wigner: If True, use the JW-transformed Pauli MPO; otherwise use
+                the fermionic operator MPO.
+
+        Returns:
+            An MPO representing the 1D Fermi-Hubbard Hamiltonian.
+
+        Raises:
+            ValueError: If ``length`` is invalid for the chosen representation.
+        """
+        if jordan_wigner:
+            if length % 2 != 0 or length < 2:
+                msg = "length must be an even integer ≥ 2 (ordering: 1↑,1↓,2↑,2↓,...)."
+                raise ValueError(msg)
+            return cls._fermi_hubbard_1d_jordan_wigner(length=length, t=t, u=u)
+        return cls._fermi_hubbard_1d_fermionic(length=length, t=t, u=u)
+
+    @classmethod
+    def _fermi_hubbard_1d_fermionic(cls, length: int, t: float, u: float) -> MPO:
+        if length <= 0:
+            msg = "length must be positive."
+            raise ValueError(msg)
+
+        physical_dimension = 4
+        identity = np.eye(physical_dimension, dtype=complex)
+        zero = np.zeros_like(identity, dtype=complex)
+        c = np.array([[0, 1], [0, 0]], dtype=complex)
+        c_dag = np.array([[0, 0], [1, 0]], dtype=complex)
+        c_up = np.kron(c, np.eye(2, dtype=complex))
+        c_down = np.kron(np.eye(2, dtype=complex), c)
+        c_up_dag = np.kron(c_dag, np.eye(2, dtype=complex))
+        c_down_dag = np.kron(np.eye(2, dtype=complex), c_dag)
+        n_up = np.kron(c_dag @ c, np.eye(2, dtype=complex))
+        n_down = np.kron(np.eye(2, dtype=complex), c_dag @ c)
+        onsite = u * n_up @ n_down
+
+        # Bond layout matches ``bose_hubbard``: channels
+        # 0=identity, 1=c↑†, 2=c↓†, 3=c↑, 4=c↓, 5=accumulator.
+        tensor = np.empty((6, 6, physical_dimension, physical_dimension), dtype=object)
+        tensor[:, :] = [[zero for _ in range(6)] for _ in range(6)]
+        tensor[0, 0] = identity
+        tensor[0, 1] = c_up_dag
+        tensor[0, 2] = c_down_dag
+        tensor[0, 3] = c_up
+        tensor[0, 4] = c_down
+        tensor[0, 5] = onsite
+        tensor[1, 5] = -t * c_up
+        tensor[2, 5] = -t * c_down
+        tensor[3, 5] = -t * c_up_dag
+        tensor[4, 5] = -t * c_down_dag
+        tensor[5, 5] = identity
+
+        tensors = [np.transpose(tensor.copy(), (2, 3, 0, 1)).astype(np.complex128) for _ in range(length)]
+        tensors[0] = tensors[0][:, :, 0:1, :]
+        if length == 1:
+            tensors[0] = tensors[0][:, :, :, 5:6]
+        else:
+            tensors[-1] = tensors[-1][:, :, :, 5:6]
+
+        mpo = cls()
+        mpo.tensors = tensors
+        mpo.length = length
+        mpo.physical_dimension = physical_dimension
+        assert mpo.check_if_valid_mpo(), "MPO initialized wrong"
+        return mpo
+
+    @classmethod
+    def _fermi_hubbard_1d_jordan_wigner(cls, length: int, t: float, u: float) -> MPO:
+        num_sites = length // 2
+        terms: list[tuple[complex | float, str]] = []
+        for site in range(num_sites):
+            up, down = 2 * site, 2 * site + 1
+            terms.extend([
+                (u / 4, ""),
+                (-u / 4, f"Z{up}"),
+                (-u / 4, f"Z{down}"),
+                (u / 4, f"Z{up} Z{down}"),
+            ])
+        for site in range(num_sites - 1):
+            up, down = 2 * site, 2 * site + 1
+            up_next = 2 * (site + 1)
+            down_next = 2 * (site + 1) + 1
+            terms.extend([
+                (-t / 2, f"X{up} Z{down} X{up_next}"),
+                (-t / 2, f"Y{up} Z{down} Y{up_next}"),
+                (-t / 2, f"X{down} Z{up_next} X{down_next}"),
+                (-t / 2, f"Y{down} Z{up_next} Y{down_next}"),
+            ])
+
+        mpo = cls()
+        mpo.from_pauli_sum(terms=terms, length=length, n_sweeps=0)
+        return mpo
+
     @classmethod
     def coupled_transmon(
         cls,
@@ -1692,12 +1833,11 @@ def bose_hubbard(
 
         # build the full tensor list
         tensors = [np.transpose(tensor.copy(), (2, 3, 0, 1)).astype(np.complex128) for _ in range(length)]
-
-        # Left boundary: take only row 0
-        tensors[0] = np.transpose(tensor.copy(), (2, 3, 0, 1))[:, :, 0:1, :].astype(np.complex128)
-
-        # Right boundary: take only col 3
-        tensors[-1] = np.transpose(tensor.copy(), (2, 3, 0, 1))[:, :, :, 3:4].astype(np.complex128)
+        tensors[0] = tensors[0][:, :, 0:1, :]
+        if length == 1:
+            tensors[0] = tensors[0][:, :, :, 3:4]
+        else:
+            tensors[-1] = tensors[-1][:, :, :, 3:4]
 
         mpo = cls()
         mpo.tensors = tensors