monoprop
Tutorials

Fermi-Hubbard Model via Majorana Propagation

This notebook demonstrates how to simulate the real-time quantum dynamics of the 1D Fermi-Hubbard model using the Majorana Propagation (MP) simulator. The goal is to track how the spin-up occupancy at a given site evolves in time starting from a simple, unentangled initial state.

This 60-site example comes from the following paper, using the same model parameters, the same half-filled Néel initial state, and the spin-up occupancy at the central site as the observable. The expected runtime for the whole MP-based simulation in this notebook is around 10-15 seconds.

G. S. Hartnett, K. S. Najafi, A. Khindanov, H. Liao, M. Schutzman, M. R. Hush, M. J. Biercuk, Y. Baum,
Fast, accurate, high-resolution simulation of large-scale Fermi-Hubbard models on a digital quantum processor,
arXiv:2605.04025 (2026). https://arxiv.org/abs/2605.04025

from __future__ import annotations

import matplotlib.pyplot as plt
import numpy as np

from monoprop import Circuit, ExpGate, MajoranaPropagator
from monoprop.fermi import FermiOperator

1. Define the model

The 1D Fermi-Hubbard Hamiltonian is:

H=−t∑i,σ(ci,σ† ci+1,σ+h.c.)+U∑ini,↑ni,↓−μ∑i,σni,σH = -t \sum_{i,\sigma} \left(c^\dagger_{i,\sigma}\, c_{i+1,\sigma} + \text{h.c.}\right) + U \sum_i n_{i,\uparrow} n_{i,\downarrow} - \mu \sum_{i,\sigma} n_{i,\sigma}

The fermionic modes are indexed by site with spin states interleaved: at each site ii, spin-up occupies mode 2i2i and spin-down occupies mode 2i+12i+1. The function below returns the ordered list of local terms for the first-order Trotter decomposition — hopping bonds first, then on-site interactions, then chemical potential terms.

def mode(site, spin):
    """Return the interleaved mode index for a given site and spin (even index for spin-up, odd for spin-down)."""
    return 2 * site if spin == "up" else 2 * site + 1


def hubbard_fermion_terms(num_sites, hopping, interaction, chemical_potential):
    """Return the ordered list of local FermionOperator terms for the first-order Trotter decomposition of the 1D Hubbard model."""
    terms = []

    # nearest-neighbour hopping (both spin species)
    for site in range(num_sites - 1):
        for spin in ("up", "down"):
            left, right = mode(site, spin), mode(site + 1, spin)
            op_terms = [((left, "+"), (right, "-")), ((right, "+"), (left, "-"))]
            terms.append(
                FermiOperator(terms=op_terms, coefficients=[-hopping, -hopping])
            )

    # on-site Hubbard interaction
    for site in range(num_sites):
        up, down = mode(site, "up"), mode(site, "down")
        terms.append(
            FermiOperator(
                terms=[((up, "+"), (up, "-"), (down, "+"), (down, "-"))],
                coefficients=[interaction],
            )
        )

    # chemical potential (one term per spin-orbital)
    for site in range(num_sites):
        for spin in ("up", "down"):
            m = mode(site, spin)
            terms.append(
                FermiOperator(
                    terms=[((m, "+"), (m, "-"))],
                    coefficients=[-chemical_potential],
                )
            )

    return terms

Here are the key physical parameters for the model and computational parameters for the Trotter evolution:

  • num_sites — number of lattice sites NN (2N2N fermionic modes total).
  • hopping — hopping amplitude t=1t = 1 (energy scale).
  • interaction — on-site interaction U=−2U = -2 (attractive).
  • chemical_potential — set to 0 for half filling.
  • trotter_dt — Trotter time step Δt\Delta t.
  • trotter_steps — number of Trotter steps.
# physical parameters
num_sites = 60
hopping = 1.0
interaction = -2.0
chemical_potential = 0.0

# time-evolution parameters
trotter_dt = 0.2
trotter_steps = 30

