Przejdź do głównej treści

Obserwacja odpornej i spójnej dynamiki nieabelowych hadronów na zaszumionych procesorach kwantowych

Szacowany czas: 6 minut na procesorze Heron (ibm_boston lub równoważny) (UWAGA: to jedynie szacunek. Twój rzeczywisty czas może się różnić.)

Efekty kształcenia​

  • Jak nieabelowe teorie cechowania na sieci (konkretnie SU(2)) można przeformułować za pomocą frameworku Loop-String-Hadron (LSH) do efektywnej symulacji kwantowej

  • Jak skonstruować obwody trotteryzowanej ewolucji w czasie dla przybliżonego hamiltonianu teorii cechowania SU(2) i zmapować je na kubity

  • Jak uruchomić te obwody na sprzęcie IBM Quantum®, używając prymitywu Qiskit Estimator z łagodzeniem błędów odczytu

Wymagania wstępne​

Tło​

Motywacja​

Chromodynamika kwantowa (QCD), teoria cechowania SU(3) oddziaływania silnego, wiąże kwarki w hadrony i rządzi uwięzieniem oraz zrywaniem strun. Klasyczne metody QCD na sieci doskonale radzą sobie z właściwościami statycznymi, ale nie mogą symulować dynamiki w czasie rzeczywistym z powodu problemu znaku. Komputery kwantowe oferują sposób ominięcia tej bariery poprzez kodowanie stopni swobody pól cechowania bezpośrednio na kubitach.

Ten tutorial demonstruje taką symulację: użycie sprzętu IBM Quantum do symulacji propagacji hadronów w czasie rzeczywistym w (1+1)-wymiarowej teorii cechowania SU(2) na sieci — najprostszej nieabelowej teorii cechowania i kroku na drodze do pełnej QCD.

Hamiltonian Koguta-Susskinda​

Teoria jest sformułowana na jednowymiarowej sieci przestrzennej ze staggerowanymi fermionami (materią) na węzłach i polami cechowania SU(2) na łączach. Po przeskalowaniu do postaci bezwymiarowej hamiltonian ma postać:

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

gdzie HEH_E to energia chromoelektrycznego pola, HMH_M to staggerowany wyraz masowy, HIH_I to wyraz oddziaływania materia-cechowanie (przeskakiwania), μ=2mgx\mu = 2\frac{m}{g}\sqrt{x} koduje masę fermionu, a x=1g2a2x = \frac{1}{g^2 a^2} to siła oddziaływania. Granica continuum teorii leży przy N→∞N \to \infty i x→∞x \to \infty.

Framework Loop-String-Hadron (LSH)​

Kluczowym wyzwaniem jest to, że przestrzeń Hilberta pola cechowania na każdym łączu jest nieskończenie wymiarowa. Framework Loop-String-Hadron (LSH) rozwiązuje ten problem, przeformułowując teorię w terminach zmiennych niezmienniczych względem cechowania — pętli strumienia, strun łączących rozdzielone ładunki oraz hadronów (par fermionów będących singletami cechowania w węźle). W bazie LSH prawo Gaussa jest spełnione automatycznie z konstrukcji, więc każdy stan bazowy jest fizyczny. Każdy węzeł sieci jest charakteryzowany przez trzy liczby kwantowe (nl,ni,no)(n_l, n_i, n_o) reprezentujące liczbę pętli, strumień wchodzący i strumień wychodzący, gdzie ni,no∈{0,1}n_i, n_o \in \{0,1\} są fermionowe, a nl≥0n_l \geq 0 jest bozonowe. Lokalna liczba fermionów jest zdefiniowana na ich podstawie jako nf(r)=ni(r)+no(r)n_f(r) = n_i(r) + n_o(r) dla węzłów parzystych i nf(r)=2−[ni(r)+no(r)]n_f(r) = 2 - [n_i(r) + n_o(r)] dla węzłów nieparzystych.

Od pełnego hamiltonianu do obwodu kwantowego: trzy kluczowe przybliżenia​

Obwód kwantowy nie symuluje pełnego hamiltonianu SU(2) dokładnie. Zamiast tego implementuje kontrolowaną serię przybliżeń, które są ważne w reżimie słabego sprzężenia (x≫1x \gg 1). Zrozumienie tego, co jest, a co nie jest przybliżone, jest kluczowe:

Przybliżenie 1 — Granica słabego sprzężenia dla HIH_I: Pełny hamiltonian oddziaływania HI(LSH)H_I^{\text{(LSH)}} (rów. 16 w [1]) zawiera prefaktory zależne od bozonowej liczby kwantowej nln_l poprzez wyrazy takie jak 1/nl+11/\sqrt{n_l+1}. W reżimie słabego sprzężenia (x≫1x \gg 1) dynamika jest zdominowana przez wyraz elektryczny HEH_E, który faworyzuje stany o dużym nln_l. Dla nl≫1n_l \gg 1 stosunek nl/(nl+1)→1n_l/(n_l+1) \to 1, a wszystkie te prefaktory upraszczają się do jedności. Hamiltonian oddziaływania redukuje się wówczas do czysto lokalnego przeskakiwania między najbliższymi sąsiadami:

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

