Spin and charge dynamics in the Fermi-Hubbard model

Usage estimate: 5 minutes on a Heron r3 processor. (NOTE: This is an estimate only. Your runtime might vary.)

Background

The Fermi-Hubbard model describes electrons hopping on a lattice with an energy penalty for double occupancy. On a one-dimensional chain of $L$ sites with open boundaries, the Hamiltonian is defined as follows:

H = -\sum_{i=0}^{L-2} \sum_{\sigma \in \{\uparrow,\downarrow\}} t_{h,i} \left( c^\dagger_{i\sigma} c_{i+1,\sigma} + \mathrm{h.c.} \right) + \sum_{i=0}^{L-1} U_i\, n_{i\uparrow} n_{i\downarrow} - \sum_{i=0}^{L-1} \sum_\sigma \mu_i\, n_{i\sigma},

where $c^\dagger_{i\sigma}$ creates a fermion of spin $\sigma$ on site $i$ and $n_{i\sigma} = c^\dagger_{i\sigma} c_{i\sigma}$ counts them. The three terms are the hopping amplitude $t_{h,i}$ on the bond between sites $i$ and $i+1$ , the on-site interaction $U_i$ , and the chemical potential $\mu_i$ . The tutorial measures energy in units of $t_h$ and time in units of $1/t_h$ .

Real-time dynamics of this model are hard to compute classically. Exact diagonalization stops near $L \lesssim 14$ sites, and after a quench, the entanglement entropy of a tensor-network solver grows linearly in time, so the bond dimension needed for fixed accuracy grows exponentially.

In one dimension, the elementary excitations are spinons, which carry spin and no charge, and holons, which carry charge and no spin. The two propagate at different speeds, so a localized disturbance splits into a fast front and a slow one. This effect is called spin-charge separation, and this tutorial measures it [1].

This tutorial uses the Dynamics Simulator Qiskit Function, which takes the Hamiltonian, an initial state, a Trotter schedule, and a list of observables, and returns expectation values. It handles the Jordan-Wigner mapping, Trotterization, qubit layout, and error suppression internally, so there are no circuits to build or transpile in this tutorial.

Requirements

Before you start this tutorial, make sure that you have the following packages installed:

Qiskit SDK v2.3.0 or later (pip install qiskit)
Qiskit Runtime v0.50.0 or later (pip install qiskit-ibm-runtime)
Qiskit Functions Catalog v0.21.0 or later (pip install qiskit-ibm-catalog)
NumPy v2.5.1 or later (pip install numpy)
Matplotlib v3.10.3 or later (pip install matplotlib)
Q-CTRL Visualizer v11.0.0 or later (pip install qctrl-visualizer)

Setup

First, authenticate using your IBM Quantum API key. Then, select the Qiskit Function as follows. (This code assumes that you already saved your account to your local environment.)

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
from qiskit_ibm_catalog import QiskitFunctionsCatalog
from qiskit_ibm_runtime import QiskitRuntimeService
import qctrlvisualizer as qv

plt.style.use(qv.get_qctrl_style())

catalog = QiskitFunctionsCatalog(channel="ibm_quantum_platform")

# Access Function
sim_qf = catalog.load("q-ctrl/dynamics-simulator")

Step 1: Map classical inputs to a quantum problem

Define the Hamiltonian

The function takes the model as a dictionary with three keys. Each one accepts either a single float, which applies uniformly along the chain, or a list that specifies the value site by site. Lists have $L-1$ entries for hopping_amplitude, because an open chain has one fewer bond than sites, and $L$ entries for on_site_interaction and chemical_potential.

This tutorial uses a uniform chain of 31 sites with a repulsive interaction $U/t_h = 8$ . A uniform chemical potential contributes only a global phase when the particle number is fixed, so this example sets it to zero. The chemical potential matters only when it varies from site to site, for example, to add disorder or a single impurity.

site_count = 31
hopping_amplitude = 1.0
on_site_interaction = 8.0
chemical_potential = 0.0

