Dinamik spin dan cas dalam model Fermi-Hubbard

Anggaran penggunaan: 5 minit pada pemproses Heron r3. (NOTA: Ini hanya anggaran. Masa larian anda mungkin berbeza.)

Latar belakang

Model Fermi-Hubbard menerangkan elektron yang melompat pada kekisi dengan penalti tenaga bagi penghunian berganda. Pada rantai satu dimensi dengan $L$ tapak dan sempadan terbuka, Hamiltonian ditakrifkan seperti 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},$$

di mana $c^\dagger_{i\sigma}$ mencipta fermion berspin $\sigma$ pada tapak $i$ dan $n_{i\sigma} = c^\dagger_{i\sigma} c_{i\sigma}$ membilangnya. Ketiga-tiga sebutan ialah amplitud lompatan $t_{h,i}$ pada ikatan antara tapak $i$ dan $i+1$, interaksi pada tapak $U_i$, dan keupayaan kimia $\mu_i$. Tutorial ini mengukur tenaga dalam unit $t_h$ dan masa dalam unit $1/t_h$.

Dinamik masa nyata model ini sukar dikira secara klasik. Pepenjuruan tepat terhenti pada kira-kira $L \lesssim 14$ tapak, dan selepas quench, entropi belitan penyelesai rangkaian tensor bertambah secara linear dengan masa, jadi dimensi ikatan yang diperlukan untuk ketepatan tetap bertambah secara eksponen.

Dalam satu dimensi, pengujaan asas ialah spinon, yang membawa spin tanpa cas, dan holon, yang membawa cas tanpa spin. Kedua-duanya merambat pada kelajuan berbeza, jadi gangguan setempat berpecah kepada muka hadapan yang pantas dan satu lagi yang perlahan. Kesan ini dipanggil pemisahan spin-cas, dan tutorial ini mengukurnya [1].

Tutorial ini menggunakan Qiskit Function Dynamics Simulator, yang menerima Hamiltonian, keadaan awal, jadual Trotter, dan senarai observable, serta mengembalikan nilai jangkaan. Ia mengendalikan pemetaan Jordan-Wigner, Trotterization, susun atur qubit, dan penindasan ralat secara dalaman, jadi tiada litar untuk dibina atau di-transpile dalam tutorial ini.

Keperluan

Sebelum anda memulakan tutorial ini, pastikan anda telah memasang pakej berikut:

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

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

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

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

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

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

Persediaan

Pertama, sahkan identiti menggunakan kunci API IBM Quantum anda. Kemudian, pilih Qiskit Function seperti berikut. (Kod ini mengandaikan anda telah menyimpan akaun anda ke persekitaran setempat anda.)

# 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 kepada masalah kuantum

Takrifkan Hamiltonian

Fungsi ini menerima model sebagai kamus dengan tiga kunci. Setiap satu menerima sama ada satu float tunggal, yang digunakan secara seragam sepanjang rantai, atau senarai yang menentukan nilai tapak demi tapak. Senarai mempunyai $L-1$ entri untuk hopping_amplitude, kerana rantai terbuka mempunyai satu ikatan kurang daripada bilangan tapak, dan $L$ entri untuk on_site_interaction dan chemical_potential.

Tutorial ini menggunakan rantai seragam dengan 31 tapak dan interaksi tolakan $U/t_h = 8$. Keupayaan kimia seragam hanya menyumbang fasa global apabila bilangan zarah tetap, jadi contoh ini menetapkannya kepada sifar. Keupayaan kimia hanya penting apabila ia berbeza dari tapak ke tapak, contohnya, untuk menambah gangguan atau satu kekotoran tunggal.

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

Tetapkan jadual Trotter

Kamus simulation memegang langkah masa ($\Delta t$) dan bilangan langkah Trotter, dan jumlah masa evolusi ialah hasil darab keduanya. Langkah yang lebih besar mencapai masa yang lebih lewat pada kedalaman litar tetap tetapi menimbulkan lebih banyak ralat Trotter, yang bertambah dengan $U$ kerana sebutan Hamiltonian kurang berkomut. Kajian perkakasan yang diterbitkan menggunakan $\Delta t = 0.15$ untuk $U/t_h = 8$ [2], dan tutorial ini menggunakan nilai yang sama. Kedalaman litar bertambah secara linear dengan bilangan langkah dan langsung 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)

Sediakan keadaan awal