# one mode per spin-orbital
num_qubits = 2 * num_sites

Each local Hamiltonian term is converted to Majorana generators and packed into our internal representations of the circuit, being a list of ExpGate objects.

def build_trotter_gates(num_sites, hopping, interaction, chemical_potential):
    """Convert each local Hubbard term into a Majorana generator gate for the MP simulator."""
    ferm_ops = hubbard_fermion_terms(
        num_sites, hopping, interaction, chemical_potential
    )
    return [ExpGate(term) for term in ferm_ops]


trotter_gates = build_trotter_gates(num_sites, hopping, interaction, chemical_potential)
# each gate shares the same Trotter time step
trotter_parameters = [trotter_dt for _ in trotter_gates]

We start with the half-filled Néel state as in the original paper. With start_spin="down", even sites begin spin-down occupied and odd sites spin-up occupied; the function below returns the corresponding list of occupied mode indices.

def neel_occupied_modes(num_sites, start_spin="up"):
    """Return the list of occupied mode indices for the half-filled Néel state, alternating spin between even and odd sites."""
    other = "down" if start_spin == "up" else "up"
    return [
        mode(site, start_spin if site % 2 == 0 else other) for site in range(num_sites)
    ]


occupied = neel_occupied_modes(num_sites, start_spin="down")
fermi_circuit = Circuit(
    gates=trotter_gates,
    parameters=trotter_parameters,
    initial_state=occupied,
)

We measure the spin-up occupancy observable ⟨nj,↑⟩\langle n_{j,\uparrow} \rangle at the central site j=⌊N/2⌋j = \lfloor N/2 \rfloor. In the Majorana basis:

nj,↑=12I−i2 γ2j γ2j+1n_{j,\uparrow} = \frac{1}{2}\mathbf{I} - \frac{i}{2}\,\gamma_{2j}\,\gamma_{2j+1}

The MP simulator tracks the monomial part; the constant shift +12+\tfrac{1}{2} is accounted for analytically.

def number_operator_majorana(site, spin):
    """Return the Majorana-basis representation of the number operator n_{site, spin}."""
    m = mode(site, spin)
    return FermiOperator(
        terms=[((m, "+"), (m, "-"))], coefficients=[1.0], num_modes=num_qubits
    )


obs_site = num_sites // 2
obs_spin = "up"
observable = number_operator_majorana(obs_site, obs_spin)

3. Run the simulation

MajoranaPropagator maintains the observable as a sum of weighted Majorana monomials. Each call to propagate applies one Trotter step, conjugating each monomial by the local generators and generating higher-weight monomials through operator spreading.

After each step, simulator.expectation_value() returns the observable ⟨nj,↑⟩\langle n_{j,\uparrow} \rangle and simulator.size() the number of active monomials.

simulator = MajoranaPropagator(
    observable,
    fermi_circuit.initial_state,
    cutoff=4,  # maximum Majorana monomial length kept
    cutoff_type="length",
    lower_atol=1e-4,
)

# check the initial state has 0 expectation value
print("Initial expectation value =", simulator.expectation_value())
Initial expectation value = 0.0
times = np.arange(trotter_steps + 1) * trotter_dt
values = np.empty(trotter_steps + 1)
term_counts = np.empty(trotter_steps + 1, dtype=int)

values[0] = simulator.expectation_value()
term_counts[0] = simulator.size()

