Lewati ke konten utama

Observasi dinamika hadron non-Abelian yang robust dan koheren pada prosesor kuantum yang noisy

Estimasi penggunaan: 6 menit pada prosesor Heron (ibm_boston atau setara) (CATATAN: Ini hanyalah estimasi. Runtime kamu mungkin bervariasi.)

Hasil pembelajaran​

  • Bagaimana teori gauge kisi non-Abelian (khususnya SU(2)) dapat diformulasikan ulang menggunakan framework Loop-String-Hadron (LSH) untuk simulasi kuantum yang efisien

  • Bagaimana membangun sirkuit evolusi waktu yang di-Trotterisasi untuk Hamiltonian teori gauge SU(2) yang diaproksimasi dan memetakannya ke qubit

  • Bagaimana menjalankan sirkuit-sirkuit ini pada hardware IBM Quantum® menggunakan primitive Qiskit Estimator dengan mitigasi error pembacaan

Prasyarat​

Latar Belakang​

Motivasi​

Quantum Chromodynamics (QCD), teori gauge SU(3) dari gaya kuat, mengikat quark menjadi hadron dan mengatur confinement dan string breaking. Metode lattice QCD klasik unggul dalam properti statis tapi tidak bisa mensimulasikan dinamika waktu-nyata karena sign problem. Komputer kuantum menawarkan jalan untuk melewati hambatan ini dengan mengkodekan derajat kebebasan medan gauge langsung ke qubit.

Tutorial ini mendemonstrasikan simulasi semacam itu: menggunakan hardware IBM Quantum untuk mensimulasikan propagasi hadron waktu-nyata dalam teori gauge lattice SU(2) (1+1)-dimensi — teori gauge non-Abelian paling sederhana dan batu loncatan menuju QCD penuh.

Hamiltonian Kogut-Susskind​

Teori ini diformulasikan pada lattice spasial 1D dengan fermion staggered (materi) pada situs dan medan gauge SU(2) pada link. Setelah penskalaan ulang ke bentuk tak berdimensi, Hamiltonian-nya adalah:

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

di mana HEH_E adalah energi medan chromoelectric, HMH_M adalah suku massa staggered, HIH_I adalah suku interaksi materi-gauge (hopping), μ=2mgx\mu = 2\frac{m}{g}\sqrt{x} mengkodekan massa fermion, dan x=1g2a2x = \frac{1}{g^2 a^2} adalah kekuatan interaksi. Limit kontinum dari teori ini berada pada N→∞N \to \infty dan x→∞x \to \infty.

Framework Loop-String-Hadron (LSH)​

Tantangan utamanya adalah ruang Hilbert medan gauge pada setiap link berdimensi tak hingga. Framework Loop-String-Hadron (LSH) mengatasi ini dengan memformulasikan ulang teori dalam variabel-variabel invarian gauge — loop fluks, string yang menghubungkan muatan-muatan terpisah, dan hadron (pasangan fermion gauge-singlet pada suatu situs). Dalam basis LSH, hukum Gauss terpenuhi secara otomatis berdasarkan konstruksinya, jadi setiap state basis bersifat fisis. Setiap situs lattice dikarakterisasi oleh tiga bilangan kuantum (nl,ni,no)(n_l, n_i, n_o) yang mewakili bilangan loop, string masuk, dan string keluar, di mana ni,no∈{0,1}n_i, n_o \in \{0,1\} bersifat fermionik dan nl≥0n_l \geq 0 bersifat bosonik. Bilangan fermion lokal didefinisikan dari nilai-nilai ini sebagai nf(r)=ni(r)+no(r)n_f(r) = n_i(r) + n_o(r) untuk situs genap dan nf(r)=2−[ni(r)+no(r)]n_f(r) = 2 - [n_i(r) + n_o(r)] untuk situs ganjil.

Dari Hamiltonian penuh ke sirkuit kuantum: tiga aproksimasi utama​

Sirkuit kuantum tidak mensimulasikan Hamiltonian SU(2) penuh secara eksak. Sebaliknya, sirkuit ini mengimplementasikan serangkaian aproksimasi terkontrol yang valid dalam rezim kopling lemah (x≫1x \gg 1). Memahami apa yang diaproksimasikan dan apa yang tidak sangatlah penting:

Aproksimasi 1 — Limit kopling lemah untuk HIH_I: Hamiltonian interaksi penuh HI(LSH)H_I^{\text{(LSH)}} (Persamaan 16 dalam [1]) berisi prefaktor yang bergantung pada bilangan kuantum bosonik nln_l melalui suku-suku seperti 1/nl+11/\sqrt{n_l+1}. Dalam rezim kopling lemah (x≫1x \gg 1), dinamika didominasi oleh suku elektrik HEH_E, yang menguntungkan state dengan nln_l besar. Untuk nl≫1n_l \gg 1, rasio nl/(nl+1)→1n_l/(n_l+1) \to 1 dan semua prefaktor ini tersederhanakan menjadi satu. Hamiltonian interaksi kemudian tereduksi menjadi hopping tetangga terdekat murni lokal:

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

