Spin at charge dynamics sa Fermi-Hubbard model

Tantiya ng paggamit: 5 minuto sa isang Heron r3 processor. (PAALALA: Tantiya lang ito. Maaaring iba ang runtime mo.)

Background

Inilalarawan ng Fermi-Hubbard model ang mga electron na tumatalon sa isang lattice na may energy penalty para sa double occupancy. Sa isang one-dimensional na chain ng $L$ site na may open boundary, ganito ang depinisyon ng Hamiltonian:

$$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},$$

kung saan ang $c^\dagger_{i\sigma}$ ay lumilikha ng fermion na may spin $\sigma$ sa site $i$ at ang $n_{i\sigma} = c^\dagger_{i\sigma} c_{i\sigma}$ ay nagbibilang sa kanila. Ang tatlong term ay ang hopping amplitude $t_{h,i}$ sa bond sa pagitan ng mga site na $i$ at $i+1$, ang on-site interaction na $U_i$, at ang chemical potential na $\mu_i$. Sinusukat ng tutorial ang energy sa yunit ng $t_h$ at ang oras sa yunit ng $1/t_h$.

Mahirap kuwentahin nang classical ang real-time dynamics ng model na ito. Humihinto ang exact diagonalization malapit sa $L \lesssim 14$ na site, at pagkatapos ng quench, linear na lumalaki sa oras ang entanglement entropy ng isang tensor-network solver, kaya exponential na lumalaki ang bond dimension na kailangan para sa nakatakdang accuracy.

Sa isang dimensyon, ang mga elementary excitation ay ang mga spinon, na may dalang spin at walang charge, at ang mga holon, na may dalang charge at walang spin. Magkaiba ang bilis ng pagkalat ng dalawa, kaya nahahati ang isang naka-localize na gambala sa mabilis na front at mabagal na front. Tinatawag itong spin-charge separation, at sinusukat ito ng tutorial na ito [1].

Gumagamit ang tutorial na ito ng Dynamics Simulator Qiskit Function, na tumatanggap ng Hamiltonian, isang initial state, isang Trotter schedule, at isang listahan ng mga observable, at nagbabalik ng mga expectation value. Ito na ang bahala sa Jordan-Wigner mapping, Trotterization, qubit layout, at error suppression sa loob, kaya walang circuit na bubuuin o ita-transpile sa tutorial na ito.

Requirements

Bago mo simulan ang tutorial na ito, siguraduhing naka-install ang mga sumusunod na package:

• Qiskit SDK v2.3.0 o mas bago (pip install qiskit)

• Qiskit Runtime v0.50.0 o mas bago (pip install qiskit-ibm-runtime)

• Qiskit Functions Catalog v0.21.0 o mas bago (pip install qiskit-ibm-catalog)

• NumPy v2.5.1 o mas bago (pip install numpy)

• Matplotlib v3.10.3 o mas bago (pip install matplotlib)

• Q-CTRL Visualizer v11.0.0 o mas bago (pip install qctrl-visualizer)

Setup

Una, mag-authenticate gamit ang iyong IBM Quantum API key. Pagkatapos, piliin ang Qiskit Function tulad ng sumusunod. (Ipinapalagay ng code na ito na na-save mo na ang account mo sa local environment mo.)

# Added by doQumentation — installs the packages this notebook needs if they are missing
import importlib.util

_needed = {"matplotlib": "matplotlib", "numpy": "numpy", "qctrlvisualizer": "qctrlvisualizer", "qiskit_ibm_catalog": "qiskit-ibm-catalog", "qiskit_ibm_runtime": "qiskit-ibm-runtime", "scipy": "scipy"}
_missing = [pip for module, pip in _needed.items()
            if importlib.util.find_spec(module) is None]
# One at a time, so a package that fails to install does not block the others
for _pip in _missing:
    %pip install -q {_pip}
if not _missing:
    print("\u2713 All required packages are installed")

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: I-map ang mga classical input sa isang quantum problem

I-define ang Hamiltonian

Tinatanggap ng function ang model bilang dictionary na may tatlong key. Tumatanggap ang bawat isa ng isang float, na pantay na ilalapat sa buong chain, o isang listahan na tumutukoy sa halaga kada site. May $L-1$ na entry ang mga listahan para sa hopping_amplitude, dahil may isang bond na mas kaunti ang open chain kaysa sa bilang ng site, at $L$ na entry para sa on_site_interaction at chemical_potential.

Gumagamit ang tutorial na ito ng uniform na chain na 31 site na may repulsive interaction na $U/t_h = 8$. Global phase lang ang naiaambag ng uniform na chemical potential kapag nakapirmi ang bilang ng particle, kaya ginawa itong zero ng halimbawang ito. Mahalaga lang ang chemical potential kapag nag-iiba ito sa bawat site, halimbawa para magdagdag ng disorder o ng isang 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}"
)

31 sites, 62 qubits, U/t_h = 8

Itakda ang Trotter schedule

