Spinová a nábojová dynamika v Fermiho-Hubbardově modelu

Odhad využití: 5 minut na procesoru Heron r3. (POZNÁMKA: Jde pouze o odhad. Tvoje doba běhu se může lišit.)

Pozadí

Fermiho-Hubbardův model popisuje elektrony poskakující po mřížce s energetickou penalizací za dvojí obsazení. Na jednorozměrném řetězci o $L$ místech s otevřenými okraji je Hamiltonián definován takto:

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

kde $c^\dagger_{i\sigma}$ vytváří fermion se spinem $\sigma$ na místě $i$ a $n_{i\sigma} = c^\dagger_{i\sigma} c_{i\sigma}$ je počítá. Tři členy jsou amplituda přeskoku $t_{h,i}$ na vazbě mezi místy $i$ a $i+1$, interakce na místě $U_i$ a chemický potenciál $\mu_i$. Tutoriál měří energii v jednotkách $t_h$ a čas v jednotkách $1/t_h$.

Dynamiku tohoto modelu v reálném čase je klasicky obtížné spočítat. Přesná diagonalizace končí zhruba u $L \lesssim 14$ míst a po kvenči entanglementová entropie řešiče založeného na tenzorových sítích roste lineárně s časem, takže vazební dimenze potřebná pro pevnou přesnost roste exponenciálně.

V jedné dimenzi jsou elementárními excitacemi spinony, které nesou spin a žádný náboj, a holony, které nesou náboj a žádný spin. Obojí se šíří různými rychlostmi, takže lokalizovaná porucha se rozdělí na rychlou a pomalou vlnoplochu. Tento jev se nazývá separace spinu a náboje a tento tutoriál ho měří [1].

Tento tutoriál používá Qiskit Function Dynamics Simulator, která přijímá Hamiltonián, počáteční stav, Trotterův rozvrh a seznam pozorovatelných a vrací střední hodnoty. Interně zajišťuje mapování Jordan-Wigner, trotterizaci, rozložení qubitů a potlačení chyb, takže v tomto tutoriálu není třeba sestavovat ani transpilovat žádné obvody.

Požadavky

Než začneš s tímto tutoriálem, ujisti se, že máš nainstalované následující balíčky:

• Qiskit SDK v2.3.0 nebo novější (pip install qiskit)

• Qiskit Runtime v0.50.0 nebo novější (pip install qiskit-ibm-runtime)

• Katalog Qiskit Functions v0.21.0 nebo novější (pip install qiskit-ibm-catalog)

• NumPy v2.5.1 nebo novější (pip install numpy)

• Matplotlib v3.10.3 nebo novější (pip install matplotlib)

• Q-CTRL Visualizer v11.0.0 nebo novější (pip install qctrl-visualizer)

Nastavení

Nejprve se ověř pomocí svého API klíče IBM Quantum. Poté vyber Qiskit Function následujícím způsobem. (Tento kód předpokládá, že jsi již uložil/a svůj účet do místního prostředí.)

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

Krok 1: Mapování klasických vstupů na kvantový problém

Definice Hamiltoniánu

Funkce přijímá model jako slovník se třemi klíči. Každý z nich přijímá buď jedno číslo typu float, které se aplikuje rovnoměrně na celý řetězec, nebo seznam, který určuje hodnotu pro každé místo zvlášť. Seznamy mají $L-1$ položek pro hopping_amplitude, protože otevřený řetězec má o jednu vazbu méně než míst, a $L$ položek pro on_site_interaction a chemical_potential.

Tento tutoriál používá rovnoměrný řetězec o 31 místech s repulzivní interakcí $U/t_h = 8$. Rovnoměrný chemický potenciál přispívá pouze globální fází, když je počet částic pevný, proto ho tento příklad nastavuje na nulu. Chemický potenciál má význam pouze tehdy, když se mění od místa k místu, například pro přidání poruchy nebo jedné příměsi.

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

Nastavení Trotterova rozvrhu

