Přeskočit na hlavní obsah

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

Odhadované využití: 6 minut na procesoru Heron (ibm_boston nebo ekvivalent) (POZNÁMKA: Jedná se pouze o odhad. Skutečná doba běhu se může lišit.)

Výsledky učení​

  • Jak lze neabelovské teorie kalibračního pole (konkrétně SU(2)) přeformulovat pomocí frameworku Loop-String-Hadron (LSH) pro efektivní kvantovou simulaci

  • Jak sestavit trotterizované obvody časového vývoje pro přibližný Hamiltonián teorie kalibračního pole SU(2) a namapovat je na qubity

  • Jak spustit tyto obvody na hardwaru IBM Quantum® pomocí primitivy Qiskit Estimator s mitigací chyb čtení

Předpoklady​

Pozadí​

Motivace​

Kvantová chromodynamika (QCD), teorie kalibračního pole SU(3) silné interakce, váže kvarky do hadronů a řídí uvěznění (confinement) a přetržení struny. Klasické metody mřížkové QCD vynikají ve statických vlastnostech, ale nedokážou simulovat dynamiku v reálném čase kvůli problému znaménka. Kvantové počítače nabízejí cestu obejít tuto bariéru zakódováním stupňů volnosti kalibračního pole přímo do qubitů.

Tento tutoriál demonstruje takovou simulaci: použij hardware IBM Quantum k simulaci šíření hadronu v reálném čase v (1+1)-dimenzionální teorii kalibračního pole mřížky SU(2) — nejjednodušší neabelovské teorii kalibračního pole a odrazovém můstku směrem k plné QCD.

Kogut-Susskindův Hamiltonián​

Teorie je formulována na 1D prostorové mřížce se staggerovanými fermiony (hmota) na uzlech a kalibračními poli SU(2) na spojnicích. Po přeškálování do bezrozměrného tvaru je Hamiltonián:

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

kde HEH_E je energie chromoelektrického pole, HMH_M je staggerovaný hmotový člen, HIH_I je člen interakce hmoty s kalibračním polem (hopping), μ=2mgx\mu = 2\frac{m}{g}\sqrt{x} kóduje hmotnost fermionu a x=1g2a2x = \frac{1}{g^2 a^2} je síla interakce. Kontinuum limita teorie leží při N→∞N \to \infty a x→∞x \to \infty.

Framework Loop-String-Hadron (LSH)​

Klíčovou výzvou je, že Hilbertův prostor kalibračního pole na každé spojnici je nekonečně-dimenzionální. Framework Loop-String-Hadron (LSH) tento problém řeší přeformulováním teorie v termínech kalibračně invariantních proměnných — smyček toku, strun spojujících oddělené náboje a hadronů (kalibračně singletních párů fermionů v uzlu). V bázi LSH je Gaussův zákon splněn automaticky konstrukcí, takže každý bázový stav je fyzikální. Každý mřížkový uzel je charakterizován třemi kvantovými čísly (nl,ni,no)(n_l, n_i, n_o) reprezentujícími počet smyček, vstupní strunu a výstupní strunu, kde ni,no∈{0,1}n_i, n_o \in \{0,1\} jsou fermionová a nl≥0n_l \geq 0 je bosonové. Lokální fermionové číslo je z nich definováno jako nf(r)=ni(r)+no(r)n_f(r) = n_i(r) + n_o(r) pro sudé uzly a nf(r)=2−[ni(r)+no(r)]n_f(r) = 2 - [n_i(r) + n_o(r)] pro liché uzly.

Od plného Hamiltoniánu ke kvantovému obvodu: tři klíčové aproximace​

Kvantový obvod nesimuluje plný Hamiltonián SU(2) přesně. Místo toho implementuje řízenou sérii aproximací, které jsou platné v režimu slabé vazby (x≫1x \gg 1). Je nezbytné rozumět tomu, co je a co není aproximováno:

Aproximace 1 — limita slabé vazby pro HIH_I: Plný interakční Hamiltonián HI(LSH)H_I^{\text{(LSH)}} (rov. 16 v [1]) obsahuje předfaktory, které závisí na bosonovém kvantovém čísle nln_l prostřednictvím členů jako 1/nl+11/\sqrt{n_l+1}. V režimu slabé vazby (x≫1x \gg 1) je dynamika dominována elektrickým členem HEH_E, který upřednostňuje stavy s velkým nln_l. Pro nl≫1n_l \gg 1 platí nl/(nl+1)→1n_l/(n_l+1) \to 1 a všechny tyto předfaktory se zjednoduší na jednotku. Interakční Hamiltonián se pak redukuje na čistě lokální hopping mezi nejbližšími sousedy:

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

