Zum Hauptinhalt springen

Beobachtung robuster und kohärenter nicht-abelscher Hadronen-Dynamik auf verrauschten Quantenprozessoren

Nutzungsschätzung: 6 Minuten auf einem Heron-Prozessor (ibm_boston oder gleichwertig) (HINWEIS: Dies ist nur eine Schätzung. Deine Laufzeit kann variieren.)

Lernziele​

  • Wie sich nicht-abelsche Gittereichtheorien (insbesondere SU(2)) mithilfe des Loop-String-Hadron-(LSH)-Frameworks für eine effiziente Quantensimulation umformulieren lassen

  • Wie man trotterisierte Zeitentwicklungs-Circuits für einen approximierten SU(2)-Eichtheorie-Hamiltonian konstruiert und auf Qubits abbildet

  • Wie man diese Circuits auf IBM-Quantum®-Hardware mit der Qiskit-Estimator-Primitive und Readout-Fehlerminderung ausführt

Voraussetzungen​

Hintergrund​

Motivation​

Die Quantenchromodynamik (QCD), die SU(3)-Eichtheorie der starken Kraft, bindet Quarks zu Hadronen und bestimmt Confinement und String-Breaking. Klassische Gitter-QCD-Methoden sind bei statischen Eigenschaften hervorragend, können aber aufgrund des Vorzeichenproblems keine Echtzeitdynamik simulieren. Quantencomputer bieten einen Weg um diese Barriere herum, indem sie die Freiheitsgrade des Eichfelds direkt auf Qubits kodieren.

Dieses Tutorial demonstriert eine solche Simulation: Verwende IBM-Quantum-Hardware, um die Echtzeit-Ausbreitung von Hadronen in einer (1+1)-dimensionalen SU(2)-Gittereichtheorie zu simulieren — der einfachsten nicht-abelschen Eichtheorie und einem Sprungbrett hin zur vollständigen QCD.

Der Kogut-Susskind-Hamiltonian​

Die Theorie ist auf einem eindimensionalen Raumgitter formuliert, mit gestaffelten Fermionen (Materie) auf den Gitterplätzen und SU(2)-Eichfeldern auf den Verbindungen. Nach Umskalierung in dimensionslose Form lautet der Hamiltonian:

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

wobei HEH_E die chromoelektrische Feldenergie ist, HMH_M der gestaffelte Massenterm, HIH_I der Materie-Eich-Wechselwirkungsterm (Hopping-Term), μ=2mgx\mu = 2\frac{m}{g}\sqrt{x} die Fermionenmasse kodiert und x=1g2a2x = \frac{1}{g^2 a^2} die Wechselwirkungsstärke ist. Der Kontinuumslimes der Theorie liegt bei N→∞N \to \infty und x→∞x \to \infty.

Das Loop-String-Hadron-(LSH)-Framework​

Eine zentrale Herausforderung besteht darin, dass der Hilbertraum des Eichfelds an jeder Verbindung unendlichdimensional ist. Das Loop-String-Hadron-(LSH)-Framework begegnet dem, indem es die Theorie in gauge-invarianten Variablen umformuliert — Flussschleifen, Strings, die getrennte Ladungen verbinden, und Hadronen (eichinvariante Fermionpaare an einem Gitterplatz). In der LSH-Basis ist das Gauß'sche Gesetz konstruktionsbedingt automatisch erfüllt, sodass jeder Basiszustand physikalisch ist. Jeder Gitterplatz wird durch drei Quantenzahlen (nl,ni,no)(n_l, n_i, n_o) charakterisiert, die die Schleifenzahl, den eingehenden String und den ausgehenden String repräsentieren, wobei ni,no∈{0,1}n_i, n_o \in \{0,1\} fermionisch und nl≥0n_l \geq 0 bosonisch ist. Die lokale Fermionenzahl ist daraus definiert als nf(r)=ni(r)+no(r)n_f(r) = n_i(r) + n_o(r) für gerade Gitterplätze und nf(r)=2−[ni(r)+no(r)]n_f(r) = 2 - [n_i(r) + n_o(r)] für ungerade Gitterplätze.

