Dinamika spin dan muatan dalam model Fermi-Hubbard

Perkiraan penggunaan: 5 menit pada prosesor Heron r3. (CATATAN: Ini hanya perkiraan. Waktu jalan aktualnya bisa berbeda.)

Latar belakang

Model Fermi-Hubbard menggambarkan elektron yang melompat pada sebuah kisi dengan penalti energi untuk okupansi ganda. Pada rantai satu dimensi dengan $L$ situs dan batas terbuka, Hamiltonian didefinisikan sebagai berikut:

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

dengan $c^\dagger_{i\sigma}$ menciptakan fermion berspin $\sigma$ pada situs $i$ dan $n_{i\sigma} = c^\dagger_{i\sigma} c_{i\sigma}$ menghitungnya. Ketiga sukunya adalah amplitudo hopping $t_{h,i}$ pada ikatan antara situs $i$ dan $i+1$, interaksi on-site $U_i$, dan potensial kimia $\mu_i$. Tutorial ini mengukur energi dalam satuan $t_h$ dan waktu dalam satuan $1/t_h$.

Dinamika real-time model ini sulit dihitung secara klasik. Diagonalisasi eksak berhenti di sekitar $L \lesssim 14$ situs, dan setelah quench, entropi keterikatan pada solver tensor-network tumbuh secara linear terhadap waktu, sehingga dimensi ikatan yang dibutuhkan untuk akurasi tetap tumbuh secara eksponensial.

Dalam satu dimensi, eksitasi elementernya adalah spinon, yang membawa spin tanpa muatan, dan holon, yang membawa muatan tanpa spin. Keduanya merambat dengan kecepatan berbeda, sehingga gangguan yang terlokalisasi terpecah menjadi front yang cepat dan front yang lambat. Efek ini disebut pemisahan spin-muatan, dan tutorial ini mengukurnya [1].

Tutorial ini memakai Qiskit Function Dynamics Simulator, yang menerima Hamiltonian, keadaan awal, jadwal Trotter, dan daftar observabel, lalu mengembalikan nilai ekspektasi. Fungsi ini menangani pemetaan Jordan-Wigner, Trotterisasi, tata letak qubit, dan penekanan error secara internal, jadi tidak ada circuit yang perlu dibangun atau di-transpile di tutorial ini.

Persyaratan

Sebelum memulai tutorial ini, pastikan kamu sudah menginstal paket-paket berikut:

• Qiskit SDK v2.3.0 atau lebih baru (pip install qiskit)

• Qiskit Runtime v0.50.0 atau lebih baru (pip install qiskit-ibm-runtime)

• Qiskit Functions Catalog v0.21.0 atau lebih baru (pip install qiskit-ibm-catalog)

• NumPy v2.5.1 atau lebih baru (pip install numpy)

• Matplotlib v3.10.3 atau lebih baru (pip install matplotlib)

• Q-CTRL Visualizer v11.0.0 atau lebih baru (pip install qctrl-visualizer)

Penyiapan

Pertama, lakukan autentikasi dengan IBM Quantum API key kamu. Lalu, pilih Qiskit Function seperti berikut. (Kode ini mengasumsikan kamu sudah menyimpan akunmu di lingkungan lokal.)

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

Langkah 1: Petakan input klasik ke masalah kuantum

Definisikan Hamiltonian

Fungsi ini menerima model sebagai dictionary dengan tiga key. Masing-masing menerima satu float tunggal, yang berlaku seragam di sepanjang rantai, atau sebuah list yang menentukan nilai per situs. List berisi $L-1$ entri untuk hopping_amplitude, karena rantai terbuka punya satu ikatan lebih sedikit daripada jumlah situs, dan $L$ entri untuk on_site_interaction dan chemical_potential.

Tutorial ini memakai rantai seragam dengan 31 situs dan interaksi tolak-menolak $U/t_h = 8$. Potensial kimia yang seragam hanya menyumbang fase global ketika jumlah partikel tetap, jadi contoh ini menetapkannya ke nol. Potensial kimia baru berpengaruh ketika nilainya berbeda dari situs ke situs, misalnya untuk menambahkan disorder atau satu pengotor.

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

Atur jadwal Trotter

Dictionary simulation berisi langkah waktu ($\Delta t$) dan jumlah langkah Trotter, dan total waktu evolusi adalah hasil kali keduanya. Langkah yang lebih besar mencapai waktu yang lebih lama pada kedalaman circuit tetap, tetapi menimbulkan lebih banyak error Trotter, yang tumbuh seiring $U$ karena suku-suku Hamiltonian makin kurang komut. Studi hardware yang dipublikasikan memakai $\Delta t = 0.15$ untuk $U/t_h = 8$ [2], dan tutorial ini memakai nilai yang sama. Kedalaman circuit tumbuh secara linear terhadap jumlah langkah dan sama sekali tidak bergantung pada $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)

