Перейти до основного вмісту

Observation of robust and coherent non-Abelian hadron dynamics on noisy quantum processors

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

Результати навчання​

  • Як неабелеві теорії ґраткової калібрувальної симетрії (зокрема SU(2)) можна переформулювати за допомогою фреймворку Loop-String-Hadron (LSH) для ефективного квантового моделювання

  • Як побудувати тротеризовані схеми часової еволюції для наближеного гамільтоніана SU(2) калібрувальної теорії та відобразити їх на кубіти

  • Як запускати ці схеми на апаратному забезпеченні IBM Quantum® за допомогою примітиву Qiskit Estimator зі зменшенням похибок зчитування

Передумови​

Передумови​

Мотивація​

Квантова хромодинаміка (КХД), SU(3) калібрувальна теорія сильної взаємодії, зв'язує кварки в адрони та керує утриманням і розривом струни. Класичні методи ґраткової КХД чудово підходять для статичних властивостей, але не можуть моделювати динаміку в реальному часі через проблему знаку. Квантові комп'ютери пропонують шлях подолання цієї перешкоди, кодуючи ступені свободи калібрувального поля безпосередньо на кубітах.

Цей туторіал демонструє таке моделювання: використати апаратне забезпечення IBM Quantum для моделювання поширення адронів у реальному часі в (1+1)-вимірній SU(2) теорії ґраткової калібрувальної симетрії — найпростішій неабелевій калібрувальній теорії та кроці до повної КХД.

Гамільтоніан Когута-Сасскінда​

Теорія формулюється на 1D просторовій ґратці зі ступінчастими ферміонами (матерією) на вузлах і 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.

Фреймворк Loop-String-Hadron (LSH)​

Ключовий виклик полягає в тому, що гільбертів простір калібрувального поля на кожному зв'язку є нескінченновимірним. Фреймворк Loop-String-Hadron (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 — Тротеризація: Оператор часової еволюції для кроку тривалістю δτ\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). Це розкладання Троттера першого порядку вносить похибку, яка зникає при δτ→0\delta_\tau \to 0. Ми фіксуємо δτ=0.0015\delta_\tau = 0.0015 протягом усього туторіалу.

Результатом цих трьох наближень є те, що динамічними є лише два ферміонні кубіти на вузол (ni,no)(n_i, n_o) — бозонний ступінь свободи nln_l поглинутий ефективними параметрами. Це дає компактну схему з 2N2N кубітами для NN вузлів ґратки, де кожен крок Троттера має сталу глибину двокубітних вентилів (13 на крок).

Що моделює цей туторіал​

Туторіал моделює поширення адронів: починаючи з вакууму сильного зв'язку (добуткового стану), розміщує мезон у центрі ґратки та еволюціонує у часі. Диференціальний протокол вимірювання — запуск схеми з центральним мезоном і без нього, а потім віднімання — виділяє когерентний сигнал адрона як від апаратного шуму, так і від граничних ефектів. Результатом є світлоконусний патерн осциляцій густини ферміонів, характерний для утримуваного дихального моду мезона.

Вимоги​

Перед початком цього туторіалу встанови наступне:

  • Qiskit SDK v2.0 або новіше, з підтримкою візуалізації

  • Qiskit Runtime v0.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: Збирає повну тротеризовану схему, накладаючи члени взаємодії, електричного поля та маси разом з вентилями 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 (параметр маси)

Для кожної кількості кроків Троттера побудуй дві схеми: одну, що ініціалізує мезон у центрі (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 на осі x, крок Троттера (час) tt на осі y, і 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 кроках Троттера містить понад 3400 двокубітних вентилів і 14 000 однокубітних вентилів.

Кроки 1-4 (стиснуті в один блок коду)​

Ключові аспекти апаратного робочого процесу:

  • 10 кроків Троттера для схем мезона та вакууму (чергуються для мінімального дрейфу)

  • Транспіляція з optimization_level=1 — розташування схеми вже ізоморфне до топології пристрою (лінійний ланцюг), тому вентилі SWAP для маршрутизації не потрібні. Транспілятор використовується лише для вибору низькошумного ланцюга фізичних кубітів і розкладання вентилів у нативний набір вентилів.

  • EstimatorV2 зі зменшенням похибок зчитування TREX та скручуванням Паулі

  • Сесія 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​

Метод поширення Паулі (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 et al. "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