انتقل إلى المحتوى الرئيسي

رصد ديناميكيات هادرونية غير أبيلية قوية ومتماسكة على معالجات كمومية صاخبة

الاستخدام التقديري: 6 دقائق على معالج Heron (ibm_boston أو ما يعادله) (ملاحظة: هذا تقدير فحسب. قد يختلف وقت التشغيل الفعلي.)

مخرجات التعلّم​

  • كيف يمكن إعادة صياغة نظريات المعايرة الشبكية غير الأبيلية (وتحديدًا SU(2)) باستخدام إطار عمل الحلقة-السلسلة-الهادرون (LSH) من أجل محاكاة كمومية فعالة

  • كيفية بناء دارات تطور زمني من نوع Trotter لهاملتوني نظرية معايرة SU(2) تقريبية وتخطيطها على الكيوبتات

  • كيفية تشغيل هذه الدارات على أجهزة IBM Quantum® باستخدام أداة Qiskit الأساسية Estimator مع تخفيف أخطاء القراءة

المتطلبات الأساسية​

الخلفية​

الدافع​

الديناميكا اللونية الكمومية (QCD)، وهي نظرية معايرة SU(3) للقوة النووية الشديدة، تربط الكواركات في هادرونات وتحكم عمليتي الحجز (confinement) وكسر السلسلة. تتفوق طرق QCD الشبكية الكلاسيكية في وصف الخصائص الساكنة لكنها لا تستطيع محاكاة الديناميكيات في الزمن الحقيقي بسبب مشكلة الإشارة. توفر الحواسيب الكمومية طريقًا للالتفاف حول هذا الحاجز عبر ترميز درجات حرية حقل المعايرة مباشرة على الكيوبتات.

يوضح هذا الدرس التعليمي محاكاة من هذا القبيل: استخدام أجهزة IBM Quantum لمحاكاة انتشار الهادرونات في الزمن الحقيقي ضمن نظرية معايرة شبكية SU(2) ذات بعدين (1+1) — وهي أبسط نظرية معايرة غير أبيلية وخطوة تمهيدية نحو QCD الكاملة.

هاملتوني كوغوت-سسكيند​

تُصاغ النظرية على شبكة مكانية أحادية البعد تحتوي على فرميونات متعرجة (staggered) (المادة) عند المواقع وحقول معايرة SU(2) على الروابط. وبعد إعادة القياس إلى صورة عديمة الأبعاد، يكون الهاملتوني:

W=HE(KS)+μHM+xHI(KS),W = H_E^{\text{(KS)}} + \mu H_M + x H_I^{\text{(KS)}},

حيث HEH_E هي طاقة الحقل الكهربائي اللوني، وHMH_M حد الكتلة المتعرجة، وHIH_I حد تفاعل المادة والمعايرة (القفز)، وμ=2mgx\mu = 2\frac{m}{g}\sqrt{x} يرمّز كتلة الفرميون، وx=1g2a2x = \frac{1}{g^2 a^2} هو قوة التفاعل. يقع نهاية الوسط المتصل للنظرية عند N→∞N \to \infty وx→∞x \to \infty.

إطار عمل الحلقة-السلسلة-الهادرون (LSH)​

يتمثل أحد التحديات الرئيسية في أن فضاء هيلبرت لحقل المعايرة عند كل رابط ذو بُعد لانهائي. يعالج إطار عمل الحلقة-السلسلة-الهادرون (LSH) هذا الأمر عبر إعادة صياغة النظرية بدلالة متغيرات ثابتة تحت المعايرة — حلقات من التدفق، وسلاسل تصل بين شحنات منفصلة، وهادرونات (أزواج فرميونية أحادية تحت المعايرة عند موقع ما). في أساس LSH، يُستوفى قانون غاوس تلقائيًا بحكم البناء، لذا فإن كل حالة أساس فيزيائية. يتميز كل موقع في الشبكة بثلاثة أعداد كمومية (nl,ni,no)(n_l, n_i, n_o) تمثل عدد الحلقات، والسلسلة الداخلة، والسلسلة الخارجة، حيث ni,no∈{0,1}n_i, n_o \in \{0,1\} فرميونية وnl≥0n_l \geq 0 بوزونية. يُعرَّف عدد الفرميونات المحلي من هذه الكميات بأنه nf(r)=ni(r)+no(r)n_f(r) = n_i(r) + n_o(r) للمواقع الزوجية وnf(r)=2−[ni(r)+no(r)]n_f(r) = 2 - [n_i(r) + n_o(r)] للمواقع الفردية.

من الهاملتوني الكامل إلى الدارة الكمومية: ثلاثة تقريبات رئيسية​

لا تحاكي الدارة الكمومية هاملتوني SU(2) الكامل بدقة إطلاقًا. بل تنفّذ سلسلة تقريبات مضبوطة صالحة في نظام الاقتران الضعيف (x≫1x \gg 1). من الضروري فهم ما يُقرَّب وما لا يُقرَّب:

التقريب 1 — نهاية الاقتران الضعيف لـ HIH_I: يحتوي هاملتوني التفاعل الكامل HI(LSH)H_I^{\text{(LSH)}} (المعادلة 16 في [1]) على عوامل مسبقة تعتمد على العدد الكمومي البوزوني nln_l عبر حدود مثل 1/nl+11/\sqrt{n_l+1}. في نظام الاقتران الضعيف (x≫1x \gg 1)، تهيمن على الديناميكيات الحد الكهربائي HEH_E، الذي يفضّل الحالات ذات nln_l الكبير. من أجل nl≫1n_l \gg 1، فإن النسبة nl/(nl+1)→1n_l/(n_l+1) \to 1 وتُبسَّط جميع هذه العوامل المسبقة إلى الواحد. يختزل هاملتوني التفاعل عندئذ إلى قفز محلي بحت بين أقرب الجيران:

HIapprox=−∑r[σ−(r)σ+(r+1)+σ+(r)σ−(r+1)],H_I^{\text{approx}} = -\sum_r \left[\sigma^-(r)\sigma^+(r+1) + \sigma^+(r)\sigma^-(r+1)\right],

وهو مستقل عن nln_l ويؤثر فقط على الكيوبتات الفرميونية (ni,no)(n_i, n_o).

التقريب 2 — متوسط التدفق العام لـ HEH_E: تعتمد الطاقة الكهربائية على nln_l عند كل رابط. في فراغ الاقتران الضعيف، تكون nln_l كبيرة ومنتظمة تقريبًا. استبدل قيم nln_l المعتمدة على الموقع بمتوسط عام واحد nˉl\bar{n}_l، مما يجعل HEH_E طورًا قطريًا يتناسب مع تشكيل الفرميونات عند كل موقع:

HEapprox=NhE0+∑{r′}(nˉl2+34)H_E^{\text{approx}} = N h_E^0 + \sum_{\{r'\}} \left(\frac{\bar{n}_l}{2} + \frac{3}{4}\right)

حيث تُجمَع {r′}\{r'\} على المواقع في تشكيل الفرميونات (ni=0,no=1)(n_i=0, n_o=1)، وhE0h_E^0 طور عام يمكنك تجاهله.

التقريب 3 — تحليل Trotter: يُحلَّل مؤثر التطور الزمني لخطوة مدتها δτ\delta_\tau على النحو التالي:

e−iδτW≈e−im~HM e−iδτHEapprox e−icHIapproxe^{-i\delta_\tau W} \approx e^{-i\tilde{m} H_M} \, e^{-i\delta_\tau H_E^{\text{approx}}} \, e^{-ic H_I^{\text{approx}}}

حيث c=δτxc = \delta_\tau x، وm~=δτμ\tilde{m} = \delta_\tau \mu، وθ=−δτ(nˉl/2+3/4)\theta = -\delta_\tau(\bar{n}_l/2 + 3/4). يُدخل هذا التحليل من الرتبة الأولى لـ Trotter خطأً يتلاشى مع δτ→0\delta_\tau \to 0. نُثبّت δτ=0.0015\delta_\tau = 0.0015 طوال الوقت.

النتيجة التي نحصل عليها من هذه التقريبات الثلاثة هي أن الكيوبتين الفرميونيين لكل موقع (ni,no)(n_i, n_o) فقط هما الديناميكيان — إذ استُوعبت درجة حرية البوزون nln_l في معاملات فعّالة. ينتج عن هذا دارة مدمجة تحتوي على 2N2N كيوبت لعدد NN من مواقع الشبكة، حيث يكون عمق بوابات الكيوبتين ثابتًا في كل خطوة Trotter (13 لكل خطوة).

ما يحاكيه هذا الدرس التعليمي​

يحاكي هذا الدرس التعليمي انتشار الهادرونات: بدءًا من فراغ الاقتران الشديد (حالة جداء)، ضع ميزونًا في مركز الشبكة وطوّره زمنيًا. يعزل بروتوكول القياس التفاضلي — تشغيل الدارة مع الميزون المركزي وبدونه ثم الطرح — الإشارة الهادرونية المتماسكة من كل من ضوضاء الجهاز والتأثيرات الحدية. والنتيجة هي نمط مخروط ضوئي من تذبذبات كثافة الفرميونات المميزة لنمط تنفس ميزون محصور.

المتطلبات​

قبل بدء هذا الدرس التعليمي، ثبّت ما يلي:

  • Qiskit SDK الإصدار 2.0 أو أحدث، مع دعم التصور

  • Qiskit Runtime الإصدار 0.22 أو أحدث (pip install qiskit-ibm-runtime)

  • حزمة Pauli Propagation (pip install pauli-prop)

  • NumPy (pip install numpy)

  • Matplotlib (pip install matplotlib)

الإعداد​

ابدأ باستيراد المكتبات اللازمة وتعريف الدوال المساعدة التي تبني الدارات الكمومية للتطور الزمني في إطار LSH. توجد ثلاث دوال أساسية لبناء الدارات:

  1. pair_hamiltonian_circuit: تنفّذ المؤثر الأحادي ثنائي الكيوبت UIU_I لهاملتوني التفاعل التقريبي بين المواقع المتجاورة. تحليل البوابات هو: CNOT→H→Rz(−c)→CNOT→Rz(c)→CNOT→H→CNOT\text{CNOT} \to H \to R_z(-c) \to \text{CNOT} \to R_z(c) \to \text{CNOT} \to H \to \text{CNOT}.

  2. electric_hamiltonian_circuit: تنفّذ المؤثر الأحادي ثنائي الكيوبت UEU_E لطاقة الحقل الكهربائي التقريبية عند كل موقع. تحليل البوابات هو: X→Rz(θ/2)→CNOT→Rz(−θ/2)→CNOT→Rz(θ/2)→XX \to R_z(\theta/2) \to \text{CNOT} \to R_z(-\theta/2) \to \text{CNOT} \to R_z(\theta/2) \to X.

  3. construct_circuit: تجمّع الدارة الكاملة المحلَّلة بطريقة Trotter، من خلال ترتيب حدود التفاعل والحقل الكهربائي والكتلة في طبقات مع بوابات SWAP لإدارة اتصال الكيوبتات.

# Added by doQumentation — required packages for this notebook
!pip install -q matplotlib numpy pauli-prop qiskit qiskit-ibm-runtime
# Import libraries

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import TwoSlopeNorm
from qiskit.circuit import QuantumCircuit
from qiskit.quantum_info import SparsePauliOp
from typing import Optional

import warnings

warnings.filterwarnings("ignore")
def pair_hamiltonian_circuit(c: float) -> QuantumCircuit:
    """Two-qubit unitary for the approximate interaction Hamiltonian H_I.

    Implements exp(-i * c * H_I^approx) for one pair of neighboring sites,
    where c = delta_tau * x.
    """
    qc_temp = QuantumCircuit(2)
    qc_temp.cx(1, 0)
    qc_temp.h(1)
    qc_temp.rz(-c, 1)
    qc_temp.cx(0, 1)
    qc_temp.rz(c, 1)
    qc_temp.cx(0, 1)
    qc_temp.h(1)
    qc_temp.cx(1, 0)
    return qc_temp

def electric_hamiltonian_circuit(theta: float) -> QuantumCircuit:
    """Two-qubit unitary for the approximate electric field Hamiltonian H_E.

    Implements exp(-i * theta * H_E^approx) for one lattice site,
    where theta = -delta_tau * (n_bar_l / 2 + 3/4).
    """
    qc_temp = QuantumCircuit(2)
    qc_temp.x(0)
    qc_temp.rz(theta / 2, 0)
    qc_temp.cx(0, 1)
    qc_temp.rz(-theta / 2, 1)
    qc_temp.cx(0, 1)
    qc_temp.rz(theta / 2, 1)
    qc_temp.x(0)
    return qc_temp

def construct_circuit(
    num_lattice_point: int,
    num_trotter_steps: int,
    c: float,
    theta: float,
    m: float,
    theory: Optional[int] = 2,
    barriers: Optional[bool] = False,
    measurement: Optional[bool] = False,
    add_init_state: Optional[bool] = True,
    inverse_mid: Optional[bool] = False,
) -> QuantumCircuit:
    """Construct the full Trotterized time-evolution circuit.

    Builds a circuit implementing n Trotter steps of the approximate SU(2)
    LSH Hamiltonian evolution. The qubit layout uses a zigzag ordering:
    n_i(0), n_i(1), n_o(0), n_o(1), n_i(2), n_i(3), n_o(2), n_o(3), ...
    which minimizes the number of SWAP layers needed.

    Args:
        num_lattice_point: Number of lattice sites
        (num_qubits = 2 * num_lattice_point).
        num_trotter_steps: Number of Trotter steps.
        c: Interaction parameter (delta_tau * x).
        theta: Electric field phase parameter.
        m: Mass parameter (m_tilde = delta_tau * mu).
        theory: 1 for single chain, 2 for SU(2). Default 2.
        barriers: Insert barriers between Trotter layers for
        visualization.
        measurement: Append measurements at the end.
        add_init_state: Prepare the half-filled (strong-coupling vacuum)
        initial state.
        inverse_mid: Swap the central sites
        (for differential measurement protocol).
    """
    num_qubits = theory * num_lattice_point
    qc = QuantumCircuit(num_qubits)

    if num_trotter_steps <= 0:
        return qc

    # --- Initial state preparation ---
    if add_init_state:
        i = 1
        while i < num_lattice_point:
            for j in range(theory):
                qc.x(i + j * num_lattice_point)
            i = i + 2
        if inverse_mid:
            mid_lattice_qubits = [num_qubits // 2 - 1, num_qubits // 2]
            qc.x(mid_lattice_qubits)
    else:
        i = 1
        while i < num_qubits - 1:
            qc.swap(i, i + 1)
            i = i + 4

    # --- Trotter steps ---
    for step in range(num_trotter_steps):
        if barriers:
            qc.barrier()

        # First SWAP layer (skipped at step 0 — absorbed into initial state mapping)
        if step > 0:
            i = 1
            while i < num_qubits - 1:
                qc.swap(i, i + 1)
                i = i + 4

        # First layer of pair interactions
        j = 0
        while j < num_qubits - 2:
            circ = pair_hamiltonian_circuit(c)
            qc.compose(circ, [j, j + 1], inplace=True)
            j = j + 2
        if num_lattice_point % 2 == 0:
            circ = pair_hamiltonian_circuit(c)
            qc.compose(circ, [j, j + 1], inplace=True)

        # Second SWAP layer
        i = 1
        while i < num_qubits - 1:
            qc.swap(i, i + 1)
            i = i + theory

        # Second layer of pair interactions
        j = 2
        while j < num_qubits - 3:
            circ = pair_hamiltonian_circuit(c)
            qc.compose(circ, [j, j + 1], inplace=True)
            j = j + 2
        if num_lattice_point % 2 != 0:
            circ = pair_hamiltonian_circuit(c)
            qc.compose(circ, [j, j + 1], inplace=True)

        # Third SWAP layer
        i = 3
        while i < num_qubits - 1:
            qc.swap(i, i + 1)
            i = i + 2 * theory

        # Electric field term
        if theta != 0:
            e_circ = electric_hamiltonian_circuit(theta)
            for j in range(num_lattice_point):
                qc.compose(e_circ, [2 * j, 2 * j + 1], inplace=True)

        # Mass term: Rz(-m_tilde) for even sites, Rz(m_tilde) for odd sites
        for q in range(num_qubits):
            if q % 2 == 0:
                qc.rz(-1 * m, q)
            else:
                qc.rz(m, q)

    if measurement:
        qc.measure_all()

    return qc
def get_probabilities(expval: float):
    """Convert a Z-expectation value to site occupation probability.

    Since <Z> = p(0) - p(1), the occupation probability is p(1) = (1 - <Z>) / 2.
    """
    p1 = round((1 - expval) / 2, 3)
    return p1

def get_number(expval_data, num_lattice_point):
    """Convert raw Z-expectation values to staggered fermion number n_f at each site.

    n_f(r) = n_i(r) + n_o(r)           for even r
    n_f(r) = 2 - [n_i(r) + n_o(r)]     for odd r

    The two qubits per site encode (n_i, n_o), and occupation probabilities
    give us <n_i> and <n_o>.
    """
    N = []
    for expvals in expval_data:
        Pstep = [get_probabilities(expval) for expval in expvals]
        Nstep = []
        for k in range(num_lattice_point):
            val = Pstep[2 * k] + Pstep[2 * k + 1]
            a = 2 * (k % 2) + (1 - 2 * (k % 2)) * val
            Nstep.append(float(a))
        N.append(Nstep)
    return N

def calculate_difference(N, N_mid, num_lattice_point):
    """Differential measurement protocol: |n_f(meson) - n_f(vacuum)|.

    Subtracting the vacuum (SCV) evolution from the meson evolution
    isolates the coherent hadron signal from symmetric noise and boundary effects.
    """
    N_diff = []
    for i in range(len(N)):
        Nstep_diff = []
        for j in range(num_lattice_point):
            Nstep_diff.append(abs(N[i][j] - N_mid[i][j]))
        N_diff.append(Nstep_diff)
    return N_diff

مثال محاكي صغير النطاق​

أولًا، وضّح سير العمل على نطاق صغير باستخدام شبكة من ستة مواقع (12 كيوبت)، حتى تتمكن من التحقق من بناء الدارة وفهم المُشاهَدات الفيزيائية قبل التشغيل على الجهاز.

الخطوة 1: تخطيط المدخلات الكلاسيكية إلى مسألة كمومية​

عرّف المعاملات الفيزيائية المطابقة لنظام الاقتران الضعيف المدروس في الورقة البحثية (x=100x = 100، m/g=1m/g = 1). معاملات الدارة المشتقة هي:

  • c=δτ⋅x=0.15c = \delta_\tau \cdot x = 0.15 (معامل التفاعل)

  • θ=−δτ(nˉl/2+3/4)=0.01\theta = -\delta_\tau (\bar{n}_l/2 + 3/4) = 0.01 (طور الحقل الكهربائي)

  • m~=δτ⋅μ=0.03\tilde{m} = \delta_\tau \cdot \mu = 0.03 (معامل الكتلة)

لكل عدد من خطوات Trotter، ابنِ دارتين: واحدة تُهيّئ ميزونًا في المركز (inverse_mid=True) وأخرى تُحضّر فراغ الاقتران الشديد (inverse_mid=False). يطرح بروتوكول القياس التفاضلي تطور الفراغ لعزل الإشارة الهادرونية.

# Physical / circuit parameters
num_lattice_point = 6  # 6 lattice sites -> 12 qubits for SU(2)
num_qubits = 2 * num_lattice_point
c = 0.15  # delta_tau * x
theta = 0.01  # electric field phase
m = 0.03  # m_tilde = delta_tau * mu
trotter_steps = range(1, 11)  # 10 Trotter steps

print(f"Lattice sites: {num_lattice_point}, Qubits: {num_qubits}")
print(f"Parameters: c={c}, theta={theta}, m_tilde={m}")
Lattice sites: 6, Qubits: 12
Parameters: c=0.15, theta=0.01, m_tilde=0.03
# Build circuits: meson initial state and vacuum (SCV) initial state
circuits_mid = [
    construct_circuit(
        num_lattice_point,
        d,
        c,
        theta,
        m,
        barriers=False,
        measurement=False,
        add_init_state=True,
        inverse_mid=True,
    )
    for d in trotter_steps
]

circuits = [
    construct_circuit(
        num_lattice_point,
        d,
        c,
        theta,
        m,
        barriers=False,
        measurement=False,
        add_init_state=True,
        inverse_mid=False,
    )
    for d in trotter_steps
]

# Visualize a single Trotter step
print(
    f"Circuit for 1 Trotter step: {circuits[0].num_qubits} qubits, depth {circuits[0].depth()}"
)
circuits[0].draw("mpl", fold=-1)
Circuit for 1 Trotter step: 12 qubits, depth 26

Output of the previous code cell

الخطوة 2: تحسين المسألة لتنفيذها على الجهاز الكمومي​

عرّف المُشاهَدات: قياسات ZZ لكيوبت واحد على كل كيوبت. يمكنك من ⟨Z⟩\langle Z \rangle استخراج احتمالات الإشغال ثم عدد الفرميونات المتعرج nf(r)n_f(r) عند كل موقع شبكي rr.

# Z observable on each qubit
observables = [
    SparsePauliOp("I" * i + "Z" + "I" * (num_qubits - i - 1))
    for i in range(num_qubits)
]

print(f"Number of observables: {len(observables)}")
Number of observables: 12

الخطوة 3: التنفيذ باستخدام أدوات Qiskit الأساسية​

استخدم StatevectorEstimator للمحاكاة الدقيقة الخالية من الضوضاء على نطاق صغير.

from qiskit.primitives import StatevectorEstimator

estimator = StatevectorEstimator()

# Run meson circuits
pubs_mid = [(circuit, observables) for circuit in circuits_mid]
result_mid = estimator.run(pubs_mid).result()

# Run vacuum (SCV) circuits
pubs = [(circuit, observables) for circuit in circuits]
result = estimator.run(pubs).result()

# Extract expectation values
raw_expvals_mid = [
    result_mid[i].data.evs[::-1] for i in range(len(circuits_mid))
]
raw_expvals = [result[i].data.evs[::-1] for i in range(len(circuits))]

print(f"Computed expectation values for {len(raw_expvals)} Trotter steps")
Computed expectation values for 10 Trotter steps

الخطوة 4: المعالجة اللاحقة وإرجاع النتيجة بالصيغة الكلاسيكية المطلوبة​

حوّل قيم التوقع إلى عدد الفرميونات المتعرج nf(r,t)n_f(r, t) وطبّق بروتوكول القياس التفاضلي (ميزون −- فراغ) لإنتاج الخريطة الحرارية لانتشار الهادرون. يعيد هذا إنتاج بنية الشكل 3 من الورقة البحثية المرجعية: موقع الشبكة rr على المحور السيني، وخطوة Trotter (الزمن) tt على المحور الصادي، وnf(r,t)n_f(r,t) كمقياس لوني.

# Compute fermion numbers
N_mid_sim = get_number(raw_expvals_mid, num_lattice_point)
N_sim = get_number(raw_expvals, num_lattice_point)
N_diff_sim = calculate_difference(N_mid_sim, N_sim, num_lattice_point)

# --- Reproduce Figure 3 style: Staggered Fermionic Occupation Number Dynamics ---
fig, axes = plt.subplots(1, 3, figsize=(18, 5))

# Convert to numpy arrays for plotting
N_mid_arr = np.array(N_mid_sim)
N_arr = np.array(N_sim)
N_diff_arr = np.array(N_diff_sim)

# Color scheme
vmax = max(max(sublist) for sublist in N_arr)
vmin = -vmax

# Meson evolution
norm1 = TwoSlopeNorm(vmin=vmin, vcenter=0, vmax=vmax)
im0 = axes[0].imshow(
    N_mid_arr,
    aspect="auto",
    origin="lower",
    cmap="RdBu_r",
    norm=norm1,
    extent=[0, num_lattice_point, 0.5, len(trotter_steps) + 0.5],
)
axes[0].set_xlabel("Lattice site $r$", fontsize=12)
axes[0].set_ylabel("Trotter step $t$", fontsize=12)
axes[0].set_title("$n_f(r,t)$ — Meson initial state", fontsize=12)
plt.colorbar(im0, ax=axes[0], label="$n_f(r,t)$")

# Vacuum (SCV) evolution
im1 = axes[1].imshow(
    N_arr,
    aspect="auto",
    origin="lower",
    cmap="RdBu_r",
    norm=norm1,
    extent=[0, num_lattice_point, 0.5, len(trotter_steps) + 0.5],
)
axes[1].set_xlabel("Lattice site $r$", fontsize=12)
axes[1].set_ylabel("Trotter step $t$", fontsize=12)
axes[1].set_title("$n_f(r,t)$ — Vacuum (SCV)", fontsize=12)
plt.colorbar(im1, ax=axes[1], label="$n_f(r,t)$")

# Differential: meson - vacuum
norm2 = TwoSlopeNorm(vmin=vmin, vcenter=0, vmax=vmax)
im2 = axes[2].imshow(
    N_diff_arr,
    aspect="auto",
    origin="lower",
    cmap="RdBu_r",
    norm=norm2,
    extent=[0, num_lattice_point, 0.5, len(trotter_steps) + 0.5],
)
axes[2].set_xlabel("Lattice site $r$", fontsize=12)
axes[2].set_ylabel("Trotter step $t$", fontsize=12)
axes[2].set_title(
    "Staggered Fermionic Occupation Number Dynamics\n$|n_f^{\\mathrm{meson}} - n_f^{\\mathrm{vacuum}}|$",
    fontsize=12,
)
plt.colorbar(im2, ax=axes[2], label="$n_f(r,t)$")

plt.suptitle(
    f"Hadron propagation — {num_lattice_point}-site lattice (StatevectorEstimator)",
    fontsize=14,
    y=1.02,
)
plt.tight_layout()
plt.show()

Output of the previous code cell

مثال على جهاز كبير النطاق​

ننتقل الآن إلى شبكة مكونة من 30 موقعًا (60 كيوبت) على جهاز IBM Quantum. عند هذا النطاق، تتكون الدارة عند 10 خطوات Trotter من أكثر من 3400 بوابة ثنائية الكيوبت و14000 بوابة أحادية الكيوبت.

الخطوات 1-4 (مضغوطة في كتلة برمجية واحدة)​

الجوانب الرئيسية لسير العمل على الجهاز:

  • 10 خطوات Trotter لدارتي الميزون والفراغ (متداخلة لتقليل الانحراف)

  • التصريف باستخدام optimization_level=1 — مخطط الدارة متماثل بالفعل مع طوبولوجيا الجهاز (سلسلة خطية)، لذا لا حاجة إلى بوابات SWAP للتوجيه. يُستخدم الـ Transpiler فقط لاختيار سلسلة منخفضة الضوضاء من الكيوبتات الفعلية وتحليل البوابات إلى مجموعة البوابات الأصلية.

  • EstimatorV2 مع تخفيف أخطاء القراءة TREX وتدوير باولي (Pauli twirling)

  • جلسة Batch لإرسال جميع المهام معًا

# -------------------------Step 1: Define parameters & build circuits-------------------------

from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
from qiskit_ibm_runtime import EstimatorV2, Batch
from qiskit_ibm_runtime.options import (
    EstimatorOptions,
    ResilienceOptionsV2,
    TwirlingOptions,
    DynamicalDecouplingOptions,
)

service = QiskitRuntimeService()

num_lattice_point_hw = 30
num_qubits_hw = 2 * num_lattice_point_hw  # 60 qubits
c_hw = 0.15
theta_hw = 0.01
m_hw = 0.03
trotter_steps_hw = range(1, 11)  # 10 Trotter steps

# Build meson and vacuum circuits
circuits_mid_hw = [
    construct_circuit(
        num_lattice_point_hw,
        d,
        c_hw,
        theta_hw,
        m_hw,
        barriers=False,
        measurement=False,
        add_init_state=True,
        inverse_mid=True,
    )
    for d in trotter_steps_hw
]

circuits_hw = [
    construct_circuit(
        num_lattice_point_hw,
        d,
        c_hw,
        theta_hw,
        m_hw,
        barriers=False,
        measurement=False,
        add_init_state=True,
        inverse_mid=False,
    )
    for d in trotter_steps_hw
]

print(f"Built {len(circuits_hw)} circuit pairs for {num_qubits_hw} qubits")

# -------------------------Step 2: Transpile for hardware-------------------------
# The circuit topology is a linear chain, isomorphic to the device topology.
# We use optimization_level=1 since no routing SWAPs are needed — the transpiler
# only needs to select a low-noise qubit chain and decompose to native gates.

backend = service.backend("ibm_boston")

layout = [
    140,
    141,
    142,
    143,
    136,
    123,
    122,
    121,
    116,
    101,
    102,
    103,
    96,
    83,
    82,
    81,
    76,
    61,
    62,
    63,
    64,
    65,
    66,
    67,
    68,
    69,
    78,
    89,
    88,
    87,
    97,
    107,
    106,
    105,
    117,
    125,
    126,
    127,
    137,
    147,
    148,
    149,
    150,
    151,
    152,
    153,
    154,
    155,
    139,
    135,
    134,
    133,
    132,
    131,
    130,
    129,
    118,
    109,
    110,
    111,
]

pm = generate_preset_pass_manager(
    optimization_level=1, backend=backend, initial_layout=layout
)

isa_circuits_mid = pm.run(circuits_mid_hw)
isa_circuits = pm.run(circuits_hw)

print(f"Transpiled circuits. Example depth: {isa_circuits[0].depth()}")

# Define and layout-map observables
observables_hw = [
    SparsePauliOp("I" * i + "Z" + "I" * (num_qubits_hw - i - 1))
    for i in range(num_qubits_hw)
]

isa_observables_mid = [
    [obs.apply_layout(isa_circuits_mid[i].layout) for obs in observables_hw]
    for i in range(len(isa_circuits_mid))
]
isa_observables = [
    [obs.apply_layout(isa_circuits[i].layout) for obs in observables_hw]
    for i in range(len(isa_circuits))
]

# Build PUBs — interleave meson and vacuum for each Trotter step
isa_pubs_mid = [
    (circ, obs) for circ, obs in zip(isa_circuits_mid, isa_observables_mid)
]
isa_pubs = [(circ, obs) for circ, obs in zip(isa_circuits, isa_observables)]

pubs_to_execute = [
    [isa_pubs_mid[i], isa_pubs[i]] for i in range(len(isa_pubs))
]

# -------------------------Step 3: Execute on hardware-------------------------

twirling_options = TwirlingOptions(
    enable_gates=True,
    enable_measure=True,
    shots_per_randomization="auto",
    strategy="active-circuit",
)

resilience_options = ResilienceOptionsV2(
    measure_mitigation=True,  # TREX readout error mitigation
    zne_mitigation=False,  # ZNE turned off
)

dd_options = DynamicalDecouplingOptions(
    enable=False  # Circuit is sufficiently dense
)

options = EstimatorOptions(
    resilience=resilience_options,
    twirling=twirling_options,
    dynamical_decoupling=dd_options,
    default_shots=10_000,
)

ids = []
with Batch(backend=backend) as batch:
    for idx, pub in enumerate(pubs_to_execute):
        print(f"Submitting job for Trotter step {idx + 1}")
        estimator = EstimatorV2(mode=batch, options=options)
        estimator.skip_transpilation = True
        job = estimator.run(pub)
        ids.append(job.job_id())
    batch_id = batch.session_id

job_info = {"ids": ids, "batch_id": batch_id}
print(f"Submitted {len(ids)} jobs. Batch ID: {batch_id}")
print(ids)
# -------------------------Step 4: Post-process results-------------------------

jobs = [service.job(job_id) for job_id in ids]
results = [job.result() for job in jobs]

# Extract expectation values (index 0 = meson, index 1 = vacuum)
raw_expvals_mid_hw = [result[0].data.evs[::-1] for result in results]
raw_expvals_hw = [result[1].data.evs[::-1] for result in results]

# Compute fermion numbers and differential
N_mid_hw = get_number(raw_expvals_mid_hw, num_lattice_point_hw)
N_hw = get_number(raw_expvals_hw, num_lattice_point_hw)
N_diff_hw = calculate_difference(N_mid_hw, N_hw, num_lattice_point_hw)
N_diff_hw_arr = np.array(N_diff_hw)

fig, ax = plt.subplots(figsize=(10, 6))
vmax = np.abs(N_hw).max()
vmin = -vmax
norm = TwoSlopeNorm(vmin=vmin, vcenter=0, vmax=vmax)
im = ax.imshow(
    N_diff_hw_arr,
    aspect="auto",
    origin="lower",
    cmap="RdBu_r",
    norm=norm,
    extent=[0, 8, 0.5, len(trotter_steps_hw) + 0.5],
)
ax.set_xlabel("Lattice site $r$", fontsize=13)
ax.set_ylabel("Trotter step $t$", fontsize=13)
ax.set_title(
    "Staggered Fermionic Occupation Number Dynamics\nQuantum Simulation on IBM Hardware — 30-site lattice (60 qubits)",
    fontsize=13,
)
cbar = plt.colorbar(im, ax=ax)
cbar.set_label("$n_f(r,t)$", fontsize=12)
plt.tight_layout()
plt.show()

Output of the previous code cell

قياس الأداء الكلاسيكي عبر Pauli Propagation​

توفر طريقة انتشار باولي (Pauli Propagation Method, PPM) محاكاة كلاسيكية خالية من الضوضاء للدارة الكمومية من خلال النشر العكسي للمُشاهَدات المقيسة عبر الدارة في صورة هايزنبرغ. تحت طبقات كليفورد (بوابات CNOT وH وS وX)، تُخطَّط مؤثرات باولي إلى مؤثرات باولي أخرى دون زيادة عدد الحدود. يمكن للطبقات غير الكليفوردية (بوابات RzR_z في الدارة) أن تسبب تفرعًا — في أسوأ الحالات، تضاعف عدد الحدود — لكن العديد من الفروع لها معاملات صغيرة ويمكن اقتطاعها.

سير العمل باستخدام pauli-prop هو:

  1. قسّم الدارة إلى جزأيها الكليفوردي وغير الكليفوردي باستخدام evolve_through_cliffords.

  2. انشر كل مُشاهَد عبر الجزء غير الكليفوردي باستخدام propagate_through_circuit، مع الاحتفاظ بحد أقصى قدره max_terms من حدود باولي وإسقاط الحدود ذات المعاملات الأقل من عتبة الاقتطاع atol.

  3. طوّر النتيجة عبر الجزء الكليفوردي باستخدام دعم كليفورد المدمج في Qiskit.

  4. استخرج قيمة التوقع بجمع معاملات حدود باولي القطرية (التي تحتوي فقط على II وZZ).

عتبة الاقتطاع​

يتحكم المعامل atol في propagate_through_circuit في مدى قوة تقليم الفروع الصغيرة لباولي. تحتفظ عتبة ضيقة جدًا (على سبيل المثال، 1e-12) بجميع الفروع تقريبًا وتعطي نتائج دقيقة، لكن زمن المحاكاة يزداد بشدة مع عمق الدارة؛ استغرقت محاكاة 120 كيوبت في الورقة البحثية نحو 8.5 ساعة بالإعدادات الافتراضية. يؤدي رفع العتبة (على سبيل المثال، إلى 1e-6 أو 1e-3) إلى إسقاط الحدود التي تقل معاملاتها عن تلك القيمة، مما يقلل بشكل كبير من عدد الحدود المتتبَّعة ويسرّع الحساب. والمقايضة هي خطأ تقريب صغير يمكن التحكم فيه ويمكنك التحقق منه بمقارنة النتائج عند عتبات مختلفة.

import time
from pauli_prop import evolve_through_cliffords, propagate_through_circuit

# ── PPM Configuration ──
# Truncation threshold: controls the speed/accuracy trade-off.
PPM_THRESHOLD = 1e-3

# Maximum Pauli terms to track per observable (hard cap on memory/time)
PPM_MAX_TERMS = 66_000

print(f"PPM settings: atol={PPM_THRESHOLD}, max_terms={PPM_MAX_TERMS}")

# We propagate each single-qubit Z observable through each circuit.
# For PPM, we work with the un-transpiled circuits (ideal noiseless simulation).

observables_pp = [
    SparsePauliOp("I" * i + "Z" + "I" * (num_qubits_hw - i - 1))
    for i in range(num_qubits_hw)
]

def ppm_expectation_values(
    circuit, observables, max_terms=PPM_MAX_TERMS, atol=PPM_THRESHOLD
):
    """Compute expectation values of single-qubit Z observables
    via Pauli propagation.

    Args:
        circuit: The quantum circuit to simulate.
        observables: List of single-qubit Z observables.
        max_terms: Maximum number of Pauli terms to retain (hard cap).
        atol: Absolute tolerance — Pauli terms with coefficients below this
              value are discarded during propagation. Larger values give
              faster simulation at the cost of approximation accuracy.
    """
    circuit = circuit.decompose(["swap"])  # decompose SWAPs into 3 CX gates
    cliff, non_cliff = evolve_through_cliffords(circuit)

    evs = []
    for obs in observables:
        evolved_obs = propagate_through_circuit(
            obs, non_cliff, max_terms=max_terms, atol=atol, frame="h"
        )[0]
        evolved_obs.paulis = evolved_obs.paulis.evolve(cliff, frame="h")
        diagonal_mask = ~evolved_obs.paulis.x.any(axis=1)
        ev = float(evolved_obs.coeffs[diagonal_mask].sum().real)
        evs.append(ev)
    return np.array(evs)

# Run PPM for each Trotter step and record wall-clock time
pp_expvals_mid = []
pp_expvals = []
pp_times = []

for idx, d in enumerate(trotter_steps_hw):
    t_start = time.perf_counter()

    # Meson circuit
    evs_mid = ppm_expectation_values(circuits_mid_hw[idx], observables_pp)

    # Vacuum circuit
    evs_vac = ppm_expectation_values(circuits_hw[idx], observables_pp)

    elapsed = time.perf_counter() - t_start
    pp_times.append(elapsed)

    pp_expvals_mid.append(evs_mid[::-1])
    pp_expvals.append(evs_vac[::-1])

    print(f"Trotter step {d:2d}: {elapsed:.1f} s")

print(f"\nTotal PPM simulation time: {sum(pp_times):.1f} s")
print(f"Truncation threshold used: {PPM_THRESHOLD}")
PPM settings: atol=0.001, max_terms=66000
Trotter step  1: 5.0 s
Trotter step  2: 7.5 s
Trotter step  3: 11.2 s
Trotter step  4: 14.7 s
Trotter step  5: 18.3 s
Trotter step  6: 22.1 s
Trotter step  7: 25.6 s
Trotter step  8: 29.4 s
Trotter step  9: 33.2 s
Trotter step 10: 36.6 s

Total PPM simulation time: 203.6 s
Truncation threshold used: 0.001
# --- PPM simulation time vs. Trotter steps ---
fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(
    list(trotter_steps_hw),
    pp_times,
    "o-",
    color="tab:blue",
    linewidth=2,
    markersize=6,
)
ax.set_xlabel("Trotter step", fontsize=13)
ax.set_ylabel("Wall-clock time (s)", fontsize=13)
ax.set_title(
    "Pauli Propagation simulation time vs. Trotter steps\n(30-site lattice, 60 qubits)",
    fontsize=13,
)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

Output of the previous code cell

# --- PPM heatmap and comparison with hardware ---
N_mid_pp = get_number(pp_expvals_mid, num_lattice_point_hw)
N_pp = get_number(pp_expvals, num_lattice_point_hw)
N_diff_pp = calculate_difference(N_mid_pp, N_pp, num_lattice_point_hw)

N_diff_pp_arr = np.array(N_diff_pp)

fig, axes = plt.subplots(1, 2, figsize=(18, 6))
vmax = np.abs(N_hw).max()
vmin = -vmax
norm = TwoSlopeNorm(vmin=vmin, vcenter=0, vmax=vmax)

# PPM result
im0 = axes[0].imshow(
    N_diff_pp_arr,
    aspect="auto",
    origin="lower",
    cmap="RdBu_r",
    norm=norm,
    extent=[0, num_lattice_point_hw, 0.5, len(trotter_steps_hw) + 0.5],
)
axes[0].set_xlabel("Lattice site $r$", fontsize=12)
axes[0].set_ylabel("Trotter step $t$", fontsize=12)
axes[0].set_title(
    "Pauli Propagation\n(classical noiseless simulation)", fontsize=12
)
plt.colorbar(im0, ax=axes[0], label="$n_f(r,t)$")

# Hardware result
im1 = axes[1].imshow(
    N_diff_hw_arr,
    aspect="auto",
    origin="lower",
    cmap="RdBu_r",
    norm=norm,
    extent=[0, num_lattice_point_hw, 0.5, len(trotter_steps_hw) + 0.5],
)
axes[1].set_xlabel("Lattice site $r$", fontsize=12)
axes[1].set_ylabel("Trotter step $t$", fontsize=12)
axes[1].set_title(
    "Quantum Simulation\n(IBM Hardware, readout error mitigation only)",
    fontsize=12,
)
plt.colorbar(im1, ax=axes[1], label="$n_f(r,t)$")

plt.suptitle(
    "Staggered Fermionic Occupation Number Dynamics — 30-site lattice",
    fontsize=14,
    y=1.02,
)
plt.tight_layout()
plt.show()

Output of the previous code cell

الخطوات التالية​

إذا وجدت هذا العمل مثيرًا للاهتمام، ففكر في استكشاف المواد التالية:

التوصيات

المراجع​

[1] الورقة البحثية الأصلية: Ilčić، Majumdar، Mathew وآخرون. "Observation of Robust and Coherent Non-Abelian Hadron Dynamics on Noisy Quantum Processors" arXiv:2602.18080 (2026)

المصدر: IBM Quantum docs — تم التحديث في 9 سبتمبر 2026
النسخة الإنجليزية على doQumentation — تم التحديث في 9 سبتمبر 2026
مشكلة في الموقع أو الترجمة؟مشكلة في المحتوى؟ حرّر في IBM Quantum docs