yang tidak bergantung pada nln_l dan hanya bekerja pada qubit fermionik (ni,no)(n_i, n_o).

Aproksimasi 2 — Fluks rata-rata global untuk HEH_E: Energi elektrik bergantung pada nln_l di setiap link. Dalam vakum kopling lemah, nln_l besar dan kira-kira seragam. Gantikan nilai nln_l yang bergantung situs dengan satu rata-rata global nˉl\bar{n}_l, sehingga HEH_E menjadi fase diagonal yang proporsional terhadap konfigurasi fermion di setiap situs:

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)

di mana {r′}\{r'\} menjumlahkan situs-situs dalam konfigurasi fermionik (ni=0,no=1)(n_i=0, n_o=1), dan hE0h_E^0 adalah fase global yang bisa kamu abaikan.

Aproksimasi 3 — Trotterisasi: Operator evolusi-waktu untuk satu langkah berdurasi δτ\delta_\tau diuraikan sebagai:

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

di mana c=δτxc = \delta_\tau x, m~=δτμ\tilde{m} = \delta_\tau \mu, dan θ=−δτ(nˉl/2+3/4)\theta = -\delta_\tau(\bar{n}_l/2 + 3/4). Dekomposisi Trotter orde-pertama ini menimbulkan error yang menghilang seiring δτ→0\delta_\tau \to 0. Kita menetapkan δτ=0.0015\delta_\tau = 0.0015 sepanjang tutorial ini.

Hasil dari tiga aproksimasi ini adalah hanya dua qubit fermionik per situs (ni,no)(n_i, n_o) yang dinamis — derajat kebebasan bosonik nln_l telah diserap ke dalam parameter efektif. Ini menghasilkan sirkuit kompak dengan 2N2N qubit untuk NN situs lattice, di mana setiap langkah Trotter memiliki kedalaman gate dua-qubit konstan (13 per langkah).

Apa yang disimulasikan tutorial ini​

Tutorial ini mensimulasikan propagasi hadron: dimulai dari vakum kopling-kuat (state produk), letakkan meson di tengah lattice dan evolusikan seiring waktu. Protokol pengukuran diferensial — menjalankan sirkuit dengan dan tanpa meson di tengah, lalu mengurangkannya — mengisolasi sinyal hadron koheren dari noise hardware maupun efek batas. Hasilnya adalah pola light-cone dari osilasi densitas fermion yang khas dari mode breathing meson terkurung.

Persyaratan​

Sebelum memulai tutorial ini, instal berikut ini:

  • Qiskit SDK v2.0 atau lebih baru, dengan dukungan visualisasi

  • Qiskit Runtime v0.22 atau lebih baru (pip install qiskit-ibm-runtime)

  • Paket Pauli Propagation (pip install pauli-prop)

  • NumPy (pip install numpy)

  • Matplotlib (pip install matplotlib)

Persiapan​

Mulai dengan mengimpor library yang diperlukan dan mendefinisikan fungsi-fungsi bantu yang membangun sirkuit kuantum untuk evolusi waktu LSH. Ada tiga fungsi inti pembangun sirkuit:

  1. pair_hamiltonian_circuit: Mengimplementasikan unitary dua-qubit UIU_I untuk Hamiltonian interaksi aproksimasi antara situs-situs tetangga. Dekomposisi gate-nya adalah: 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: Mengimplementasikan unitary dua-qubit UEU_E untuk energi medan elektrik aproksimasi di setiap situs. Dekomposisi gate-nya adalah: 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: Merakit sirkuit Trotterisasi penuh, menumpuk suku interaksi, elektrik, dan massa dengan gate SWAP untuk mengelola konektivitas 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

Contoh simulator skala-kecil​

Pertama, demonstrasikan alur kerja pada skala kecil menggunakan lattice enam-situs (12 qubit), sehingga kamu bisa memverifikasi konstruksi sirkuit dan memahami observabel fisis sebelum menjalankannya pada hardware.

Langkah 1: Memetakan input klasik ke masalah kuantum​

Definisikan parameter fisis yang sesuai dengan rezim kopling lemah yang dipelajari dalam paper (x=100x = 100, m/g=1m/g = 1). Parameter sirkuit yang diturunkan adalah:

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

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

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

Untuk setiap jumlah langkah Trotter, buat dua sirkuit: satu menginisialisasi meson di tengah (inverse_mid=True) dan satu menyiapkan vakum kopling-kuat (inverse_mid=False). Protokol pengukuran diferensial mengurangkan evolusi vakum untuk mengisolasi sinyal hadron.

# 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

Langkah 2: Optimalkan masalah untuk eksekusi pada hardware kuantum​

Definisikan observabel: pengukuran ZZ satu-qubit pada setiap qubit. Dari ⟨Z⟩\langle Z \rangle kamu bisa mengekstrak probabilitas okupasi lalu bilangan fermion staggered nf(r)n_f(r) pada setiap situs lattice 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

Langkah 3: Eksekusi menggunakan primitif Qiskit​

Gunakan StatevectorEstimator untuk simulasi eksak tanpa noise pada skala kecil.

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

Langkah 4: Post-proses dan kembalikan hasil dalam format klasik yang diinginkan​

Konversikan nilai ekspektasi ke bilangan fermion staggered nf(r,t)n_f(r, t) dan terapkan protokol pengukuran diferensial (meson −- vakum) untuk menghasilkan heatmap propagasi hadron. Ini mereproduksi struktur Gambar 3 dari paper referensi: situs lattice rr pada sumbu-x, langkah Trotter (waktu) tt pada sumbu-y, dan nf(r,t)n_f(r,t) sebagai skala warna.

# 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

Contoh hardware skala-besar​

Sekarang kita naik skala ke lattice 30-situs (60 qubit) pada hardware IBM Quantum. Pada skala ini, sirkuit pada 10 langkah Trotter terdiri dari lebih dari 3400 gate dua-qubit dan 14.000 gate satu-qubit.

Langkah 1-4 (dipadatkan menjadi satu blok kode)​

Aspek-aspek utama dari alur kerja hardware:

  • 10 langkah Trotter untuk sirkuit meson dan vakum (diselang-seling untuk meminimalkan drift)

  • Transpilasi dengan optimization_level=1 — layout sirkuit sudah isomorfik dengan topologi perangkat (rantai linear), sehingga tidak diperlukan SWAP routing. Transpiler digunakan semata-mata untuk memilih rantai qubit fisik ber-noise-rendah dan mendekomposisi gate ke dalam set gate native.

  • EstimatorV2 dengan mitigasi error pembacaan TREX dan Pauli twirling

  • Sesi Batch untuk mengirimkan semua job bersamaan

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

Benchmarking klasik melalui Pauli Propagation​

Metode Pauli Propagation (PPM) menyediakan simulasi klasik tanpa noise dari sirkuit kuantum dengan mem-back-propagate observabel yang diukur melalui sirkuit dalam gambaran Heisenberg. Di bawah lapisan Clifford (gate CNOT, H, S, X), operator Pauli dipetakan ke operator Pauli lain tanpa menambah jumlah suku. Lapisan non-Clifford (gate RzR_z dalam sirkuit) bisa menyebabkan percabangan — dalam kasus terburuk, menggandakan jumlah suku — tapi banyak cabang memiliki koefisien kecil dan bisa dipotong.

Alur kerja dengan pauli-prop adalah:

  1. Pisahkan sirkuit menjadi bagian Clifford dan non-Clifford-nya menggunakan evolve_through_cliffords.

  2. Propagasikan setiap observabel melalui bagian non-Clifford menggunakan propagate_through_circuit, mempertahankan hingga max_terms suku Pauli dan membuang suku dengan koefisien di bawah ambang pemotongan atol.

  3. Evolusikan hasil melalui bagian Clifford menggunakan dukungan Clifford bawaan Qiskit.

  4. Ekstrak nilai ekspektasi dengan menjumlahkan koefisien suku Pauli diagonal (yang hanya berisi II dan ZZ).

Ambang pemotongan​

Parameter atol dalam propagate_through_circuit mengontrol seberapa agresif cabang Pauli kecil dipangkas. Ambang yang sangat ketat (misalnya, 1e-12) mempertahankan hampir semua cabang dan memberi hasil eksak, tapi waktu simulasi meningkat tajam seiring kedalaman sirkuit; simulasi 120-qubit dalam paper memakan waktu sekitar 8,5 jam dengan pengaturan default. Menaikkan ambang (misalnya, ke 1e-6 atau 1e-3) membuang suku-suku yang koefisiennya berada di bawah nilai tersebut, secara dramatis mengurangi jumlah suku yang dilacak dan mempercepat komputasi. Trade-off-nya adalah error aproksimasi kecil yang terkontrol, yang bisa kamu validasi dengan membandingkan hasil pada ambang yang berbeda.

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

Langkah selanjutnya​

Jika kamu merasa karya ini menarik, pertimbangkan untuk menjelajahi materi berikut:

Rekomendasi

Referensi​

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

Sumber: IBM Quantum docs — diperbarui 9 Sep 2026
Versi bahasa Inggris di doQumentation — diperbarui 9 Sep 2026
Masalah situs atau terjemahan?Masalah konten? Edit di IBM Quantum docs