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 four parallel arrays required by `propagate`: Majorana index tuples (`majoranas`), real prefactors (`gen_coeffs`), a group index mapping each generator to its local term (`param_inds`), and the shared time step per group (`parameters`). The real prefactor is obtained by factoring out the imaginary phase $(-i)^{w(w-1)/2}$ from each complex Majorana coefficient, where $w$ is the monomial weight."
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,205
  step  4/30  t=0.8  <n>=0.6151  terms=2,455
  step  5/30  t=1.0  <n>=0.6467  terms=4,515
  step  6/30  t=1.2  <n>=0.5934  terms=7,449


  step  7/30  t=1.4  <n>=0.5140  terms=11,434


  step  8/30  t=1.6  <n>=0.4603  terms=16,849
  step  9/30  t=1.8  <n>=0.4507  terms=23,542
  step 10/30  t=2.0  <n>=0.4706  terms=31,293


  step 11/30  t=2.2  <n>=0.4938  terms=40,275


  step 12/30  t=2.4  <n>=0.5041  terms=51,001
  step 13/30  t=2.6  <n>=0.5032  terms=63,024
  step 14/30  t=2.8  <n>=0.5015  terms=76,534


  step 15/30  t=3.0  <n>=0.5054  terms=91,052


  step 16/30  t=3.2  <n>=0.5119  terms=106,738


  step 17/30  t=3.4  <n>=0.5135  terms=124,028
  step 18/30  t=3.6  <n>=0.5061  terms=141,868
  step 19/30  t=3.8  <n>=0.4937  terms=160,579
  step 20/30  t=4.0  <n>=0.4846  terms=180,293


  step 21/30  t=4.2  <n>=0.4846  terms=200,343
  step 22/30  t=4.4  <n>=0.4931  terms=220,705
  step 23/30  t=4.6  <n>=0.5040  terms=241,125
  step 24/30  t=4.8  <n>=0.5111  terms=261,654


  step 25/30  t=5.0  <n>=0.5109  terms=281,095
  step 26/30  t=5.2  <n>=0.5050  terms=300,441


  step 27/30  t=5.4  <n>=0.4974  terms=318,630


  step 28/30  t=5.6  <n>=0.4923  terms=335,810


  step 29/30  t=5.8  <n>=0.4921  terms=351,551


  step 30/30  t=6.0  <n>=0.4959  terms=365,808

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

Majoranic chemistry simulations

Previous Page

Kicked Ising Model on IBM Eagle (127 Qubits)

Next Page

On this page

1. Define the model3. Run the simulation4. Results