Siapkan keadaan awal

Setiap situs menampung dua orbital spin, jadi rantai dengan $L$ situs memakai $2L$ qubit. Keadaan awalnya adalah string okupansi Fock dengan satu karakter per situs:

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

String dibaca dengan situs paling tidak signifikan lebih dulu, jadi karakter paling kanan adalah situs 0 dan indeks situs fisik dari karakter $k$ yang dihitung dari kiri adalah $i = L - 1 - k$. Jadi string udud berarti situs 0 adalah d, situs 1 adalah u, situs 2 adalah d, dan situs 3 adalah u.

Selanjutnya, siapkan keadaan Néel half-filled, dengan u dan d bergantian, dan situs tengah diganti dengan kekosongan (vacancy). Kekosongan ini adalah gangguan yang penyebarannya dilacak.

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

Bangun observabel

Setiap observabel adalah dictionary yang merepresentasikan jumlah berbobot dari hasil kali operator okupansi:

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

Setiap entri operators adalah satu suku $T_{\alpha}$, ditulis sebagai list pasangan [name, site] yang dikalikan, dan entri coefficients yang bersesuaian adalah bobotnya. Dua nama yang dipakai adalah "Nup" untuk $n_{j,\uparrow}$ dan "Ndn" untuk $n_{j,\downarrow}$. Satu observabel satu-titik $n_{3,\uparrow}$ ditulis sebagai berikut:

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

Sel berikutnya mencetak himpunan input lengkap terkecil yang diterima fungsi, yang berguna sebagai acuan bentuk saat membaca kode selanjutnya.

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

Dibutuhkan tiga kelompok observabel. Kerapatan muatan $n_i = n_{i\uparrow} + n_{i\downarrow}$ dan spin lokal $S^z_i = (n_{i\uparrow} - n_{i\downarrow})/2$ keduanya adalah jumlah dua suku dari operator yang sama, hanya berbeda pada koefisiennya. Korelator spin $S^z_i S^z_{i_*}$ antara situs $i$ dan situs kekosongan diuraikan menjadi empat suku dua-titik dengan koefisien $\pm 1/4$.

Operator okupansi saling komut, jadi urutan faktor dalam satu suku tidak berpengaruh, tetapi urutan list observables menentukan urutan nilai yang dikembalikan. Kode mencatat sebuah indeks untuk tiap kelompok saat membangun list. Korelator melewati situs kekosongan, karena di sana hasil kali akan mengulang satu faktor.

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

Langkah 2: Optimalkan masalah untuk eksekusi di hardware kuantum

Tutorial ini tidak punya langkah transpilasi. Fungsi ini melakukan langkah-langkah berikut secara internal [2], dan kamu tidak bisa mengonfigurasinya dari notebook:

• Memetakan mode fermionik ke qubit dengan urutan yang dipilih agar bobot Pauli tetap rendah pada kisi heavy-hex

• Memakai jaringan SWAP fermionik sehingga setiap suku hopping menjadi rotasi dua-qubit

• Menanamkan rantai di sepanjang jalur yang menghindari qubit lemah

• Menerapkan dynamical decoupling dan Pauli twirling selama kompilasi

Kamu memilih backend dan jumlah shot. Fungsi ini menerima nama backend sebagai string, jadi tutorial ini memilih perangkat paling sepi antrean yang cukup besar untuk rantai tersebut dan mengirimkan namanya.

Sinyal terkecil yang perlu kamu resolusikan menentukan jumlah shot. Tracer muatan mencapai amplitudo sekitar $0.15$, tetapi tracer spin terhubung memuncak sekitar $0.04$, kira-kira empat kali lebih kecil. Studi yang dipublikasikan memakai 20.000 shot per circuit karena alasan ini [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

Langkah 3: Eksekusi dengan Qiskit primitives

Selanjutnya, kirim simulasinya. run langsung mengembalikan job handle tanpa menunggu.

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

Pakai job.status() untuk mengikuti progres dan job.result() untuk menunggu sampai job selesai. Hasilnya adalah mapping yang entri observable_expectation_values-nya diindeks pertama menurut observabel, sesuai urutan saat diberikan, lalu menurut waktu. Ada step_count + 1 titik waktu karena nilai pada $t=0$ ikut disertakan.

result = job.result()

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

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

Langkah 4: Pasca-proses dan kembalikan hasil dalam format klasik yang diinginkan

Rekonstruksi korelator

Pisahkan array datar kembali menjadi tiga kelompok memakai indeks yang dicatat sebelumnya, lalu bentuk dua korelator tracer. Tracer muatan adalah perubahan kerapatan per situs relatif terhadap $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.$$

Tracer spin adalah korelator terhubung antara tiap situs dan situs kekosongan $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).$$