model = {
    "hopping_amplitude": [hopping_amplitude] * (site_count - 1),
    "on_site_interaction": [on_site_interaction] * site_count,
    "chemical_potential": [chemical_potential] * site_count,
}

print(
    f"{site_count} sites, {2 * site_count} qubits, U/t_h = {on_site_interaction:.0f}"
)

Output:

31 sites, 62 qubits, U/t_h = 8

Set the Trotter schedule

The simulation dictionary holds the time step ( $\Delta t$ ) and the number of Trotter steps, and the total evolution time is their product. Larger steps reach later times at fixed circuit depth but incur more Trotter error, which grows with $U$ because the Hamiltonian terms commute less. The published hardware study used $\Delta t = 0.15$ for $U/t_h = 8$ [2], and this tutorial uses the same value. Circuit depth grows linearly in the step count and does not depend on $L$ at all.

time_step = 0.15
step_count = 30

simulation = {"time_step": time_step, "step_count": step_count}

print(
    f"Total evolution time: t = {time_step * step_count:.2f} (units of 1/t_h)"
)

Output:

Total evolution time: t = 4.50 (units of 1/t_h)

Prepare the initial state

Each site holds two spin orbitals, so an $L$ -site chain uses $2L$ qubits. The initial state is a Fock occupation string with one character per site:

\begin{aligned} \texttt{empty (e)} &\longleftrightarrow \left(n_{i,\uparrow}, n_{i,\downarrow}\right) = (0,0), \\ \texttt{up (u)} &\longleftrightarrow \left(n_{i,\uparrow}, n_{i,\downarrow}\right) = (1,0), \\ \texttt{down (d)} &\longleftrightarrow \left(n_{i,\uparrow}, n_{i,\downarrow}\right) = (0,1), \\ \texttt{both (b)} &\longleftrightarrow \left(n_{i,\uparrow}, n_{i,\downarrow}\right) = (1,1). \end{aligned}

The string is read least significant site first, so the rightmost character is site 0 and the physical site index of character $k$ counted from the left is $i = L - 1 - k$ . The string udud therefore means site 0 is d, site 1 is u, site 2 is d, and site 3 is u.

Next, prepare a half-filled Néel state, alternating u and d, with the central site replaced by a vacancy. The vacancy is the disturbance whose spreading is tracked.

vacancy_site = site_count // 2

# Build in increasing site order, then reverse: the leftmost character of the
# string is the highest site index.
site_occupations = [
    "u" if site % 2 == 0 else "d" for site in range(site_count)
]
site_occupations[vacancy_site] = "e"
initial_state = "".join(reversed(site_occupations))

print(f"Vacancy at site {vacancy_site}")
print(initial_state)

Output:

Vacancy at site 15
udududududududueudududududududu

Build the observables

Each observable is a dictionary that represents a weighted sum of products of occupation operators:

O = \sum_{\alpha} c_{\alpha} \prod_{(\sigma, j) \in T_{\alpha}} n_{j,\sigma}.

Every entry of operators is one term $T_{\alpha}$ , written as a list of [name, site] pairs that multiply together, and the matching entry of coefficients is its weight. The two names are "Nup" for $n_{j,\uparrow}$ and "Ndn" for $n_{j,\downarrow}$ . A single one-point observable $n_{3,\uparrow}$ is written as follows:

{"operators": [[["Nup", 3]]], "coefficients": [1.0]}

The next cell prints the smallest complete set of inputs the function accepts, which is useful as a shape reference while reading the code that follows.

minimal_inputs = {
    "model": {
        "hopping_amplitude": 1.0,
        "on_site_interaction": [2.0, 2.5],
        "chemical_potential": 0.5,
    },
    "simulation": {"step_count": 5, "time_step": 0.1},
    "observables": [{"operators": [[["Nup", 0]]], "coefficients": [1.0]}],
    "initial_state": "ud",
    "shots": 4096,
}