Setiap tapak memegang dua orbital spin, jadi rantai $L$ tapak menggunakan $2L$ qubit. Keadaan awal ialah rentetan penghunian Fock dengan satu aksara bagi setiap tapak:

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

Rentetan dibaca bermula dengan tapak paling tidak signifikan, jadi aksara paling kanan ialah tapak 0 dan indeks tapak fizikal bagi aksara $k$ yang dikira dari kiri ialah $i = L - 1 - k$. Oleh itu rentetan udud bermaksud tapak 0 ialah d, tapak 1 ialah u, tapak 2 ialah d, dan tapak 3 ialah u.

Seterusnya, sediakan keadaan Néel separuh penuh, berselang-seli u dan d, dengan tapak tengah digantikan dengan kekosongan. Kekosongan ialah gangguan yang penyebarannya dijejaki.

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

Bina observable

Setiap observable ialah kamus yang mewakili hasil tambah berpemberat bagi hasil darab operator penghunian:

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

Setiap entri operators ialah satu sebutan $T_{\alpha}$, ditulis sebagai senarai pasangan [name, site] yang didarabkan bersama, dan entri coefficients yang sepadan ialah pemberatnya. Kedua-dua nama ialah "Nup" untuk $n_{j,\uparrow}$ dan "Ndn" untuk $n_{j,\downarrow}$. Satu observable satu titik tunggal $n_{3,\uparrow}$ ditulis seperti berikut:

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

Sel seterusnya mencetak set input lengkap terkecil yang diterima oleh fungsi, yang berguna sebagai rujukan bentuk semasa membaca kod yang berikut.

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

Tiga kumpulan observable diperlukan. Ketumpatan cas $n_i = n_{i\uparrow} + n_{i\downarrow}$ dan spin setempat $S^z_i = (n_{i\uparrow} - n_{i\downarrow})/2$ kedua-duanya ialah hasil tambah dua sebutan ke atas operator yang sama, hanya berbeza pada pekalinya. Korelator spin $S^z_i S^z_{i_*}$ antara tapak $i$ dan tapak kekosongan berkembang kepada empat sebutan dua titik dengan pekali $\pm 1/4$.

Operator penghunian berkomut, jadi susunan faktor dalam sesuatu sebutan tidak relevan, tetapi susunan senarai observables menetapkan susunan nilai yang dikembalikan. Kod merekodkan indeks bagi setiap kumpulan semasa ia membina senarai. Korelator melangkau tapak kekosongan, di mana hasil darab akan mengulangi 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: Optimumkan masalah untuk pelaksanaan pada perkakasan kuantum

Tutorial ini tidak mempunyai langkah transpilasi. Fungsi ini melakukan langkah berikut secara dalaman [2], dan anda tidak boleh mengkonfigurasinya daripada notebook:

• Memetakan mod fermion kepada qubit di bawah susunan yang dipilih untuk mengekalkan pemberat Pauli rendah pada kekisi heavy-hex

• Menggunakan rangkaian SWAP fermion supaya setiap sebutan lompatan menjadi putaran dua qubit

• Membenamkan rantai di sepanjang laluan yang mengelakkan qubit lemah

• Menggunakan dynamical decoupling dan Pauli twirling semasa kompilasi

Anda memilih backend dan bilangan shot. Fungsi ini menerima nama backend sebagai rentetan, jadi tutorial ini memilih peranti paling kurang sibuk yang cukup besar untuk rantai dan menghantar namanya.

Isyarat terkecil yang perlu anda leraikan menentukan bilangan shot. Pengesan cas mencapai amplitud sekitar $0.15$, tetapi pengesan spin bersambung memuncak sekitar $0.04$, kira-kira empat kali lebih kecil. Kajian yang diterbitkan menggunakan 20,000 shot bagi setiap litar atas sebab 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: Laksanakan dengan menggunakan primitif Qiskit

Seterusnya, hantar simulasi. run mengembalikan pemegang job serta-merta dan bukannya menyekat.

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

Gunakan job.status() untuk mengikuti kemajuan dan job.result() untuk menyekat sehingga job selesai. Hasilnya ialah pemetaan yang entri observable_expectation_valuesnya diindeks dahulu mengikut observable, mengikut susunan ia dibekalkan, dan kemudian mengikut masa. Terdapat step_count + 1 titik masa kerana nilai pada $t=0$ 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 pulangkan hasil dalam format klasik yang dikehendaki

Bina semula korelator