Vom vollständigen Hamiltonian zum Quanten-Circuit: drei zentrale Näherungen​

Der Quanten-Circuit simuliert nicht den vollständigen SU(2)-Hamiltonian exakt. Stattdessen implementiert er eine kontrollierte Reihe von Näherungen, die im Regime schwacher Kopplung (x≫1x \gg 1) gültig sind. Es ist wichtig zu verstehen, was genähert wird und was nicht:

Näherung 1 — Grenzfall schwacher Kopplung für HIH_I: Der vollständige Wechselwirkungs-Hamiltonian HI(LSH)H_I^{\text{(LSH)}} (Gl. 16 in [1]) enthält Vorfaktoren, die über Terme wie 1/nl+11/\sqrt{n_l+1} von der bosonischen Quantenzahl nln_l abhängen. Im Regime schwacher Kopplung (x≫1x \gg 1) wird die Dynamik vom elektrischen Term HEH_E dominiert, der Zustände mit großem nln_l bevorzugt. Für nl≫1n_l \gg 1 geht das Verhältnis nl/(nl+1)→1n_l/(n_l+1) \to 1, und alle diese Vorfaktoren vereinfachen sich zu eins. Der Wechselwirkungs-Hamiltonian reduziert sich dann auf ein rein lokales Nächste-Nachbar-Hopping:

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

der unabhängig von nln_l ist und nur auf die fermionischen (ni,no)(n_i, n_o)-Qubits wirkt.

Näherung 2 — Globaler gemittelter Fluss für HEH_E: Die elektrische Energie hängt an jeder Verbindung von nln_l ab. Im Vakuum schwacher Kopplung ist nln_l groß und näherungsweise gleichförmig. Ersetze die ortsabhängigen nln_l-Werte durch einen einzigen globalen Mittelwert nˉl\bar{n}_l, wodurch HEH_E zu einer diagonalen Phase wird, die proportional zur Fermionenkonfiguration an jedem Gitterplatz ist:

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)

