Спінова та зарядова динаміка в моделі Фермі-Габбарда

Оцінка використання: 5 хвилин на процесорі Heron r3. (ПРИМІТКА: це лише оцінка. Фактичний час виконання може відрізнятися.)

Передісторія

Модель Фермі-Габбарда описує електрони, що стрибають по ґратці, з енергетичним штрафом за подвійне заповнення. На одновимірному ланцюжку з $L$ вузлів з відкритими межами гамільтоніан визначається так:

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

де $c^\dagger_{i\sigma}$ створює ферміон зі спіном $\sigma$ на вузлі $i$, а $n_{i\sigma} = c^\dagger_{i\sigma} c_{i\sigma}$ рахує їх. Три доданки — це амплітуда стрибка $t_{h,i}$ на зв'язку між вузлами $i$ та $i+1$, локальна взаємодія $U_i$ і хімічний потенціал $\mu_i$. Посібник вимірює енергію в одиницях $t_h$, а час — в одиницях $1/t_h$.

Динаміку цієї моделі в реальному часі важко обчислювати класично. Точна діагоналізація зупиняється приблизно на $L \lesssim 14$ вузлах, а після квенчу ентропія заплутаності розв'язувача на тензорних мережах зростає лінійно за часом, тому розмірність зв'язку, потрібна для фіксованої точності, зростає експоненційно.

В одному вимірі елементарними збудженнями є спіноні, які несуть спін і не несуть заряду, та холони, які несуть заряд і не несуть спіну. Вони поширюються з різними швидкостями, тому локалізоване збурення розщеплюється на швидкий і повільний фронти. Цей ефект називається розділенням спіну та заряду, і цей посібник вимірює його [1].

Цей посібник використовує Qiskit Function Dynamics Simulator, яка приймає гамільтоніан, початковий стан, розклад Троттера та список спостережуваних і повертає середні значення. Вона самостійно виконує відображення Йордана–Вігнера, тротеризацію, розташування кубітів і придушення помилок, тому в цьому посібнику не потрібно будувати чи транспілювати жодних схем.

Вимоги

Перш ніж почати цей навчальний посібник, переконайся, що встановлено такі пакети:

• Qiskit SDK v2.3.0 або новіший (pip install qiskit)

• Qiskit Runtime v0.50.0 або новіший (pip install qiskit-ibm-runtime)

• Каталог Qiskit Functions v0.21.0 або новіший (pip install qiskit-ibm-catalog)

• NumPy v2.5.1 або новіший (pip install numpy)

• Matplotlib v3.10.3 або новіший (pip install matplotlib)

• Q-CTRL Visualizer v11.0.0 або новіший (pip install qctrl-visualizer)

Налаштування

Спочатку автентифікуйся за допомогою свого ключа API IBM Quantum. Потім вибери Qiskit Function так. (Цей код припускає, що ти вже зберіг свій обліковий запис у локальному середовищі.)

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

Крок 1: Відобрази класичні входи на квантову задачу

Визначення гамільтоніана

Функція приймає модель як словник із трьома ключами. Кожен приймає або одне дійсне число, яке застосовується рівномірно вздовж ланцюжка, або список, що задає значення вузол за вузлом. Списки мають $L-1$ елементів для hopping_amplitude, оскільки відкритий ланцюжок має на один зв'язок менше, ніж вузлів, і $L$ елементів для on_site_interaction та chemical_potential.

Цей посібник використовує однорідний ланцюжок із 31 вузла з відштовхувальною взаємодією $U/t_h = 8$. Однорідний хімічний потенціал дає лише глобальну фазу, коли кількість частинок фіксована, тому в цьому прикладі його встановлено в нуль. Хімічний потенціал має значення лише тоді, коли він змінюється від вузла до вузла, наприклад, щоб додати безлад або одну домішку.

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

Задання розкладу Троттера

