Dinamica spinului și a sarcinii în modelul Fermi-Hubbard

Estimare de utilizare: 5 minute pe un procesor Heron r3. (NOTĂ: Aceasta este doar o estimare. Timpul tău de execuție poate varia.)

Context

Modelul Fermi-Hubbard descrie electroni care sar pe o rețea, cu o penalizare energetică pentru ocuparea dublă. Pe un lanț unidimensional de $L$ situri cu frontiere deschise, hamiltonianul este definit astfel:

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

unde $c^\dagger_{i\sigma}$ creează un fermion de spin $\sigma$ pe situl $i$, iar $n_{i\sigma} = c^\dagger_{i\sigma} c_{i\sigma}$ îi numără. Cei trei termeni sunt amplitudinea de salt $t_{h,i}$ pe legătura dintre situl $i$ și $i+1$, interacțiunea locală $U_i$ și potențialul chimic $\mu_i$. Tutorialul măsoară energia în unități de $t_h$ și timpul în unități de $1/t_h$.

Dinamica în timp real a acestui model este greu de calculat clasic. Diagonalizarea exactă se oprește aproape de $L \lesssim 14$ situri, iar după un quench, entropia de încâlcire a unui solver cu rețele tensoriale crește liniar în timp, astfel că dimensiunea legăturii necesară pentru o acuratețe fixă crește exponențial.

Într-o dimensiune, excitațiile elementare sunt spinonii, care poartă spin și nu poartă sarcină, și holonii, care poartă sarcină și nu poartă spin. Cei doi se propagă cu viteze diferite, astfel că o perturbație localizată se împarte într-un front rapid și unul lent. Acest efect se numește separare spin-sarcină, iar acest tutorial o măsoară [1].

Acest tutorial folosește Qiskit Function Dynamics Simulator, care primește hamiltonianul, o stare inițială, o schemă Trotter și o listă de observabile și returnează valori medii. Gestionează intern maparea Jordan-Wigner, trotterizarea, aranjamentul qubiților și suprimarea erorilor, deci în acest tutorial nu ai circuite de construit sau de transpilat.

Cerințe

Înainte de a începe acest tutorial, asigură-te că ai instalate următoarele pachete:

• Qiskit SDK v2.3.0 sau mai nou (pip install qiskit)

• Qiskit Runtime v0.50.0 sau mai nou (pip install qiskit-ibm-runtime)

• Qiskit Functions Catalog v0.21.0 sau mai nou (pip install qiskit-ibm-catalog)

• NumPy v2.5.1 sau mai nou (pip install numpy)

• Matplotlib v3.10.3 sau mai nou (pip install matplotlib)

• Q-CTRL Visualizer v11.0.0 sau mai nou (pip install qctrl-visualizer)

Configurare

Mai întâi, autentifică-te folosind cheia API IBM Quantum. Apoi selectează Qiskit Function după cum urmează. (Acest cod presupune că ți-ai salvat deja contul în mediul local.)

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

Pasul 1: Mapează intrările clasice la o problemă cuantică

Definește hamiltonianul

Funcția primește modelul ca dicționar cu trei chei. Fiecare acceptă fie un singur float, care se aplică uniform pe lanț, fie o listă care specifică valoarea sit cu sit. Listele au $L-1$ intrări pentru hopping_amplitude, deoarece un lanț deschis are cu o legătură mai puțin decât situri, și $L$ intrări pentru on_site_interaction și chemical_potential.

Acest tutorial folosește un lanț uniform de 31 de situri cu o interacțiune repulsivă $U/t_h = 8$. Un potențial chimic uniform contribuie doar cu o fază globală când numărul de particule este fix, așa că acest exemplu îl setează la zero. Potențialul chimic contează doar când variază de la un sit la altul, de exemplu pentru a adăuga dezordine sau o singură impuritate.

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

Stabilește schema Trotter

Dicționarul simulation conține pasul de timp ($\Delta t$) și numărul de pași Trotter, iar timpul total de evoluție este produsul lor. Pașii mai mari ajung la momente de timp mai târzii la adâncime fixă a circuitului, dar produc o eroare Trotter mai mare, care crește cu $U$ deoarece termenii hamiltonianului comută mai puțin. Studiul hardware publicat a folosit $\Delta t = 0.15$ pentru $U/t_h = 8$ [2], iar acest tutorial folosește aceeași valoare. Adâncimea circuitului crește liniar cu numărul de pași și nu depinde deloc de $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)