Keduanya dibangun dari operator basis Fock, jadi berasal dari bitstring terukur yang sama dan yang kedua tidak memakan eksekusi hardware tambahan.

Saat kamu membangun observabel, kode mencatat posisi tiap operator muatan (charge_indices), spin (spin_indices), dan korelasi spin (correlator_indices) di dalam list. Karena simulator mengembalikan nilai ekspektasi dengan urutan yang sama, kamu bisa memakai posisi ini untuk mengambil hasil yang sesuai. Untuk korelasi spin, kode juga memakai correlator_sites untuk menempatkan tiap hasil pada posisi kisi yang benar, dan membiarkan kekosongan sebagai 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

Visualisasikan front gelombang

Kalau kamu memplot tiap korelator sebagai heatmap ruang dan waktu, light cone langsung terlihat. Waktu mengalir ke bawah, dan sumbu horizontal adalah posisi di sepanjang rantai.

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

Ekstrak kecepatan front gelombang

Ubah heatmap menjadi angka dengan mengikuti protokol deteksi dari studi yang dipublikasikan [2]. Pada tiap irisan waktu, simetrikan profil spasial terhadap kekosongan, yang merata-ratakan bagian antisimetris dari noise hardware, normalisasi dengan amplitudo terbesar di luar pusat, dan catat situs terluar yang berada di atas ambang. Deteksi lalu dipotong menjadi run dari langkah waktu berurutan, dan hanya run terpanjang yang di-fit, sehingga transien waktu awal tidak ikut memengaruhi hasil.

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)

Kecepatan hasil fit bergantung pada ambang deteksi dan kuantil pemangkasan, dan sensitivitas itu mendominasi shot noise. Pindai keduanya dan laporkan median beserta envelope penuhnya, seperti yang dilakukan analisis yang dipublikasikan.

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

Menggambar balik deteksi ke heatmap menunjukkan apakah algoritma melacak tepi terdepan. Penanda terisi adalah deteksi yang masuk ke fit. Penanda kosong dibuang oleh aturan masking, entah karena posisinya terlalu dekat dengan kekosongan, terlalu dekat dengan ujung rantai, atau cukup terlambat sehingga front sudah melebar dan posisinya tidak lagi tajam. Garis putus-putus adalah fit pada waktu yang dipakainya, dan kemiringannya adalah kecepatan pada judul tiap panel.

Kolom vakansi ditandai di bawah sumbu. Muatan terpotong di sana karena nilainya berorde satu, sedangkan spin tidak punya nilai di sana, jadi digambar dengan warna abu-abu.

Kerucut muatan terbuka lebih lebar daripada kerucut spin. Kedua gangguan meninggalkan cacat yang sama pada saat yang sama dan merambat keluar dengan kecepatan berbeda.

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

Front muatan mendahului front spin. Sebuah hole dalam latar Mott satu dimensi tidak merambat sebagai elektron. Sebaliknya, ia terpecah menjadi holon yang membawa muatan dan spinon yang membawa gangguan spin, dan pada $U/t_h = 8$ spinon adalah yang lebih lambat dari keduanya.

Urutannya bergantung pada interaksi. Pada $U = 0$, modelnya adalah fermion bebas dengan $v_c = 2$ dan $v_s = 4$, jadi spin lebih cepat. Tolakan menekan transpor spin sementara transpor muatan nyaris tidak terpengaruh, dan pada $U/t_h = 8$ nilai perangkat keras yang dipublikasikan adalah $v_c = 1.78$ dan $v_s = 1.46$ [2]. Menjalankan ulang dengan beberapa nilai $U$ dan memplot $v_s / v_c$ akan memetakan crossover-nya.

Perhatikan beberapa catatan tentang angka-angka ini. Deteksi berbasis ambang cenderung sedikit bias ke bawah, jadi kecepatan muatan yang beberapa persen di bawah nilai yang dipublikasikan itu wajar. Envelope spin lebar karena korelator terhubung lebih dekat ke noise floor, itulah sebabnya yang dikutip adalah median dari hasil scan, bukan satu fit saja.

Referensi

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