Hawak ng simulation dictionary ang time step ($\Delta t$) at ang bilang ng Trotter step, at ang kabuuang evolution time ay ang produkto nila. Mas malalayong oras ang naaabot ng mas malalaking step sa nakapirmi na circuit depth pero mas malaki ang Trotter error, na lumalaki kasabay ng $U$ dahil mas kaunti ang nagko-commute na term ng Hamiltonian. Gumamit ang nailathalang hardware study ng $\Delta t = 0.15$ para sa $U/t_h = 8$ [2], at ang parehong halaga ang ginagamit ng tutorial na ito. Linear na lumalaki ang circuit depth sa bilang ng step at hindi ito nakadepende sa $L$.

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)"
)

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

Ihanda ang initial state

May dalawang spin orbital ang bawat site, kaya $2L$ na qubit ang ginagamit ng chain na may $L$ site. Ang initial state ay isang Fock occupation string na may isang character kada 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}$$

Binabasa ang string mula sa pinakamababang significant na site muna, kaya ang pinakakanang character ay site 0 at ang pisikal na site index ng character na $k$ na binibilang mula kaliwa ay $i = L - 1 - k$. Kaya ang string na udud ay nangangahulugang d ang site 0, u ang site 1, d ang site 2, at u ang site 3.

Susunod, ihanda ang half-filled na Néel state, na salitan ang u at d, kung saan pinalitan ng vacancy ang gitnang site. Ang vacancy ang gambalang susundan ang pagkalat.

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)

Vacancy at site 15
udududududududueudududududududu

Buuin ang mga observable

Ang bawat observable ay isang dictionary na kumakatawan sa weighted sum ng mga produkto ng occupation operator:

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

Ang bawat entry ng operators ay isang term na $T_{\alpha}$, na isinulat bilang listahan ng mga pares na [name, site] na pinaparami, at ang katumbas na entry ng coefficients ang timbang nito. Ang dalawang pangalan ay "Nup" para sa $n_{j,\uparrow}$ at "Ndn" para sa $n_{j,\downarrow}$. Ganito isinusulat ang isang one-point observable na $n_{3,\uparrow}$:

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

Ipi-print ng susunod na cell ang pinakamaliit na kumpletong set ng mga input na tinatanggap ng function, na kapaki-pakinabang bilang sanggunian sa hugis habang binabasa ang code na kasunod.

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}")

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

Tatlong grupo ng mga observable ang kailangan. Ang charge density na $n_i = n_{i\uparrow} + n_{i\downarrow}$ at ang local spin na $S^z_i = (n_{i\uparrow} - n_{i\downarrow})/2$ ay parehong two-term sum sa iisang mga operator, na magkaiba lang sa mga coefficient. Ang spin correlator na $S^z_i S^z_{i_*}$ sa pagitan ng site $i$ at ng vacancy site ay lumalawak sa apat na two-point term na may mga coefficient na $\pm 1/4$.

Nagko-commute ang mga occupation operator, kaya walang epekto ang ayos ng mga factor sa loob ng isang term, pero ang ayos ng listahan ng observables ang nagtatakda ng ayos ng mga ibinabalik na halaga. Nagre-record ang code ng index para sa bawat grupo habang binubuo ang listahan. Nilalaktawan ng correlator ang vacancy site, kung saan uulitin ng produkto ang isang 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]])

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: I-optimize ang problema para sa pag-execute sa quantum hardware

Walang transpilation step ang tutorial na ito. Ginagawa ng function ang mga sumusunod na hakbang sa loob [2], at hindi mo ito ma-configure mula sa notebook:

• Minamapa ang mga fermionic mode sa mga qubit sa ayos na pinili para panatilihing mababa ang Pauli weight sa heavy-hex lattice

• Gumagamit ng fermionic SWAP network para maging two-qubit rotation ang bawat hopping term

• Inilalagay ang chain sa isang landas na umiiwas sa mahihinang qubit

• Naglalapat ng dynamical decoupling at Pauli twirling habang nagko-compile

Ikaw ang pipili ng backend at ng bilang ng shot. Tumatanggap ang function ng pangalan ng backend bilang string, kaya pinipili ng tutorial na ito ang device na pinakamaluwag ang pila at sapat ang laki para sa chain, at ipinapasa ang pangalan nito.

Ang pinakamaliit na signal na kailangan mong makita ang nagtatakda ng bilang ng shot. Umaabot ang charge tracer sa mga amplitude na malapit sa $0.15$, pero ang connected spin tracer ay pumipinakto sa paligid ng $0.04$, halos apat na beses na mas maliit. Dahil dito, gumamit ang nailathalang pag-aaral ng 20,000 shot kada circuit [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")

ibm_boston, 20000 shots per circuit

Step 3: I-execute gamit ang Qiskit primitives

Susunod, i-submit ang simulation. Agad na nagbabalik ng job handle ang run sa halip na maghintay.

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())

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

Gamitin ang job.status() para subaybayan ang progreso at ang job.result() para maghintay hanggang matapos ang job. Ang result ay isang mapping na ang entry na observable_expectation_values ay naka-index muna ayon sa observable, sa ayos na ibinigay ang mga ito, at pagkatapos ay ayon sa oras. May step_count + 1 na time point dahil kasama ang halaga sa $t=0$.