který je nezávislý na nln_l a působí pouze na fermionové qubity (ni,no)(n_i, n_o).

Aproximace 2 — globálně zprůměrovaný tok pro HEH_E: Elektrická energie závisí na nln_l na každé spojnici. Ve vakuu slabé vazby je nln_l velké a přibližně jednotné. Nahraď hodnoty nln_l závislé na uzlu jediným globálním průměrem nˉl\bar{n}_l, čímž se HEH_E stane diagonální fází úměrnou fermionové konfiguraci na každém uzlu:

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)

kde {r′}\{r'\} sčítá přes uzly ve fermionové konfiguraci (ni=0,no=1)(n_i=0, n_o=1) a hE0h_E^0 je globální fáze, kterou lze ignorovat.

Aproximace 3 — Trotterizace: Operátor časového vývoje pro krok o délce δτ\delta_\tau se rozkládá 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}}}

kde c=δτxc = \delta_\tau x, m~=δτμ\tilde{m} = \delta_\tau \mu a θ=−δτ(nˉl/2+3/4)\theta = -\delta_\tau(\bar{n}_l/2 + 3/4). Tento Trotterův rozklad prvního řádu zavádí chybu, která mizí s δτ→0\delta_\tau \to 0. V celém textu udržujeme δτ=0.0015\delta_\tau = 0.0015.

Výsledkem těchto tří aproximací je, že dynamické jsou pouze dva fermionové qubity na uzel (ni,no)(n_i, n_o) — bosonový stupeň volnosti nln_l byl absorbován do efektivních parametrů. To dává kompaktní obvod s 2N2N qubity pro NN mřížkových uzlů, kde každý Trotterův krok má konstantní hloubku dvouqubitových hradel (13 na krok).

Co tento tutoriál simuluje​

Tutoriál simuluje šíření hadronu: počínaje vakuem silné vazby (produktový stav) umísti do středu mřížky meson a evolvuj v čase. Diferenciální měřicí protokol — spuštění obvodu s centrálním mesonem a bez něj, a následné odečtení — izoluje koherentní signál hadronu jak od hardwarového šumu, tak od okrajových efektů. Výsledkem je vzor světelného kužele oscilací hustoty fermionů, charakteristický pro uvězněný dýchací mód mesonu.

Požadavky​

Před zahájením tohoto tutoriálu nainstaluj následující:

  • Qiskit SDK v2.0 nebo novější, s podporou vizualizace

  • Qiskit Runtime v0.22 nebo novější (pip install qiskit-ibm-runtime)

  • Balíček Pauli Propagation (pip install pauli-prop)

  • NumPy (pip install numpy)

  • Matplotlib (pip install matplotlib)

Nastavení​

Začni importem potřebných knihoven a definicí pomocných funkcí, které sestavují kvantové obvody pro časový vývoj LSH. Existují tři základní funkce pro sestavení obvodu:

  1. pair_hamiltonian_circuit: Implementuje dvouqubitovou unitáru UIU_I pro přibližný interakční Hamiltonián mezi sousedními uzly. Rozklad hradel je: 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 dvouqubitovou unitáru UEU_E pro přibližnou energii elektrického pole na každém uzlu. Rozklad hradel je: 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: Sestaví plný trotterizovaný obvod, vrství interakční, elektrický a hmotový člen se SWAP hradly pro zvládnutí konektivity qubitů.

# 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

Příklad na malém měřítku v simulátoru​

Nejprve demonstruj postup v malém měřítku pomocí šestiuzlové mřížky (12 qubitů), aby ses mohl/a ověřit konstrukci obvodu a pochopit fyzikální pozorovatelné veličiny před spuštěním na hardwaru.

Krok 1: Mapování klasických vstupů na kvantový problém​

Definuj fyzikální parametry odpovídající režimu slabé vazby zkoumanému v článku (x=100x = 100, m/g=1m/g = 1). Odvozené parametry obvodu jsou:

  • c=δτ⋅x=0.15c = \delta_\tau \cdot x = 0.15 (interakční parametr)

  • θ=−δτ(nˉl/2+3/4)=0.01\theta = -\delta_\tau (\bar{n}_l/2 + 3/4) = 0.01 (fáze elektrického pole)

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

Pro každý počet Trotterových kroků sestav dva obvody: jeden inicializující meson ve středu (inverse_mid=True) a jeden připravující vakuum silné vazby (inverse_mid=False). Diferenciální měřicí protokol odečítá evoluci vakua, aby izoloval signál 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: Optimalizace problému pro spuštění na kvantovém hardwaru​