które jest niezależne od nln_l i działa jedynie na fermionowe kubity (ni,no)(n_i, n_o).

Przybliżenie 2 — Globalny uśredniony strumień dla HEH_E: Energia elektryczna zależy od nln_l na każdym łączu. W próżni słabego sprzężenia nln_l jest duże i w przybliżeniu jednorodne. Zastępujemy zależne od węzła wartości nln_l pojedynczą globalną średnią nˉl\bar{n}_l, czyniąc HEH_E diagonalną fazą proporcjonalną do konfiguracji fermionów w każdym węźle:

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)

gdzie {r′}\{r'\} sumuje po węzłach w konfiguracji fermionowej (ni=0,no=1)(n_i=0, n_o=1), a hE0h_E^0 jest globalną fazą, którą można pominąć.

Przybliżenie 3 — Trotteryzacja: Operator ewolucji w czasie dla kroku o czasie trwania δτ\delta_\tau jest rozkładany jako:

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

gdzie c=δτxc = \delta_\tau x, m~=δτμ\tilde{m} = \delta_\tau \mu, oraz θ=−δτ(nˉl/2+3/4)\theta = -\delta_\tau(\bar{n}_l/2 + 3/4). Ta dekompozycja Trottera pierwszego rzędu wprowadza błąd, który zanika wraz z δτ→0\delta_\tau \to 0. W całym tekście przyjmujemy δτ=0.0015\delta_\tau = 0.0015.

Wynikiem tych trzech przybliżeń jest to, że dynamiczne pozostają tylko dwa kubity fermionowe na węzeł (ni,no)(n_i, n_o) — bozonowy stopień swobody nln_l został wchłonięty do parametrów efektywnych. Daje to zwarty obwód z 2N2N kubitami dla NN węzłów sieci, gdzie każdy krok Trottera ma stałą głębokość bramek dwukubitowych (13 na krok).

Co symuluje ten samouczek​

Samouczek symuluje propagację hadronu: zaczynając od próżni silnego sprzężenia (stanu produktowego), umieszcza mezon w centrum sieci i ewoluuje w czasie. Protokół pomiaru różnicowego — uruchomienie obwodu z centralnym mezonem i bez niego, a następnie odjęcie wyników — izoluje spójny sygnał hadronu zarówno od szumu sprzętowego, jak i efektów brzegowych. Wynikiem jest wzór stożka świetlnego oscylacji gęstości fermionów, charakterystyczny dla ograniczonego trybu oddychania mezonu.

Wymagania​

Przed rozpoczęciem tego samouczka zainstaluj następujące elementy:

  • Qiskit SDK w wersji 2.0 lub nowszej, z obsługą wizualizacji

  • Qiskit Runtime w wersji 0.22 lub nowszej (pip install qiskit-ibm-runtime)

  • Pakiet Pauli Propagation (pip install pauli-prop)

  • NumPy (pip install numpy)

  • Matplotlib (pip install matplotlib)

Konfiguracja​

Zacznij od zaimportowania niezbędnych bibliotek i zdefiniowania funkcji pomocniczych, które budują obwody kwantowe do ewolucji czasowej LSH. Istnieją trzy podstawowe funkcje budujące obwody:

  1. pair_hamiltonian_circuit: Implementuje dwukubitową unitarną UIU_I dla przybliżonego hamiltonianu oddziaływania między sąsiednimi węzłami. Dekompozycja bramek to: 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: Implementuje dwukubitową unitarną UEU_E dla przybliżonej energii pola elektrycznego w każdym węźle. Dekompozycja bramek to: 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: Składa pełny obwód Trotterowski, warstwując człony oddziaływania, pola elektrycznego i masy z bramkami SWAP w celu zarządzania łącznością kubitów.

# 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

Przykład na symulatorze w małej skali​

Najpierw zademonstruj przepływ pracy w małej skali, używając sieci sześciowęzłowej (12 kubitów), aby móc zweryfikować konstrukcję obwodu i zrozumieć fizyczne obserwable przed uruchomieniem na sprzęcie.

Krok 1: Odwzorowanie danych klasycznych na problem kwantowy​

Zdefiniuj parametry fizyczne odpowiadające reżimowi słabego sprzężenia badanemu w artykule (x=100x = 100, m/g=1m/g = 1). Wyprowadzone parametry obwodu to:

  • c=δτ⋅x=0.15c = \delta_\tau \cdot x = 0.15 (parametr oddziaływania)

  • θ=−δτ(nˉl/2+3/4)=0.01\theta = -\delta_\tau (\bar{n}_l/2 + 3/4) = 0.01 (faza pola elektrycznego)

  • m~=δτ⋅μ=0.03\tilde{m} = \delta_\tau \cdot \mu = 0.03 (parametr masy)

Dla każdej liczby kroków Trottera zbuduj dwa obwody: jeden inicjalizujący mezon w centrum (inverse_mid=True) i jeden przygotowujący próżnię silnego sprzężenia (inverse_mid=False). Protokół pomiaru różnicowego odejmuje ewolucję próżni, aby wyizolować sygnał hadronu.

# 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

Krok 2: Optymalizacja problemu pod kątem wykonania na sprzęcie kwantowym​

Zdefiniuj obserwable: pomiary jednokubitowe ZZ na każdym kubicie. Z ⟨Z⟩\langle Z \rangle możesz wyodrębnić prawdopodobieństwa obsadzenia, a następnie naprzemienną liczbę fermionów nf(r)n_f(r) w każdym węźle sieci 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

Krok 3: Wykonanie z użyciem prymitywów Qiskit​

Użyj StatevectorEstimator do dokładnej, bezszumowej symulacji w małej skali.

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

Krok 4: Postprocessing i zwrócenie wyniku w pożądanym formacie klasycznym​

Przekonwertuj wartości oczekiwane na naprzemienną liczbę fermionów nf(r,t)n_f(r, t) i zastosuj protokół pomiaru różnicowego (mezon −- próżnia), aby wygenerować mapę cieplną propagacji hadronu. Odtwarza to strukturę Rysunku 3 z artykułu referencyjnego: węzeł sieci rr na osi x, krok Trottera (czas) tt na osi y, oraz nf(r,t)n_f(r,t) jako skala kolorów.

# 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

Przykład na sprzęcie w dużej skali​

Skalujemy teraz do sieci 30-węzłowej (60 kubitów) na sprzęcie IBM Quantum. W tej skali obwód przy 10 krokach Trottera obejmuje ponad 3400 bramek dwukubitowych i 14 000 bramek jednokubitowych.

Kroki 1-4 (skompresowane do jednego bloku kodu)​

Kluczowe aspekty przepływu pracy na sprzęcie:

  • 10 kroków Trottera dla obwodów mezonu i próżni (przeplatane dla minimalnego dryftu)

  • Transpilacja z optimization_level=1 — układ obwodu jest już izomorficzny z topologią urządzenia (łańcuch liniowy), więc żadne bramki SWAP routingu nie są potrzebne. Transpiler jest używany wyłącznie do wyboru łańcucha kubitów fizycznych o niskim szumie i dekompozycji bramek na natywny zestaw bramek.

  • EstimatorV2 z korekcją błędów odczytu TREX i twirlingiem Pauliego

  • Sesja Batch do wspólnego przesłania wszystkich zadań

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

Klasyczne testy porównawcze za pomocą Pauli Propagation​

Metoda Pauli Propagation (PPM) zapewnia bezszumową klasyczną symulację obwodu kwantowego poprzez propagację wsteczną zmierzonych obserwabli przez obwód w obrazie Heisenberga. W warstwach Clifforda (bramki CNOT, H, S, X) operatory Pauliego są odwzorowywane na inne operatory Pauliego bez zwiększania liczby członów. Warstwy nie-Clifforda (bramki RzR_z w obwodzie) mogą powodować rozgałęzianie — w najgorszym przypadku podwajając liczbę członów — ale wiele gałęzi ma małe współczynniki i można je obciąć.

Przepływ pracy z pauli-prop wygląda następująco:

  1. Podziel obwód na jego części Clifforda i nie-Clifforda za pomocą evolve_through_cliffords.

  2. Propaguj każdą obserwablę przez część nie-Clifforda za pomocą propagate_through_circuit, zachowując do max_terms członów Pauliego i odrzucając człony o współczynnikach poniżej progu obcięcia atol.

  3. Ewoluuj wynik przez część Clifforda, korzystając z wbudowanej obsługi Clifforda w Qiskit.

  4. Wyodrębnij wartość oczekiwaną, sumując współczynniki diagonalnych członów Pauliego (zawierających tylko II i ZZ).

Próg obcięcia​

Parametr atol w propagate_through_circuit kontroluje, jak agresywnie przycinane są małe gałęzie Pauliego. Bardzo ścisły próg (na przykład 1e-12) zachowuje niemal wszystkie gałęzie i daje dokładne wyniki, ale czas symulacji rośnie stromo wraz z głębokością obwodu; symulacja 120-kubitowa w artykule zajęła około 8,5 godziny przy ustawieniach domyślnych. Podniesienie progu (na przykład do 1e-6 lub 1e-3) odrzuca człony, których współczynniki spadają poniżej tej wartości, drastycznie zmniejszając liczbę śledzonych członów i przyspieszając obliczenia. Kompromisem jest niewielki, kontrolowany błąd przybliżenia, który można zweryfikować, porównując wyniki przy różnych progach.

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

Kolejne kroki​

Jeśli ta praca wydała Ci się interesująca, rozważ zapoznanie się z następującymi materiałami:

Rekomendacje

Odniesienia​

[1] Oryginalny artykuł: Ilčić, Majumdar, Mathew et al. "Observation of Robust and Coherent Non-Abelian Hadron Dynamics on Noisy Quantum Processors" arXiv:2602.18080 (2026)

Źródło: IBM Quantum docs — zaktualizowano 9 wrz 2026
Wersja angielska na doQumentation — zaktualizowano 9 wrz 2026
Problem ze stroną lub tłumaczeniem?Problem z treścią? Edytuj w IBM Quantum docs