for key, value in minimal_inputs.items():
    print(f"{key}: {value}")

Output:

model: {'hopping_amplitude': 1.0, 'on_site_interaction': [2.0, 2.5], 'chemical_potential': 0.5}
simulation: {'step_count': 5, 'time_step': 0.1}
observables: [{'operators': [[['Nup', 0]]], 'coefficients': [1.0]}]
initial_state: ud
shots: 4096

Three groups of observables are needed. The charge density $n_i = n_{i\uparrow} + n_{i\downarrow}$ and the local spin $S^z_i = (n_{i\uparrow} - n_{i\downarrow})/2$ are both two-term sums over the same operators, differing only in their coefficients. The spin correlator $S^z_i S^z_{i_*}$ between site $i$ and the vacancy site expands into four two-point terms with coefficients $\pm 1/4$ .

Occupation operators commute, so the order of factors within a term is irrelevant, but the order of the observables list fixes the order of the returned values. The code records an index for each group while it builds the list. The correlator skips the vacancy site, where the product would repeat a factor.

observables = []
charge_indices = []
spin_indices = []
correlator_indices = []

# Charge density: n_i = n_{i,up} + n_{i,down}
for site in range(site_count):
    charge_indices.append(len(observables))
    observables.append(
        {
            "operators": [[["Nup", site]], [["Ndn", site]]],
            "coefficients": [1.0, 1.0],
        }
    )

# Local spin: S^z_i = (n_{i,up} - n_{i,down}) / 2
for site in range(site_count):
    spin_indices.append(len(observables))
    observables.append(
        {
            "operators": [[["Nup", site]], [["Ndn", site]]],
            "coefficients": [0.5, -0.5],
        }
    )

# Spin correlator: S^z_i S^z_{i_*}, expanded into four two-point terms.
correlator_sites = [
    site for site in range(site_count) if site != vacancy_site
]
for site in correlator_sites:
    correlator_indices.append(len(observables))
    observables.append(
        {
            "operators": [
                [["Nup", site], ["Nup", vacancy_site]],
                [["Nup", site], ["Ndn", vacancy_site]],
                [["Ndn", site], ["Nup", vacancy_site]],
                [["Ndn", site], ["Ndn", vacancy_site]],
            ],
            "coefficients": [0.25, -0.25, -0.25, 0.25],
        }
    )

print(f"{len(observables)} observables")
print("Example correlator:", observables[correlator_indices[0]])

Output:

92 observables
Example correlator: {'operators': [[['Nup', 0], ['Nup', 15]], [['Nup', 0], ['Ndn', 15]], [['Ndn', 0], ['Nup', 15]], [['Ndn', 0], ['Ndn', 15]]], 'coefficients': [0.25, -0.25, -0.25, 0.25]}

Step 2: Optimize problem for quantum hardware execution

This tutorial has no transpilation step. The function performs the following steps internally [2], and you cannot configure them from the notebook:

Maps the fermionic modes to qubits under an ordering chosen to keep Pauli weights low on a heavy-hex lattice
Uses a fermionic SWAP network so every hopping term becomes a two-qubit rotation
Embeds the chain along a path that avoids weak qubits
Applies dynamical decoupling and Pauli twirling during compilation

You choose the backend and the shot count. The function takes a backend name as a string, so this tutorial selects the least busy device large enough for the chain and passes its name.

The smallest signal that you need to resolve sets the shot count. The charge tracer reaches amplitudes near $0.15$ , but the connected spin tracer peaks around $0.04$ , roughly four times smaller. The published study used 20,000 shots per circuit for this reason [2].

service = QiskitRuntimeService()
backend = service.backend("ibm_boston")
backend_name = backend.name
shot_count = 20_000

print(f"{backend_name}, {shot_count} shots per circuit")

Output:

ibm_boston, 20000 shots per circuit

Step 3: Execute by using Qiskit primitives