Definuj pozorovatelné veličiny: jednoqubitová měření ZZ na každém qubitu. Z ⟨Z⟩\langle Z \rangle lze získat pravděpodobnosti obsazení a poté staggerované fermionové číslo nf(r)n_f(r) na každém mřížkovém uzlu 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: Spuštění pomocí primitiv Qiskit​

Použij StatevectorEstimator pro přesnou bezšumovou simulaci v malém měřítku.

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: Následné zpracování a vrácení výsledku v požadovaném klasickém formátu​

Převeď střední hodnoty na staggerované fermionové číslo nf(r,t)n_f(r, t) a aplikuj diferenciální měřicí protokol (meson −- vakuum) k vytvoření teplotní mapy šíření hadronu. Tím se reprodukuje struktura obrázku 3 z referenčního článku: mřížkový uzel rr na ose x, Trotterův krok (čas) tt na ose y a nf(r,t)n_f(r,t) jako barevná škála.

# 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

Příklad na velkém měřítku na hardwaru​

Nyní se rozšíříme na 30-uzlovou mřížku (60 qubitů) na hardwaru IBM Quantum. V tomto měřítku obvod při 10 Trotterových krocích obsahuje přes 3400 dvouqubitových hradel a 14 000 jednoqubitových hradel.

Kroky 1-4 (zkompresované do jednoho bloku kódu)​

Klíčové aspekty hardwarového postupu:

  • 10 Trotterových kroků pro mesonový a vakuový obvod (prokládané pro minimální drift)

  • Transpilace s optimization_level=1 — rozvržení obvodu je již izomorfní s topologií zařízení (lineární řetězec), takže nejsou potřeba žádné SWAPy pro směrování. Transpiler se používá výhradně k výběru málo zašuměného řetězce fyzických qubitů a rozkladu hradel do nativní sady hradel.

  • EstimatorV2 s mitigací chyb čtení TREX a Pauliho twirlingem

  • Relace Batch pro odeslání všech úloh společně

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

Klasické benchmarkování pomocí Pauli Propagation​

Metoda Pauli Propagation (PPM) poskytuje bezšumovou klasickou simulaci kvantového obvodu zpětným šířením měřených pozorovatelných veličin obvodem v Heisenbergově obraze. Pod Cliffordovými vrstvami (hradla CNOT, H, S, X) se Pauliho operátory zobrazují na jiné Pauliho operátory bez zvýšení počtu členů. Non-Cliffordovy vrstvy (hradla RzR_z v obvodu) mohou způsobit větvení — v nejhorším případě zdvojnásobení počtu členů — ale mnoho větví má malé koeficienty a lze je oříznout.

Postup s pauli-prop je:

  1. Rozděl obvod na Cliffordovu a non-Cliffordovu část pomocí evolve_through_cliffords.

  2. Propaguj každou pozorovatelnou veličinu přes non-Cliffordovu část pomocí propagate_through_circuit, udržuj až max_terms Pauliho členů a zahazuj členy s koeficienty pod prahovou hodnotou oříznutí atol.

  3. Evolvuj výsledek přes Cliffordovu část pomocí vestavěné podpory Cliffordu v Qiskitu.

  4. Extrahuj střední hodnotu sečtením koeficientů diagonálních Pauliho členů (obsahujících pouze II a ZZ).

Prahová hodnota oříznutí​

Parametr atol v propagate_through_circuit řídí, jak agresivně jsou malé Pauliho větve prořezávány. Velmi přísná prahová hodnota (například 1e-12) zachovává téměř všechny větve a dává přesné výsledky, ale doba simulace prudce roste s hloubkou obvodu; simulace 120 qubitů v článku trvala přibližně 8,5 hodiny s výchozím nastavením. Zvýšení prahové hodnoty (například na 1e-6 nebo 1e-3) zahodí členy, jejichž koeficienty klesnou pod tuto hodnotu, což dramaticky snižuje počet sledovaných členů a urychluje výpočet. Kompromisem je malá, kontrolovatelná chyba aproximace, kterou lze ověřit porovnáním výsledků při různých prahových hodnotách.

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

Další kroky​

Pokud tě tato práce zaujala, zvaž prozkoumání následujících materiálů:

Doporučení

Reference​

[1] Původní článek: Ilčić, Majumdar, Mathew a kol. "Observation of Robust and Coherent Non-Abelian Hadron Dynamics on Noisy Quantum Processors" arXiv:2602.18080 (2026)

Zdroj: IBM Quantum docs — aktualizováno 9. 9. 2026
Anglická verze na doQumentation — aktualizováno 9. 9. 2026
Problém s webem nebo překladem?Problém v obsahu? Upravit v IBM Quantum docs