Slovník simulation obsahuje časový krok ($\Delta t$) a počet Trotterových kroků a celkový čas vývoje je jejich součin. Větší kroky dosáhnou pozdějších časů při pevné hloubce obvodu, ale způsobují větší Trotterovu chybu, která roste s $U$, protože členy Hamiltoniánu komutují méně. Publikovaná hardwarová studie použila $\Delta t = 0.15$ pro $U/t_h = 8$ [2] a tento tutoriál používá stejnou hodnotu. Hloubka obvodu roste lineárně s počtem kroků a vůbec nezávisí na $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)

Příprava počátečního stavu

Každé místo obsahuje dva spinové orbitaly, takže řetězec o $L$ místech používá $2L$ qubitů. Počáteční stav je Fockův obsazovací řetězec s jedním znakem na místo:

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

Řetězec se čte od nejméně významného místa, takže znak úplně vpravo je místo 0 a fyzický index místa pro znak $k$ počítaný zleva je $i = L - 1 - k$. Řetězec udud tedy znamená, že místo 0 je d, místo 1 je u, místo 2 je d a místo 3 je u.

Dále připrav napůl zaplněný Néelův stav se střídajícími se u a d, kde je centrální místo nahrazeno vakancí. Vakance je porucha, jejíž šíření se sleduje.

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

Sestavení pozorovatelných

Každá pozorovatelná je slovník, který představuje vážený součet součinů operátorů obsazení:

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

Každá položka operators je jeden člen $T_{\alpha}$, zapsaný jako seznam dvojic [name, site], které se násobí, a odpovídající položka coefficients je jeho váha. Dva názvy jsou "Nup" pro $n_{j,\uparrow}$ a "Ndn" pro $n_{j,\downarrow}$. Jedna jednobodová pozorovatelná $n_{3,\uparrow}$ se zapisuje takto:

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

Následující buňka vypíše nejmenší úplnou sadu vstupů, kterou funkce přijímá, což je užitečné jako vzor struktury při čtení následujícího kódu.

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

Potřebné jsou tři skupiny pozorovatelných. Hustota náboje $n_i = n_{i\uparrow} + n_{i\downarrow}$ a lokální spin $S^z_i = (n_{i\uparrow} - n_{i\downarrow})/2$ jsou obě dvoučlenné součty nad stejnými operátory, lišící se pouze koeficienty. Spinový korelátor $S^z_i S^z_{i_*}$ mezi místem $i$ a místem vakance se rozvine na čtyři dvoubodové členy s koeficienty $\pm 1/4$.

Operátory obsazení komutují, takže pořadí činitelů v rámci členu je irelevantní, ale pořadí seznamu observables určuje pořadí vrácených hodnot. Kód si při sestavování seznamu zaznamenává index pro každou skupinu. Korelátor přeskakuje místo vakance, kde by součin opakoval činitel.

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]}

Krok 2: Optimalizace problému pro spuštění na kvantovém hardwaru

Tento tutoriál nemá krok transpilace. Funkce interně provádí následující kroky [2] a z notebooku je nelze konfigurovat:

• Mapuje fermionické módy na qubity v pořadí zvoleném tak, aby Pauliho váhy zůstaly nízké na mřížce heavy-hex

• Používá fermionickou síť SWAP, takže každý člen přeskoku se stane dvouqubitovou rotací

• Vkládá řetězec podél cesty, která se vyhýbá slabým qubitům

• Během kompilace aplikuje dynamické oddělení a Pauliho twirling

Ty si vybíráš backend a počet shotů. Funkce přijímá název backendu jako řetězec, takže tento tutoriál vybere nejméně vytížené zařízení dostatečně velké pro řetězec a předá jeho název.

Počet shotů určuje nejmenší signál, který potřebuješ rozlišit. Nábojový tracer dosahuje amplitud kolem $0.15$, ale propojený spinový tracer vrcholí kolem $0.04$, tedy zhruba čtyřikrát méně. Publikovaná studie z tohoto důvodu použila 20 000 shotů na obvod [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

Krok 3: Spuštění pomocí primitiv Qiskit

Dále odešli simulaci. run okamžitě vrátí handle úlohy, aniž by blokoval.

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

Pomocí job.status() sleduj průběh a pomocí job.result() počkej, až úloha skončí. Výsledek je mapování, jehož položka observable_expectation_values je indexována nejprve podle observabilu, v pořadí, v jakém byly zadány, a poté podle času. Časových bodů je step_count + 1, protože hodnota při $t=0$ je zahrnuta.

result = job.result()

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

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

Krok 4: Následné zpracování a vrácení výsledku v požadovaném klasickém formátu

Rekonstrukce korelátorů

Rozděl ploché pole zpět do tří skupin pomocí dříve zaznamenaných indexů a poté sestav dva tracerové korelátory. Nábojový tracer je změna hustoty na jednotlivých místech vzhledem k $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.$$

Spinový tracer je propojený korelátor mezi každým místem a místem vakance $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).$$