Next, submit the simulation. run returns a job handle immediately rather than blocking.

job = sim_qf.run(
    model=model,
    simulation=simulation,
    observables=observables,
    initial_state=initial_state,
    shots=shot_count,
    backend_name=backend_name,
)

print(job.job_id, job.status())

Output:

6261b073-be19-48a1-be60-2f369b25613c QUEUED

Use job.status() to follow progress and job.result() to block until the job finishes. The result is a mapping whose observable_expectation_values entry is indexed first by observable, in the order they were supplied, and then by time. There are step_count + 1 time points because the value at $t=0$ is included.

result = job.result()

print("Result keys:", sorted(result.keys()))

Output:

Result keys: ['observable_expectation_values', 'observables', 'provider_job_ids', 'warnings']

Step 4: Post-process and return result in desired classical format

Reconstruct the correlators

Split the flat array back into the three groups using the indices recorded earlier, then form the two tracer correlators. The charge tracer is the change in per-site density relative to $t=0$ :

C^c_i(t) = \langle n_{i\uparrow}(t) + n_{i\downarrow}(t) \rangle - \langle n_{i\uparrow}(0) + n_{i\downarrow}(0) \rangle.

The spin tracer is the connected correlator between each site and the vacancy site $i_*$ :

C^s_i(t) = 4 \left( \langle S^z_i S^z_{i_*} \rangle - \langle S^z_i \rangle \langle S^z_{i_*} \rangle \right).

Both are built from Fock-basis operators, so they come from the same measured bitstrings and the second one costs no extra hardware execution.

When you constructed the observables, the code recorded the position of each charge (charge_indices), spin (spin_indices), and spin-correlation (correlator_indices) operator in the list. Because the simulator returns the expectation values in the same order, you can use these positions to extract the corresponding results. For the spin correlations, the code also uses correlator_sites to place each result at the correct lattice position, and leaves the vacancy as NaN.

expectation_values = np.asarray(
    result["observable_expectation_values"], dtype=float
)
times = np.arange(expectation_values.shape[1]) * time_step

charge_density = expectation_values[charge_indices]
local_spin = expectation_values[spin_indices]
spin_correlator = np.full_like(charge_density, np.nan)
spin_correlator[correlator_sites] = expectation_values[correlator_indices]

# Charge tracer: change in per-site density relative to t = 0.
charge_tracer = charge_density - charge_density[:, [0]]

# Spin tracer: connected correlator with the vacancy site.
spin_tracer = 4.0 * (spin_correlator - local_spin * local_spin[vacancy_site])
spin_tracer[vacancy_site] = np.nan

print(
    f"{expectation_values.shape[0]} observables, "
    f"{expectation_values.shape[1]} time points up to t = {times[-1]:.2f}"
)
print(
    f"Peak amplitudes: charge {np.nanmax(np.abs(charge_tracer)):.3f}, "
    f"spin {np.nanmax(np.abs(spin_tracer)):.3f}"
)

Output:

92 observables, 31 time points up to t = 4.50
Peak amplitudes: charge 0.978, spin 0.168

Visualize the wavefronts

If you plot each correlator as a heatmap of space and time, the light cone is directly visible. Time flows downward, and the horizontal axis is the position along the chain.

def plot_tracer(axis, tracer, title):
    """Draw one tracer correlator as a space-time heatmap, time flowing down."""
    # Scale from the off-centre sites. The vacancy site carries an order-one
    # change in charge density and is left to clip rather than set the scale.
    off_center = np.ones(site_count, dtype=bool)
    off_center[vacancy_site] = False
    limit = np.nanpercentile(np.abs(tracer[off_center]), 98)

    colormap = qv.QCTRL_DIVERGENT_COLORMAP.reversed().copy()
    colormap.set_bad("0.82")  # the spin tracer has no value at the vacancy

    image = axis.imshow(
        tracer.T,
        aspect="auto",
        origin="upper",
        cmap=colormap,
        vmin=-limit,
        vmax=limit,
        extent=(-0.5, site_count - 0.5, times[-1], times[0]),
        interpolation="nearest",
    )

    axis.annotate(
        "vacancy",
        xy=(vacancy_site, times[-1]),
        xytext=(vacancy_site, times[-1] + 0.55),
        ha="center",
        va="top",
        fontsize=8,
        annotation_clip=False,
        arrowprops=dict(arrowstyle="-", lw=0.8, color="0.4"),
    )

    axis.set_title(title)
    axis.set_xlabel("Site")
    return image