result = job.result()

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

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

Step 4: I-post-process at ibalik ang result sa nais na classical na format

Buuin muli ang mga correlator

Hatiin muli ang flat array sa tatlong grupo gamit ang mga indices na naitala kanina, tapos buuin ang dalawang tracer correlator. Ang charge tracer ay ang pagbabago sa density kada site kumpara sa $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.$$

Ang spin tracer ay ang connected correlator sa pagitan ng bawat site at ng vacancy site na $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).$$

Parehong buo mula sa mga operator sa Fock-basis, kaya galing sila sa parehong mga sinukat na bitstring at hindi na kailangan ng dagdag na hardware execution para sa pangalawa.

Noong binuo mo ang mga observable, itinala ng code ang posisyon ng bawat charge (charge_indices), spin (spin_indices), at spin-correlation (correlator_indices) na operator sa listahan. Dahil ibinabalik ng simulator ang mga expectation value sa parehong pagkakasunod-sunod, magagamit mo ang mga posisyong ito para kunin ang kaukulang mga resulta. Para sa mga spin correlation, ginagamit din ng code ang correlator_sites para ilagay ang bawat resulta sa tamang posisyon sa lattice, at iniiwang NaN ang vacancy.

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}"
)

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

I-visualize ang mga wavefront

Kapag ipinlot mo ang bawat correlator bilang heatmap ng espasyo at oras, direktang makikita ang light cone. Pababa ang daloy ng oras, at ang pahalang na axis ay ang posisyon sa kahabaan ng 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()

Kunin ang mga velocity ng wavefront

Gawing mga numero ang mga heatmap sa pamamagitan ng pagsunod sa detection protocol mula sa nailathalang pag-aaral [2]. Sa bawat time slice, i-symmetrize ang spatial profile sa paligid ng vacancy, na nag-a-average palabas sa antisymmetric na bahagi ng hardware noise, i-normalize gamit ang pinakamalaking off-center na amplitude, at itala ang pinakalabas na site na lampas sa threshold. Pagkatapos, hinihiwa ang mga detection sa mga run ng magkakasunod na time step, at ang pinakamahabang run lang ang fi-nifit, para hindi mapasama sa resulta ang early-time transient.

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)

Nakadepende ang fitted velocity sa detection threshold at sa trimming quantile, at mas malaki ang epekto ng sensitivity na iyon kaysa sa shot noise. I-scan ang dalawa at iulat ang median kasama ang buong envelope, gaya ng ginagawa sa nailathalang analysis.

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}")

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()

Kapag ibinalik sa mga heatmap ang mga detection, makikita kung nasundan ng algorithm ang leading edge. Ang mga filled marker ay ang mga detection na napasama sa fit. Ang mga open marker ay inalis ng mga masking rule, dahil masyado silang malapit sa vacancy, masyadong malapit sa dulo ng chain, o huli na kaya ang front ay lumapad na at hindi na malinaw ang posisyon nito. Ang dashed line ay ang fit sa mga oras na ginamit nito, at ang slope nito ang velocity sa title ng bawat panel.

Minarkahan sa ilalim ng axis ang column ng vacancy. Na-clip doon ang charge dahil ang halaga nito ay nasa order na isa, at walang halaga ang spin doon, kaya kulay gray ang pagkaguhit.

Mas malawak ang bukas ng charge cone kaysa sa spin cone. Iniiwan ng dalawang disturbance ang parehong defect sa parehong sandali at naglalakbay palabas sa magkaibang bilis.

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()

Nauuna ang charge front sa spin front. Ang isang hole sa one-dimensional na Mott background ay hindi naglalakbay bilang isang electron. Sa halip, nahahati ito sa isang holon na may dalang charge at isang spinon na may dalang spin disturbance, at sa $U/t_h = 8$ ang spinon ang mas mabagal sa dalawa.

Nakadepende ang pagkakasunod-sunod sa interaction. Sa $U = 0$, ang model ay mga free fermion na may $v_c = 2$ at $v_s = 4$, kaya mas mabilis ang spin. Pinipigilan ng repulsion ang spin transport habang halos hindi ginagalaw ang charge transport, at pagdating ng $U/t_h = 8$, ang mga nailathalang halaga mula sa hardware ay $v_c = 1.78$ at $v_s = 1.46$ [2]. Kapag inulit ang run gamit ang ilang halaga ng $U$ at ipinlot ang $v_s / v_c$, makikita ang crossover.

Tandaan ang mga caveat na ito sa mga numero. Bahagyang mababa ang bias ng threshold-based detection, kaya inaasahan ang charge velocity na ilang porsiyentong mas mababa sa nailathalang halaga. Malapad ang spin envelope dahil mas malapit sa noise floor ang connected correlator, kaya ang median ng scan ang binabanggit sa halip na isang fit lang.

Mga Sanggunian

[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.