Словник simulation містить часовий крок ($\Delta t$) і кількість кроків Троттера, а загальний час еволюції є їхнім добутком. Більші кроки досягають пізніших моментів часу за фіксованої глибини схеми, але спричиняють більшу похибку Троттера, яка зростає з $U$, бо члени гамільтоніана менше комутують. Опубліковане апаратне дослідження використовувало $\Delta t = 0.15$ для $U/t_h = 8$ [2], і цей посібник використовує те саме значення. Глибина схеми зростає лінійно з кількістю кроків і зовсім не залежить від $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)

Підготовка початкового стану

Кожен вузол містить дві спінові орбіталі, тому ланцюжок із $L$ вузлів використовує $2L$ кубітів. Початковий стан — це рядок заповнення Фока з одним символом на вузол:

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

Рядок читається починаючи з наймолодшого вузла, тому крайній правий символ — це вузол 0, а фізичний індекс вузла символу $k$, відрахованого зліва, дорівнює $i = L - 1 - k$. Отже, рядок udud означає, що вузол 0 — це d, вузол 1 — u, вузол 2 — d, а вузол 3 — u.

Далі підготуй напівзаповнений стан Нееля, що чергує u та d, із центральним вузлом, заміненим на вакансію. Вакансія — це збурення, поширення якого відстежується.

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

Побудова спостережуваних

Кожна спостережувана — це словник, що подає зважену суму добутків операторів заповнення:

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

Кожен елемент operators — це один доданок $T_{\alpha}$, записаний як список пар [name, site], які перемножуються, а відповідний елемент coefficients є його вагою. Два імені: "Nup" для $n_{j,\uparrow}$ та "Ndn" для $n_{j,\downarrow}$. Одна одноточкова спостережувана $n_{3,\uparrow}$ записується так:

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

Наступна комірка виводить найменший повний набір вхідних даних, який приймає функція; він корисний як зразок структури під час читання коду далі.

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

Потрібні три групи спостережуваних. Густина заряду $n_i = n_{i\uparrow} + n_{i\downarrow}$ та локальний спін $S^z_i = (n_{i\uparrow} - n_{i\downarrow})/2$ — це двочленні суми за одними й тими самими операторами, що відрізняються лише коефіцієнтами. Спінова кореляція $S^z_i S^z_{i_*}$ між вузлом $i$ та вузлом вакансії розкладається на чотири двоточкові доданки з коефіцієнтами $\pm 1/4$.

Оператори заповнення комутують, тому порядок множників у доданку не має значення, але порядок списку observables визначає порядок повернених значень. Код записує індекс для кожної групи під час побудови списку. Корелятор пропускає вузол вакансії, де добуток повторював би множник.

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

Крок 2: Оптимізуй задачу для виконання на квантовому обладнанні

У цьому навчальному посібнику немає кроку транспіляції. Функція виконує внутрішньо такі кроки [2], і налаштувати їх з ноутбука не можна:

• Відображає ферміонні моди на кубіти в порядку, обраному так, щоб ваги Паулі лишалися низькими на гратці heavy-hex

• Використовує ферміонну мережу SWAP, щоб кожен доданок переходу став дворівневим обертанням двох кубітів

• Вбудовує ланцюжок вздовж шляху, що оминає слабкі кубіти

• Застосовує динамічне розв'язування та твірлінг Паулі під час компіляції

Ти обираєш бекенд і кількість шотів. Функція приймає назву бекенда як рядок, тому в цьому посібнику обирається найменш завантажений пристрій, достатньо великий для ланцюжка, і передається його назва.

Найменший сигнал, який потрібно розрізнити, визначає кількість шотів. Зарядовий трейсер сягає амплітуд близько $0.15$, а зв'язаний спіновий трейсер досягає піку близько $0.04$, тобто приблизно в чотири рази менше. З цієї причини в опублікованому дослідженні використовували 20 000 шотів на схему [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

Крок 3: Виконання за допомогою примітивів Qiskit

Далі надішли симуляцію. run одразу повертає дескриптор завдання, не блокуючи виконання.

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

Використовуй job.status(), щоб стежити за прогресом, і job.result(), щоб заблокувати виконання до завершення завдання. Результат — це відображення, чий елемент observable_expectation_values індексується спершу за спостережуваною, у порядку, в якому їх було надано, а потім за часом. Часових точок step_count + 1, оскільки значення при $t=0$ включене.