for step in range(trotter_steps):
    simulator.propagate(fermi_circuit)
    values[step + 1] = simulator.expectation_value()
    term_counts[step + 1] = simulator.size()
    print(
        f"  step {step + 1:2d}/{trotter_steps}  "
        f"t={times[step + 1]:.1f}  "
        f"<n>={values[step + 1]:.4f}  "
        f"terms={term_counts[step + 1]:,}"
    )
  step  1/30  t=0.2  <n>=0.0759  terms=32
  step  2/30  t=0.4  <n>=0.2648  terms=460
  step  3/30  t=0.6  <n>=0.4747  terms=1,197
  step  4/30  t=0.8  <n>=0.6151  terms=2,451
  step  5/30  t=1.0  <n>=0.6467  terms=4,495
  step  6/30  t=1.2  <n>=0.5934  terms=7,344
  step  7/30  t=1.4  <n>=0.5140  terms=11,233
  step  8/30  t=1.6  <n>=0.4603  terms=16,588
  step  9/30  t=1.8  <n>=0.4506  terms=23,089
  step 10/30  t=2.0  <n>=0.4703  terms=30,643
  step 11/30  t=2.2  <n>=0.4934  terms=39,384
  step 12/30  t=2.4  <n>=0.5039  terms=49,663
  step 13/30  t=2.6  <n>=0.5031  terms=61,319
  step 14/30  t=2.8  <n>=0.5014  terms=73,899
  step 15/30  t=3.0  <n>=0.5052  terms=87,388
  step 16/30  t=3.2  <n>=0.5116  terms=101,969
  step 17/30  t=3.4  <n>=0.5130  terms=118,202
  step 18/30  t=3.6  <n>=0.5058  terms=134,297
  step 19/30  t=3.8  <n>=0.4941  terms=149,958


  step 20/30  t=4.0  <n>=0.4858  terms=166,851
  step 21/30  t=4.2  <n>=0.4860  terms=183,546
  step 22/30  t=4.4  <n>=0.4940  terms=199,652
  step 23/30  t=4.6  <n>=0.5038  terms=215,819
  step 24/30  t=4.8  <n>=0.5100  terms=230,997
  step 25/30  t=5.0  <n>=0.5102  terms=245,309
  step 26/30  t=5.2  <n>=0.5056  terms=258,395
  step 27/30  t=5.4  <n>=0.4994  terms=270,983


  step 28/30  t=5.6  <n>=0.4945  terms=281,947
  step 29/30  t=5.8  <n>=0.4926  terms=291,284
  step 30/30  t=6.0  <n>=0.4944  terms=299,435

4. Results

The observable starts at zero (the central site is spin-down occupied in the Néel state), rises as hopping spreads spin-up electrons across the chain, and oscillates around the thermal value of 0.50.5 (dashed line). The right panel shows the number of active Majorana monomials, which grows with time until truncation saturates the expansion.

fig, axes = plt.subplots(1, 2, figsize=(12, 4), constrained_layout=True)

# left: expectation value vs time
ax = axes[0]
ax.plot(
    times,
    values,
    color="#33766f",
    linewidth=2,
    marker="s",
    markersize=6,
    markerfacecolor="none",
    markeredgecolor="#33766f",
    label="Majorana Propagation",
)
ax.axhline(0.5, color="gray", linestyle="--", linewidth=1, label="thermal")
ax.set_xlabel("Time $t$", fontsize=12)
ax.set_ylabel(rf"$\langle n_{{{obs_site},\uparrow}} \rangle$", fontsize=12)
ax.set_title(
    f"{num_sites}-site Hubbard model  ($U={interaction}$, $t={hopping}$)", fontsize=11
)
ax.set_xlim(0, times[-1])
ax.set_ylim(-0.05, 0.75)
ax.spines[["top", "right"]].set_visible(False)
ax.legend(frameon=False, fontsize=9)

# right: active monomial count (proxy for computational cost)
ax2 = axes[1]
ax2.plot(
    times,
    term_counts,
    color="#7d1cb1",
    linewidth=2,
    marker="o",
    markersize=5,
    markerfacecolor="none",
)
ax2.set_xlabel("Time $t$", fontsize=12)
ax2.set_ylabel("Active Majorana monomials", fontsize=12)
ax2.set_title("Number of terms vs time", fontsize=11)
ax2.set_xlim(0, times[-1])
ax2.spines[["top", "right"]].set_visible(False)

plt.show()
png

Chemistry simulations with Majorana propagation

Previous Page

Kicked Ising Model on IBM Eagle (127 Qubits)

Next Page

On this page

1. Define the model3. Run the simulation4. Results