fig, axes = plt.subplots(1, 2, figsize=(11, 4.5), sharey=True)

charge_image = plot_tracer(
    axes[0], charge_tracer, r"Charge tracer $C^{c}_{i}(t)$"
)
spin_image = plot_tracer(axes[1], spin_tracer, r"Spin tracer $C^{s}_{i}(t)$")

axes[0].set_ylabel(r"Time $t$  $(1/t_h)$")
fig.colorbar(charge_image, ax=axes[0], fraction=0.046)
fig.colorbar(spin_image, ax=axes[1], fraction=0.046)
fig.suptitle(f"$L$ = {site_count},  $U/t_h$ = {on_site_interaction:.0f}")
fig.tight_layout()
plt.show()

Output:

Extract the wavefront velocities

Turn the heatmaps into numbers by following the detection protocol from the published study [2]. At each time slice, symmetrize the spatial profile about the vacancy, which averages away the antisymmetric part of the hardware noise, normalize by the largest off-center amplitude, and record the outermost site above a threshold. Detections are then cut into runs of consecutive time steps, and only the longest run is fit, which keeps the early-time transient out of the result.

def front_distance(profile, center, threshold):
    """Distance from `center` to the outermost site above a relative threshold."""
    values = np.nan_to_num(np.asarray(profile, dtype=float))
    length = values.size

    symmetrized = values.copy()
    for site in range(length):
        mirror = 2 * center - site
        if 0 <= mirror < length:
            symmetrized[site] = 0.5 * (values[site] + values[mirror])

    amplitude = np.abs(symmetrized)
    off_center = np.ones(length, dtype=bool)
    off_center[center] = False

    peak = amplitude[off_center].max()
    if peak <= 0.0:
        return np.nan
    above = np.flatnonzero((amplitude / peak > threshold) & off_center)
    return float(np.abs(above - center).max()) if above.size else np.nan


def longest_run(mask):
    """Indices of the longest run of consecutive True values."""
    best, current = [], []
    for index, keep in enumerate(mask):
        current = current + [index] if keep else []
        if len(current) > len(best):
            best = current
    return np.array(best, dtype=int)


def front_velocity(tracer, center, threshold=0.3, quantile=0.75, min_block=5):
    """Fit front position against time, using only the longest contiguous run."""
    distances = np.array(
        [
            front_distance(tracer[:, k], center, threshold)
            for k in range(tracer.shape[1])
        ]
    )
    edge = max(center, tracer.shape[0] - 1 - center)
    usable = (
        np.isfinite(distances) & (distances >= 2) & (distances <= edge - 2)
    )

    block = longest_run(usable)
    if block.size < min_block:
        return np.nan, distances, np.array([], dtype=int), np.nan

    candidates = [
        stats.theilslopes(distances[block[:k]], times[block[:k]]).slope
        for k in range(min_block, block.size + 1)
    ]
    cutoff = np.quantile(candidates, quantile)
    keep = max(i for i, value in enumerate(candidates) if value >= cutoff)

    fit_block = block[: keep + min_block]
    slope, intercept = np.polyfit(times[fit_block], distances[fit_block], 1)
    return float(slope), distances, fit_block, float(intercept)

The fitted velocity depends on the detection threshold and the trimming quantile, and that sensitivity dominates the shot noise. Scan both and report the median with the full envelope, as the published analysis does.

