ديناميكيات السبين والشحنة في نموذج فيرمي-هابارد

تقدير الاستخدام: 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$ موقعًا، وبعد الإخماد (quench) تنمو إنتروبيا التشابك لحلّال الشبكات التنسورية خطيًا مع الزمن، فينمو بُعد الرابطة اللازم لدقة ثابتة أُسّيًا.

في بُعد واحد، الإثارات الأولية هي السبينونات (spinons) التي تحمل سبينًا دون شحنة، والهولونات (holons) التي تحمل شحنة دون سبين. وينتشر النوعان بسرعتين مختلفتين، فينقسم الاضطراب الموضعي إلى جبهة سريعة وأخرى بطيئة. يُسمّى هذا التأثير فصل السبين والشحنة، ويقيسه هذا الدرس [1].

يستخدم هذا الدرس Qiskit Function المسماة Dynamics Simulator، والتي تأخذ الهاميلتوني وحالة ابتدائية وجدول Trotter وقائمة من الملاحظات، وتُرجع القيم المتوقعة. وهي تتولى داخليًا ربط Jordan-Wigner وتقسيم Trotter وتخطيط الكيوبتات وكبت الأخطاء، لذا لا توجد دوائر لبنائها أو تحويلها في هذا الدرس.

المتطلبات

قبل أن تبدأ هذا الدرس، تأكد من تثبيت الحزم التالية:

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

الإعداد

أولًا، سجّل الدخول باستخدام مفتاح IBM Quantum API. ثم اختر 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

ضبط جدول Trotter

يحتوي القاموس simulation على الخطوة الزمنية ($\Delta t$) وعدد خطوات Trotter، والزمن الكلي للتطور هو حاصل ضربهما. تصل الخطوات الأكبر إلى أزمنة أبعد بعمق دائرة ثابت لكنها تسبب خطأ Trotter أكبر، وهو ينمو مع $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$ كيوبت. الحالة الابتدائية سلسلة إشغال Fock بحرف واحد لكل موقع:

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

بعد ذلك، حضّر حالة Néel نصف ممتلئة، تتناوب فيها 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: تحسين المسألة للتنفيذ على العتاد الكمّي

لا توجد في هذا الدرس خطوة تحويل (transpilation). تنفّذ الدالة الخطوات التالية داخليًا [2]، ولا يمكنك ضبطها من دفتر الملاحظات:

• تربط الأنماط الفرميونية بالكيوبتات وفق ترتيب مختار للإبقاء على أوزان Pauli منخفضة على شبكة heavy-hex

• تستخدم شبكة SWAP فرميونية بحيث يصبح كل حد قفز دورانًا على كيوبتين

• تُدمج السلسلة على مسار يتجنب الكيوبتات الضعيفة

• تطبّق الفصل الديناميكي (dynamical decoupling) وتدوير Pauli (Pauli twirling) أثناء التجميع

أنت تختار الـBackend وعدد اللقطات. تأخذ الدالة اسم Backend كسلسلة نصية، لذا يختار هذا الدرس الجهاز الأقل انشغالًا والكبير بما يكفي للسلسلة ويمرّر اسمه.

أصغر إشارة تحتاج إلى تمييزها هي التي تحدد عدد اللقطات. يصل متتبّع الشحنة إلى سعات قرب $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 primitives

بعد ذلك، أرسل المحاكاة. تُرجع 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.