wobei {r′}\{r'\} über die Gitterplätze mit der fermionischen Konfiguration (ni=0,no=1)(n_i=0, n_o=1) summiert und hE0h_E^0 eine globale Phase ist, die du ignorieren kannst.

Näherung 3 — Trotterisierung: Der Zeitentwicklungsoperator für einen Schritt der Dauer δτ\delta_\tau wird zerlegt als:

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

wobei c=δτxc = \delta_\tau x, m~=δτμ\tilde{m} = \delta_\tau \mu und θ=−δτ(nˉl/2+3/4)\theta = -\delta_\tau(\bar{n}_l/2 + 3/4). Diese Trotter-Zerlegung erster Ordnung führt einen Fehler ein, der für δτ→0\delta_\tau \to 0 verschwindet. Wir setzen durchgängig δτ=0.0015\delta_\tau = 0.0015.

Das Ergebnis dieser drei Näherungen ist, dass nur die beiden fermionischen Qubits pro Gitterplatz (ni,no)(n_i, n_o) dynamisch sind — der bosonische Freiheitsgrad nln_l wurde in effektive Parameter absorbiert. Dies ergibt einen kompakten Circuit mit 2N2N Qubits für NN Gitterplätze, bei dem jeder Trotter-Schritt eine konstante Zwei-Qubit-Gate-Tiefe hat (13 pro Schritt).

Was dieses Tutorial simuliert​

Das Tutorial simuliert die Hadronen-Ausbreitung: Ausgehend vom Vakuum starker Kopplung (ein Produktzustand) wird ein Meson in der Mitte des Gitters platziert und zeitlich entwickelt. Das differenzielle Messprotokoll — Ausführen des Circuits mit und ohne zentrales Meson und anschließendes Subtrahieren — isoliert das kohärente Hadronensignal sowohl vom Hardware-Rauschen als auch von Randeffekten. Das Ergebnis ist ein Lichtkegelmuster von Fermionendichte-Oszillationen, das für die Atmungsmode eines eingeschlossenen (confined) Mesons charakteristisch ist.

Anforderungen​

Bevor du mit diesem Tutorial beginnst, installiere Folgendes:

  • Qiskit SDK v2.0 oder höher, mit Unterstützung für Visualisierung

  • Qiskit Runtime v0.22 oder höher (pip install qiskit-ibm-runtime)

  • Pauli-Propagation-Paket (pip install pauli-prop)

  • NumPy (pip install numpy)

  • Matplotlib (pip install matplotlib)

Einrichtung​

Beginne damit, die notwendigen Bibliotheken zu importieren und die Hilfsfunktionen zu definieren, die die Quanten-Circuits für die LSH-Zeitentwicklung erstellen. Es gibt drei zentrale Circuit-Erstellungsfunktionen:

  1. pair_hamiltonian_circuit: Implementiert die Zwei-Qubit-Unitäre UIU_I für den approximierten Wechselwirkungs-Hamiltonian zwischen benachbarten Gitterplätzen. Die Gate-Zerlegung lautet: 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: Implementiert die Zwei-Qubit-Unitäre UEU_E für die approximierte elektrische Feldenergie an jedem Gitterplatz. Die Gate-Zerlegung lautet: 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: Setzt den vollständigen trotterisierten Circuit zusammen, indem Wechselwirkungs-, elektrische und Massenterme mit SWAP-Gates geschichtet werden, um die Qubit-Konnektivität zu handhaben.

# 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

Kleinskaliges Simulator-Beispiel​

Demonstriere den Workflow zunächst im kleinen Maßstab anhand eines Sechs-Platz-Gitters (12 Qubits), damit du die Circuit-Konstruktion überprüfen und die physikalischen Observablen verstehen kannst, bevor du sie auf der Hardware ausführst.

Schritt 1: Klassische Eingaben auf ein Quantenproblem abbilden​

Definiere die physikalischen Parameter, die dem im Paper untersuchten Regime schwacher Kopplung entsprechen (x=100x = 100, m/g=1m/g = 1). Die abgeleiteten Circuit-Parameter sind:

  • c=δτ⋅x=0.15c = \delta_\tau \cdot x = 0.15 (Wechselwirkungsparameter)

  • θ=−δτ(nˉl/2+3/4)=0.01\theta = -\delta_\tau (\bar{n}_l/2 + 3/4) = 0.01 (Phase des elektrischen Felds)

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

Baue für jede Trotter-Schrittzahl zwei Circuits: einen, der ein Meson in der Mitte initialisiert (inverse_mid=True), und einen, der das Vakuum starker Kopplung vorbereitet (inverse_mid=False). Das differenzielle Messprotokoll subtrahiert die Vakuum-Entwicklung, um das Hadronensignal zu isolieren.

# 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

Schritt 2: Problem für die Ausführung auf Quantenhardware optimieren​

Definiere die Observablen: Einzelqubit-ZZ-Messungen an jedem Qubit. Aus ⟨Z⟩\langle Z \rangle kannst du Besetzungswahrscheinlichkeiten und anschließend die gestaffelte Fermionenzahl nf(r)n_f(r) an jedem Gitterplatz rr extrahieren.

# 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

Schritt 3: Ausführen mit Qiskit-Primitiven​

Verwende StatevectorEstimator für eine exakte, rauschfreie Simulation im kleinen Maßstab.

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

Schritt 4: Ergebnis nachbearbeiten und im gewünschten klassischen Format zurückgeben​

Wandle die Erwartungswerte in die gestaffelte Fermionenzahl nf(r,t)n_f(r, t) um und wende das differenzielle Messprotokoll (Meson −- Vakuum) an, um die Heatmap der Hadronausbreitung zu erzeugen. Dies reproduziert die Struktur von Abbildung 3 aus dem Referenzpapier: Gitterplatz rr auf der x-Achse, Trotter-Schritt (Zeit) tt auf der y-Achse und nf(r,t)n_f(r,t) als Farbskala.

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

Ausgabe der vorherigen Codezelle

Hardware-Beispiel im großen Maßstab​

Wir skalieren nun auf ein Gitter mit 30 Plätzen (60 Qubits) auf IBM-Quantum-Hardware hoch. In diesem Maßstab umfasst der Circuit bei 10 Trotter-Schritten über 3400 Zwei-Qubit-Gates und 14.000 Ein-Qubit-Gates.

Schritte 1–4 (zu einem einzigen Codeblock zusammengefasst)​

Wichtige Aspekte des Hardware-Workflows:

  • 10 Trotter-Schritte für die Meson- und Vakuum-Circuits (verschachtelt für minimale Drift)

  • Transpilation mit optimization_level=1 — das Circuit-Layout ist bereits isomorph zur Gerätetopologie (eine lineare Kette), sodass keine Routing-SWAPs erforderlich sind. Der Transpiler wird ausschließlich verwendet, um eine rauscharme Kette physischer Qubits auszuwählen und Gates in den nativen Gate-Satz zu zerlegen.

  • EstimatorV2 mit TREX-Readout-Fehlerminderung und Pauli-Twirling

  • Batch-Session, um alle Jobs gemeinsam einzureichen

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

Ausgabe der vorherigen Codezelle

Klassisches Benchmarking mittels Pauli-Propagation​

Die Pauli-Propagations-Methode (PPM) liefert eine rauschfreie klassische Simulation des Quanten-Circuits, indem gemessene Observablen im Heisenberg-Bild rückwärts durch den Circuit propagiert werden. Bei Clifford-Schichten (CNOT-, H-, S-, X-Gates) werden Pauli-Operatoren auf andere Pauli-Operatoren abgebildet, ohne dass die Termanzahl wächst. Nicht-Clifford-Schichten (die RzR_z-Gates im Circuit) können eine Verzweigung verursachen — im schlimmsten Fall eine Verdopplung der Anzahl der Terme —, aber viele Zweige haben kleine Koeffizienten und können abgeschnitten werden.

Der Workflow mit pauli-prop sieht wie folgt aus:

  1. Aufteilen des Circuits in seine Clifford- und Nicht-Clifford-Teile mithilfe von evolve_through_cliffords.

  2. Propagieren jeder Observable durch den Nicht-Clifford-Teil mithilfe von propagate_through_circuit, wobei bis zu max_terms Pauli-Terme beibehalten und Terme mit Koeffizienten unterhalb des Abschneide-Schwellenwerts atol verworfen werden.

  3. Weiterentwickeln des Ergebnisses durch den Clifford-Teil mithilfe der integrierten Clifford-Unterstützung von Qiskit.

  4. Extrahieren des Erwartungswerts durch Summieren der Koeffizienten diagonaler Pauli-Terme (die nur II und ZZ enthalten).

Abschneide-Schwellenwert​

Der Parameter atol in propagate_through_circuit steuert, wie aggressiv kleine Pauli-Zweige entfernt werden. Ein sehr enger Schwellenwert (zum Beispiel 1e-12) behält nahezu alle Zweige bei und liefert exakte Ergebnisse, aber die Simulationszeit wächst mit der Circuit-Tiefe stark an; die 120-Qubit-Simulation im Paper dauerte mit den Standardeinstellungen etwa 8,5 Stunden. Eine Erhöhung des Schwellenwerts (zum Beispiel auf 1e-6 oder 1e-3) verwirft Terme, deren Koeffizienten unter diesem Wert liegen, wodurch die Anzahl der verfolgten Terme drastisch reduziert und die Berechnung beschleunigt wird. Der Kompromiss ist ein kleiner, kontrollierbarer Approximationsfehler, den du überprüfen kannst, indem du Ergebnisse bei unterschiedlichen Schwellenwerten vergleichst.

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

Ausgabe der vorherigen Codezelle

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

Ausgabe der vorherigen Codezelle

Nächste Schritte​

Wenn dich diese Arbeit angesprochen hat, lohnt sich ein Blick auf die folgenden Materialien:

Empfehlungen

Referenzen​

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

Quelle: IBM Quantum docs — aktualisiert 9. Sept. 2026
Englische Version auf doQumentation — aktualisiert 9. Sept. 2026
Problem mit Website oder Übersetzung?Inhaltliches Problem? Auf IBM Quantum docs bearbeiten