thresholds = [0.15, 0.225, 0.30, 0.375, 0.45]
quantiles = [0.60, 0.675, 0.75, 0.825, 0.90]


def velocity_envelope(tracer, center):
    """Median, minimum, and maximum velocity over the parameter scan."""
    values = [
        front_velocity(tracer, center, threshold=p, quantile=q)[0]
        for p in thresholds
        for q in quantiles
    ]
    values = np.array([v for v in values if np.isfinite(v)])
    if values.size == 0:
        return np.nan, np.nan, np.nan
    return float(np.median(values)), float(values.min()), float(values.max())


charge_velocity, charge_low, charge_high = velocity_envelope(
    charge_tracer, vacancy_site
)
spin_velocity, spin_low, spin_high = velocity_envelope(
    spin_tracer, vacancy_site
)

_, charge_front, charge_block, charge_intercept = front_velocity(
    charge_tracer, vacancy_site
)
_, spin_front, spin_block, spin_intercept = front_velocity(
    spin_tracer, vacancy_site
)

print(f"v_c = {charge_velocity:.2f}  [{charge_low:.2f}, {charge_high:.2f}]")
print(f"v_s = {spin_velocity:.2f}  [{spin_low:.2f}, {spin_high:.2f}]")
print(f"v_s / v_c = {spin_velocity / charge_velocity:.2f}")

Output:

v_c = 1.93  [1.81, 2.18]
v_s = 1.61  [1.23, 2.12]
v_s / v_c = 0.83

PUBLISHED = {"charge": 1.78, "spin": 1.46}

edge = max(vacancy_site, site_count - 1 - vacancy_site)
fig, axis = plt.subplots(figsize=(7, 4.6))

axis.axhspan(0, 2, color="0.92", zorder=0)
axis.axhspan(edge - 2, edge + 2, color="0.92", zorder=0)

series = [
    (
        "charge",
        "$v_c$",
        charge_front,
        charge_block,
        charge_velocity,
        charge_intercept,
        charge_low,
        charge_high,
        "C0",
        "o",
    ),
    (
        "spin",
        "$v_s$",
        spin_front,
        spin_block,
        spin_velocity,
        spin_intercept,
        spin_low,
        spin_high,
        "C1",
        "s",
    ),
]

for name, sym, front, block, v, b, lo, hi, color, marker in series:
    excluded = np.setdiff1d(np.arange(front.size), block)
    axis.plot(
        times[excluded],
        front[excluded],
        marker,
        color=color,
        markersize=5,
        markerfacecolor="none",
        alpha=0.45,
        zorder=2,
    )
    axis.plot(
        times[block],
        front[block],
        marker,
        color=color,
        markersize=5,
        label=f"{name}, {sym} = {v:.2f} [{lo:.2f}, {hi:.2f}]",
        zorder=3,
    )

    span = times[block]
    anchor = v * span[0] + b
    axis.plot(
        span,
        v * (span - span[0]) + anchor,
        "-",
        color=color,
        lw=1.6,
        zorder=4,
    )
    axis.fill_between(
        span,
        lo * (span - span[0]) + anchor,
        hi * (span - span[0]) + anchor,
        color=color,
        alpha=0.15,
        zorder=1,
    )
    axis.plot(
        span,
        PUBLISHED[name] * (span - span[0]) + anchor,
        "--",
        color=color,
        lw=1.2,
        alpha=0.8,
        zorder=4,
    )

axis.plot([], [], "--", color="0.4", lw=1.2, label="published [2]")
axis.plot(
    [],
    [],
    "o",
    color="0.5",
    markerfacecolor="none",
    label="excluded from fit",
)

axis.set_xlabel(r"Time $t$  $(1/t_h)$")
axis.set_ylabel("Front distance from vacancy (sites)")
axis.set_xlim(0, times[-1])
axis.set_ylim(
    0, min(edge + 1, np.nanmax(np.r_[charge_front, spin_front]) + 3)
)
axis.legend(loc="upper left", frameon=False, fontsize=9)
axis.set_title(
    f"L = {site_count}, $U/t_h$ = {on_site_interaction:.0f}, front propagation",
    fontsize=11,
)