Oba jsou sestaveny z operátorů ve Fockově bázi, takže pocházejí ze stejných změřených bitových řetězců a druhý z nich nevyžaduje žádné další spuštění na hardwaru.

Při sestavování observabilů kód zaznamenal pozici každého nábojového (charge_indices), spinového (spin_indices) a spin-korelačního (correlator_indices) operátoru v seznamu. Protože simulátor vrací očekávané hodnoty ve stejném pořadí, můžeš tyto pozice použít k extrakci odpovídajících výsledků. U spinových korelací kód také používá correlator_sites, aby každý výsledek umístil na správnou pozici v mřížce, a vakanci ponechá jako 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

Vizualizace čel vln

Pokud každý korelátor vykreslíš jako teplotní mapu prostoru a času, světelný kužel je přímo viditelný. Čas plyne směrem dolů a vodorovná osa je poloha podél řetězce.

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

Extrakce rychlostí čel vln

Přeměň teplotní mapy na čísla podle detekčního protokolu z publikované studie [2]. V každém časovém řezu symetrizuj prostorový profil kolem vakance, čímž se vyruší antisymetrická část hardwarového šumu, normalizuj největší amplitudou mimo střed a zaznamenej nejvzdálenější místo nad prahem. Detekce se poté rozdělí na úseky po sobě jdoucích časových kroků a fituje se pouze nejdelší úsek, což udrží počáteční přechodový jev mimo výsledek.

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)

Fitovaná rychlost závisí na detekčním prahu a kvantilu ořezu a tato citlivost převažuje nad shot šumem. Proskenuj obojí a uveď medián s úplnou obálkou, stejně jako to dělá publikovaná analýza.

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

Vykreslení detekcí zpět do teplotních map ukazuje, zda algoritmus sledoval náběžnou hranu. Plné značky jsou detekce, které vstoupily do fitu. Prázdné značky byly vyřazeny pravidly maskování, a to buď proto, že leží příliš blízko vakance, příliš blízko konce řetězce, nebo byly tak pozdě, že se čelo rozšířilo a jeho poloha už není ostrá. Přerušovaná čára je fit přes časy, které použil, a její sklon je rychlost uvedená v nadpisu každého panelu.

Sloupec vakance je označen pod osou. Náboj tam saturuje, protože jeho hodnota je řádu jedna, a spin tam nemá hodnotu, proto je vykreslen šedě.

Nábojový kužel se otevírá šířeji než spinový. Obě poruchy opouštějí stejný defekt ve stejném okamžiku a šíří se ven různými rychlostmi.

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

Nábojové čelo předbíhá spinové čelo. Díla v jednorozměrném Mottově pozadí necestuje jako elektron. Místo toho se frakcionalizuje na holon nesoucí náboj a spinon nesoucí spinovou poruchu a při $U/t_h = 8$ je spinon z těch dvou pomalejší.

Pořadí závisí na interakci. Při $U = 0$ je model volné fermiony s $v_c = 2$ a $v_s = 4$, takže spin je rychlejší. Repulze potlačuje transport spinu a transport náboje ponechává téměř nedotčený a při $U/t_h = 8$ jsou publikované hodnoty z hardwaru $v_c = 1.78$ a $v_s = 1.46$ [2]. Opakované spuštění s několika hodnotami $U$ a vykreslení $v_s / v_c$ zmapuje přechod.

Ke číslům si všimni těchto výhrad. Detekce založená na prahu je mírně zkreslená směrem dolů, takže se očekává nábojová rychlost o několik procent nižší než publikovaná hodnota. Spinová obálka je široká, protože propojený korelátor leží blíž ke šumovému pozadí, a proto se uvádí medián ze skenu místo jediného fitu.

Reference

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