result = job.result()

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

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

Крок 4: Постобробка та повернення результату в потрібному класичному форматі

Відновлення корелярів

Розділи плоский масив назад на три групи, використовуючи раніше записані індекси, а потім сформуй два трейсерні корелятори. Зарядовий трейсер — це зміна густини на вузол відносно $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.$$

Спіновий трейсер — це зв'язаний корелятор між кожним вузлом і вузлом вакансії $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).$$

Обидва побудовані з операторів у базисі Фока, тож вони походять з одних і тих самих виміряних бітових рядків, і другий не потребує додаткового виконання на обладнанні.

Під час побудови спостережуваних код записав позицію кожного зарядового (charge_indices), спінового (spin_indices) та спін-кореляційного (correlator_indices) оператора у списку. Оскільки симулятор повертає математичні сподівання в тому самому порядку, ти можеш використати ці позиції, щоб витягти відповідні результати. Для спінових кореляцій код також використовує correlator_sites, щоб розмістити кожен результат у правильній позиції гратки, а вакансію лишає як 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

Візуалізація хвильових фронтів

Якщо побудувати кожен корелятор як теплову карту простору й часу, світловий конус стає добре видимим. Час тече вниз, а горизонтальна вісь — це позиція вздовж ланцюжка.

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

Витягання швидкостей хвильових фронтів

Перетвори теплові карти на числа за протоколом виявлення з опублікованого дослідження [2]. На кожному часовому зрізі симетризуй просторовий профіль відносно вакансії, що усереднює антисиметричну частину апаратного шуму, нормалізуй на найбільшу позацентральну амплітуду та запиши найдальший вузол над порогом. Потім виявлення розрізаються на серії послідовних часових кроків, і підганяється лише найдовша серія, що не пускає ранній перехідний процес у результат.

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)

Підігнана швидкість залежить від порога виявлення та квантиля обрізання, і ця чутливість переважає над шумом пострілів. Перебери обидва та повідом медіану з повною оболонкою, як це робить опублікований аналіз.

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

Накладання виявлень назад на теплові карти показує, чи алгоритм відстежив передній фронт. Заповнені маркери — це виявлення, що ввійшли в підгонку. Порожні маркери були відсіяні правилами маскування: або вони розташовані занадто близько до вакансії, занадто близько до кінця ланцюжка, або надто пізні, коли фронт вже розширився і його положення більше не чітке. Пунктирна лінія — це підгонка за часами, які вона використала, а її нахил — це швидкість у заголовку кожної панелі.

Стовпець вакансії позначено під віссю. Заряд там обрізається, бо його значення порядку одиниці, а спін там не має значення, тому його намальовано сірим.

Зарядовий конус розкривається ширше за спіновий. Обидва збурення залишають той самий дефект у той самий момент і розходяться назовні з різними швидкостями.

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

Зарядовий фронт випереджає спіновий. Дірка в одновимірному фоні Мотта не рухається як електрон. Натомість вона фракціонується на холон, що несе заряд, і спінон, що несе спінове збурення, а при $U/t_h = 8$ спінон повільніший із двох.

Порядок залежить від взаємодії. При $U = 0$ модель є вільними ферміонами з $v_c = 2$ та $v_s = 4$, тож спін швидший. Відштовхування пригнічує спіновий транспорт, майже не зачіпаючи зарядовий, і до $U/t_h = 8$ опубліковані апаратні значення становлять $v_c = 1.78$ та $v_s = 1.46$ [2]. Повторний запуск із кількома значеннями $U$ та побудова графіка $v_s / v_c$ відображає перехід.

Врахуй такі застереження щодо чисел. Виявлення за порогом дещо зміщене в бік заниження, тому очікується зарядова швидкість на кілька відсотків нижча за опубліковане значення. Спінова оболонка широка, бо зв'язаний корелятор лежить ближче до рівня шуму, тому наводиться медіана сканування, а не одна підгонка.

Посилання

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