fig.tight_layout()
plt.show()

Output:

Drawing the detections back onto the heatmaps shows whether the algorithm tracked the leading edge. Filled markers are the detections that entered the fit. Open markers were cut by the masking rules, either because they sit too close to the vacancy, too close to the chain end, or were late enough that the front has broadened and its position is no longer sharp. The dashed line is the fit over the times it used, and its slope is the velocity in each panel title.

The vacancy column is marked below the axis. The charge clips there because its value is of order one, and spin has no value there, so it is drawn in gray.

The charge cone opens wider than the spin cone. Both disturbances leave the same defect at the same moment and travel outward at different speeds.

def plot_tracer_front(axis, tracer, front, block, velocity, intercept, title):
    """Overlay the detected wavefront and the fitted cone on a tracer heatmap."""
    image = plot_tracer(axis, tracer, title)

    excluded = np.setdiff1d(np.arange(front.size), block)
    for index, filled in ((block, True), (excluded, False)):
        good = index[np.isfinite(front[index])]
        if good.size == 0:
            continue
        for sign in (-1, 1):
            axis.plot(
                vacancy_site + sign * front[good],
                times[good],
                "o",
                markersize=4,
                color="k",
                markerfacecolor="k" if filled else "none",
                markeredgewidth=0.9,
                alpha=1.0 if filled else 0.4,
            )

    span = times[block]
    for sign in (-1, 1):
        axis.plot(
            vacancy_site + sign * (velocity * span + intercept),
            span,
            "--",
            color="k",
            lw=1.5,
        )

    axis.set_title(f"{title}    $v$ = {velocity:.2f}")
    return image


fig, axes = plt.subplots(1, 2, figsize=(11, 4.5), sharey=True)

charge_image = plot_tracer_front(
    axes[0],
    charge_tracer,
    charge_front,
    charge_block,
    charge_velocity,
    charge_intercept,
    r"Charge tracer $C^{c}_{i}(t)$",
)
spin_image = plot_tracer_front(
    axes[1],
    spin_tracer,
    spin_front,
    spin_block,
    spin_velocity,
    spin_intercept,
    r"Spin tracer $C^{s}_{i}(t)$",
)

axes[0].set_ylabel(r"Time $t$  $(1/t_h)$")
fig.colorbar(charge_image, ax=axes[0], fraction=0.046)
fig.colorbar(spin_image, ax=axes[1], fraction=0.046)
fig.suptitle(f"$L$ = {site_count},  $U/t_h$ = {on_site_interaction:.0f}")
fig.tight_layout()
plt.show()

Output:

The charge front outruns the spin front. A hole in a one-dimensional Mott background does not travel as an electron. Instead, it fractionalizes into a holon carrying the charge and a spinon carrying the spin disturbance, and at $U/t_h = 8$ the spinon is the slower of the two.

The ordering depends on the interaction. At $U = 0$ , the model is free fermions with $v_c = 2$ and $v_s = 4$ , so spin is faster. Repulsion suppresses spin transport while leaving charge transport nearly untouched, and by $U/t_h = 8$ the published hardware values are $v_c = 1.78$ and $v_s = 1.46$ [2]. Rerunning with several values of $U$ and plotting $v_s / v_c$ maps the crossover.

Note these caveats on the numbers. Threshold-based detection is biased slightly low, so a charge velocity a few percent under the published value is expected. The spin envelope is wide because the connected correlator sits closer to the noise floor, which is why the median of the scan rather than a single fit is quoted.

References

[1] Essler, F. H. L., Frahm, H., Gohmann, F., Klumper, A., & Korepin, V. E. (2005). The One-Dimensional Hubbard Model. Cambridge University Press.

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

Was this page helpful?

Report a bug, typo, or request content on GitHub.