Pisahkan tatasusunan rata kembali kepada tiga kumpulan menggunakan indeks yang direkodkan sebelum ini, kemudian bentuk dua korelator pengesan. Pengesan cas ialah perubahan ketumpatan bagi setiap tapak relatif kepada $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.$$

Pengesan spin ialah korelator bersambung antara setiap tapak dan tapak 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).$$

Kedua-duanya dibina daripada operator asas Fock, jadi ia datang daripada rentetan bit terukur yang sama dan yang kedua tidak memerlukan pelaksanaan perkakasan tambahan.

Apabila anda membina observable, kod merekodkan kedudukan setiap operator cas (charge_indices), spin (spin_indices), dan korelasi spin (correlator_indices) dalam senarai. Oleh sebab simulator mengembalikan nilai jangkaan dalam susunan yang sama, anda boleh menggunakan kedudukan ini untuk mengekstrak hasil yang sepadan. Bagi korelasi spin, kod juga menggunakan correlator_sites untuk meletakkan setiap hasil pada kedudukan kekisi yang betul, 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

Visualkan muka gelombang

Jika anda memplot setiap korelator sebagai peta haba ruang dan masa, kon cahaya dapat dilihat secara langsung. Masa mengalir ke bawah, dan paksi mengufuk ialah kedudukan 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 halaju muka gelombang

Tukarkan peta haba kepada nombor dengan mengikuti protokol pengesanan daripada kajian yang diterbitkan [2]. Pada setiap hirisan masa, simetrikan profil ruang tentang kekosongan, yang memurnikan bahagian antisimetri hingar perkakasan, normalkan dengan amplitud luar pusat terbesar, dan rekodkan tapak paling luar yang melebihi ambang. Pengesanan kemudian dipotong kepada larian langkah masa berturutan, dan hanya larian terpanjang dipadankan, yang mengekalkan peralihan awal daripada 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)

Halaju yang dipadankan bergantung pada ambang pengesanan dan kuantil pemangkasan, dan kepekaan itu mengatasi hingar shot. Imbas kedua-duanya dan laporkan median dengan sampul penuh, seperti yang dilakukan oleh analisis yang diterbitkan.

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

Melukis semula pengesanan ke atas peta haba menunjukkan sama ada algoritma itu menjejaki pinggir hadapan. Penanda berisi ialah pengesanan yang digunakan dalam penyuaian. Penanda terbuka dipotong oleh peraturan penyembunyian, sama ada kerana terlalu dekat dengan kekosongan, terlalu dekat dengan hujung rantai, atau terlalu lewat sehingga bahagian hadapan telah melebar dan kedudukannya tidak lagi tajam. Garis putus-putus ialah penyuaian ke atas masa yang digunakan, dan kecerunannya ialah halaju dalam setiap tajuk panel.

Lajur kekosongan ditandakan di bawah paksi. Cas terpangkas di situ kerana nilainya berada pada susunan satu, dan spin tidak mempunyai nilai di situ, jadi ia dilukis dalam warna kelabu.

Kon cas terbuka lebih luas daripada kon spin. Kedua-dua gangguan meninggalkan kecacatan yang sama pada saat yang sama dan bergerak ke luar pada kelajuan yang berbeza.

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

Bahagian hadapan cas memintas bahagian hadapan spin. Lubang dalam latar Mott satu dimensi tidak bergerak sebagai elektron. Sebaliknya, ia berpecah kepada holon yang membawa cas dan spinon yang membawa gangguan spin, dan pada $U/t_h = 8$ spinon ialah yang lebih perlahan antara keduanya.

Susunan itu bergantung pada interaksi. Pada $U = 0$, model ialah fermion bebas dengan $v_c = 2$ dan $v_s = 4$, jadi spin lebih laju. Tolakan menyekat pengangkutan spin sambil hampir tidak menjejaskan pengangkutan cas, dan pada $U/t_h = 8$ nilai perkakasan yang diterbitkan ialah $v_c = 1.78$ dan $v_s = 1.46$ [2]. Menjalankan semula dengan beberapa nilai $U$ dan memplot $v_s / v_c$ memetakan peralihan tersebut.

Perhatikan beberapa kaveat pada nombor ini. Pengesanan berasaskan ambang berat sebelah sedikit ke bawah, jadi halaju cas yang beberapa peratus di bawah nilai yang diterbitkan adalah dijangka. Sampul spin adalah lebar kerana korelator bersambung lebih hampir dengan aras hingar, itulah sebabnya median imbasan dipetik dan bukannya satu penyuaian tunggal.

Rujukan

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