Pregătește starea inițială

Fiecare sit conține doi orbitali de spin, deci un lanț de $L$ situri folosește $2L$ qubiți. Starea inițială este un șir de ocupare Fock cu un caracter pe sit:

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

Șirul se citește începând cu situl cel mai puțin semnificativ, deci caracterul din dreapta este situl 0, iar indicele fizic al sitului pentru caracterul $k$ numărat din stânga este $i = L - 1 - k$. Prin urmare, șirul udud înseamnă că situl 0 este d, situl 1 este u, situl 2 este d, iar situl 3 este u.

Apoi, pregătește o stare Néel la semiumplere, alternând u și d, cu situl central înlocuit de o vacanță. Vacanța este perturbația a cărei propagare este urmărită.

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

Construiește observabilele

Fiecare observabilă este un dicționar care reprezintă o sumă ponderată de produse de operatori de ocupare:

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

Fiecare intrare din operators este un termen $T_{\alpha}$, scris ca o listă de perechi [name, site] care se înmulțesc, iar intrarea corespunzătoare din coefficients este ponderea sa. Cele două nume sunt "Nup" pentru $n_{j,\uparrow}$ și "Ndn" pentru $n_{j,\downarrow}$. O singură observabilă cu un punct $n_{3,\uparrow}$ se scrie astfel:

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

Celula următoare afișează cel mai mic set complet de intrări pe care îl acceptă funcția, util ca referință de formă în timp ce citești codul care urmează.

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

Sunt necesare trei grupuri de observabile. Densitatea de sarcină $n_i = n_{i\uparrow} + n_{i\downarrow}$ și spinul local $S^z_i = (n_{i\uparrow} - n_{i\downarrow})/2$ sunt ambele sume cu doi termeni peste aceiași operatori, diferind doar prin coeficienți. Corelatorul de spin $S^z_i S^z_{i_*}$ dintre situl $i$ și situl vacanței se extinde în patru termeni cu două puncte, cu coeficienți $\pm 1/4$.

Operatorii de ocupare comută, deci ordinea factorilor dintr-un termen este irelevantă, dar ordinea listei observables fixează ordinea valorilor returnate. Codul înregistrează un indice pentru fiecare grup pe măsură ce construiește lista. Corelatorul omite situl vacanței, unde produsul ar repeta un 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]}

Pasul 2: Optimizează problema pentru execuția pe hardware cuantic

Acest tutorial nu are un pas de transpilare. Funcția efectuează intern următorii pași [2], iar tu nu îi poți configura din notebook:

• Mapează modurile fermionice pe qubiți, într-o ordonare aleasă pentru a menține greutățile Pauli mici pe o rețea heavy-hex

• Folosește o rețea fermionică de SWAP-uri, astfel încât fiecare termen de salt devine o rotație pe doi qubiți

• Încorporează lanțul de-a lungul unui traseu care evită qubiții slabi

• Aplică decuplare dinamică și Pauli twirling în timpul compilării

Tu alegi backendul și numărul de shot-uri. Funcția primește numele backendului ca șir de caractere, deci acest tutorial selectează dispozitivul cel mai puțin ocupat, suficient de mare pentru lanț, și transmite numele acestuia.

Cel mai mic semnal pe care trebuie să îl rezolvi stabilește numărul de shot-uri. Trasorul de sarcină atinge amplitudini apropiate de $0.15$, dar trasorul de spin conectat atinge un maxim în jur de $0.04$, de aproximativ patru ori mai mic. Studiul publicat a folosit 20.000 de shot-uri per circuit din acest motiv [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

Pasul 3: Execută folosind primitivele Qiskit

Apoi, trimite simularea. run returnează imediat un handle al jobului, fără a bloca execuția.

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

Folosește job.status() pentru a urmări progresul și job.result() pentru a bloca execuția până când jobul se termină. Rezultatul este o mapare a cărei intrare observable_expectation_values este indexată mai întâi după observabilă, în ordinea în care au fost furnizate, apoi după timp. Există step_count + 1 puncte de timp deoarece valoarea de la $t=0$ este inclusă.

result = job.result()

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

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

Pasul 4: Postprocesează și returnează rezultatul în formatul clasic dorit

Reconstruiește corelatorii

Împarte tabloul plat înapoi în cele trei grupuri folosind indicii înregistrați anterior, apoi formează cei doi corelatori trasori. Trasorul de sarcină este variația densității pe site față de $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.$$

Trasorul de spin este corelatorul conectat dintre fiecare site și site-ul vacanței $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).$$

Ambii sunt construiți din operatori în baza Fock, deci provin din aceleași șiruri de biți măsurate, iar al doilea nu costă nicio execuție suplimentară pe hardware.

Când ai construit observabilele, codul a înregistrat poziția fiecărui operator de sarcină (charge_indices), de spin (spin_indices) și de corelație de spin (correlator_indices) în listă. Deoarece simulatorul returnează valorile medii în aceeași ordine, poți folosi aceste poziții pentru a extrage rezultatele corespunzătoare. Pentru corelațiile de spin, codul folosește și correlator_sites pentru a plasa fiecare rezultat în poziția corectă a rețelei și lasă vacanța ca 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}"
)

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

Vizualizează fronturile de undă

Dacă reprezinți fiecare corelator ca hartă termică a spațiului și timpului, conul de lumină este vizibil direct. Timpul curge în jos, iar axa orizontală este poziția de-a lungul lanțului.

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

Extrage vitezele fronturilor de undă

Transformă hărțile termice în numere urmând protocolul de detecție din studiul publicat [2]. La fiecare felie de timp, simetrizează profilul spațial în jurul vacanței, ceea ce elimină prin mediere partea antisimetrică a zgomotului hardware, normalizează după cea mai mare amplitudine din afara centrului și înregistrează cel mai exterior site aflat peste un prag. Detecțiile sunt apoi tăiate în secvențe de pași de timp consecutivi, iar doar cea mai lungă secvență este ajustată, ceea ce ține tranzitoriul de la timpuri scurte în afara rezultatului.

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)

Viteza ajustată depinde de pragul de detecție și de cuantila de tăiere, iar această sensibilitate domină zgomotul de eșantionare. Parcurge ambele și raportează mediana împreună cu întreaga înfășurătoare, așa cum face analiza publicată.

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

Redesenarea detecțiilor peste hărțile termice arată dacă algoritmul a urmărit marginea din față. Markerii plini sunt detecțiile care au intrat în ajustare. Markerii goi au fost excluși de regulile de mascare, fie pentru că sunt prea aproape de vacanță, prea aproape de capătul lanțului, fie pentru că sunt suficient de târzii încât frontul s-a lărgit și poziția lui nu mai este precisă. Linia întreruptă este ajustarea pe intervalul de timp folosit, iar panta ei este viteza din titlul fiecărui panou.

Coloana vacanței este marcată sub axă. Sarcina este tăiată acolo pentru că valoarea ei este de ordinul unu, iar spinul nu are valoare acolo, deci este desenat în gri.

Conul sarcinii se deschide mai larg decât conul spinului. Ambele perturbații pleacă din același defect în același moment și se propagă spre exterior cu viteze diferite.

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

Frontul sarcinii îl depășește pe cel al spinului. O gaură într-un fundal Mott unidimensional nu se propagă ca un electron. În schimb, se fracționalizează într-un holon care poartă sarcina și un spinon care poartă perturbația de spin, iar la $U/t_h = 8$ spinonul este cel mai lent dintre cei doi.

Ordinea depinde de interacțiune. La $U = 0$, modelul este de fermioni liberi cu $v_c = 2$ și $v_s = 4$, deci spinul este mai rapid. Repulsia suprimă transportul de spin și lasă transportul de sarcină aproape neatins, iar la $U/t_h = 8$ valorile publicate pe hardware sunt $v_c = 1.78$ și $v_s = 1.46$ [2]. Rularea repetată cu mai multe valori ale lui $U$ și reprezentarea grafică a lui $v_s / v_c$ cartografiază tranziția.

Ține cont de aceste rezerve privind valorile. Detecția bazată pe prag este ușor părtinitoare spre valori mici, deci este de așteptat o viteză a sarcinii cu câteva procente sub valoarea publicată. Înfășurătoarea spinului este largă deoarece corelatorul conectat se află mai aproape de pragul de zgomot, motiv pentru care se citează mediana scanării, nu o singură ajustare.

Referințe

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