Langkau ke kandungan utama

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

Anggaran penggunaan: 6 minit pada pemproses Heron (ibm_boston atau setara) (NOTA: Ini hanyalah anggaran. Masa jalanan anda mungkin berbeza.)

Hasil pembelajaran​

  • Bagaimana teori gauge kekisi bukan-Abelian (khususnya SU(2)) boleh dirumuskan semula menggunakan rangka kerja Loop-String-Hadron (LSH) untuk simulasi kuantum yang cekap

  • Bagaimana membina Circuit evolusi masa Trotterized untuk Hamiltonian teori gauge SU(2) anggaran dan memetakannya kepada Qubit

  • Bagaimana menjalankan Circuit ini pada perkakasan IBM Quantum® menggunakan primitif Estimator Qiskit dengan mitigasi ralat pembacaan

Prerequisites​

Background​

Motivation​

Kromodinamik Kuantum (QCD), teori gauge SU(3) bagi daya kuat, mengikat quark menjadi hadron dan mengawal kurungan dan pemecahan tali. Kaedah QCD kekisi klasik cemerlang dalam sifat statik tetapi tidak dapat mensimulasikan dinamik masa nyata disebabkan masalah tanda. Komputer kuantum menawarkan laluan untuk mengatasi halangan ini dengan mengekod darjah kebebasan medan gauge secara langsung ke atas Qubit.

Tutorial ini menunjukkan simulasi sedemikian: menggunakan perkakasan IBM Quantum untuk mensimulasikan perambatan hadron masa nyata dalam teori gauge kekisi SU(2) (1+1)-dimensi — teori gauge bukan-Abelian yang paling mudah dan batu loncatan ke arah QCD penuh.

Hamiltonian Kogut-Susskind​

Teori ini dirumuskan pada kekisi ruang 1D dengan fermion tersusun-selang-seli (jirim) pada tapak dan medan gauge SU(2) pada pautan. Selepas penskalaan semula kepada bentuk tanpa dimensi, Hamiltonian ialah:

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

di mana HEH_E ialah tenaga medan kromoelektrik, HMH_M ialah sebutan jisim tersusun-selang-seli, HIH_I ialah sebutan interaksi jirim-gauge (hopping), μ=2mgx\mu = 2\frac{m}{g}\sqrt{x} mengekod jisim fermion, dan x=1g2a2x = \frac{1}{g^2 a^2} ialah kekuatan interaksi. Had kontinum bagi teori ini terletak pada N→∞N \to \infty dan x→∞x \to \infty.

The Loop-String-Hadron (LSH) framework​

Satu cabaran utama ialah ruang Hilbert medan gauge pada setiap pautan adalah tak terhingga dimensi. Rangka kerja Loop-String-Hadron (LSH) menangani ini dengan merumuskan semula teori dari segi pemboleh ubah tak varian-gauge — gelung fluks, tali yang menghubungkan cas terpisah, dan hadron (pasangan fermion singlet-gauge pada tapak). Dalam asas LSH, hukum Gauss dipenuhi secara automatik melalui pembinaan, jadi setiap keadaan asas adalah fizikal. Setiap tapak kekisi dicirikan oleh tiga nombor kuantum (nl,ni,no)(n_l, n_i, n_o) yang mewakili nombor gelung, tali masuk, dan tali keluar, di mana ni,no∈{0,1}n_i, n_o \in \{0,1\} adalah fermionik dan nl≥0n_l \geq 0 adalah bosonik. Nombor fermion setempat ditakrifkan daripada ini sebagai nf(r)=ni(r)+no(r)n_f(r) = n_i(r) + n_o(r) untuk tapak genap dan nf(r)=2−[ni(r)+no(r)]n_f(r) = 2 - [n_i(r) + n_o(r)] untuk tapak ganjil.

Daripada Hamiltonian penuh kepada litar kuantum: tiga penghampiran utama​

Circuit kuantum tidak mensimulasikan Hamiltonian SU(2) penuh dengan tepat. Sebaliknya, ia mengimplementasikan siri anggaran terkawal yang sah dalam rejim gandingan lemah (x≫1x \gg 1). Memahami apa yang dianggarkan dan tidak dianggarkan adalah penting:

Anggaran 1 — Had gandingan lemah untuk HIH_I: Hamiltonian interaksi penuh HI(LSH)H_I^{\text{(LSH)}} (Pers. 16 dalam [1]) mengandungi prapemberat yang bergantung kepada nombor kuantum bosonik nln_l melalui sebutan seperti 1/nl+11/\sqrt{n_l+1}. Dalam rejim gandingan lemah (x≫1x \gg 1), dinamik didominasi oleh sebutan elektrik HEH_E, yang menggemari keadaan dengan nln_l yang besar. Untuk nl≫1n_l \gg 1, nisbah nl/(nl+1)→1n_l/(n_l+1) \to 1 dan semua prapemberat ini dipermudahkan kepada satu. Hamiltonian interaksi kemudian tereduksi kepada hopping jiran-terdekat setempat semata-mata:

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 kepada nln_l dan hanya bertindak pada Qubit fermionik (ni,no)(n_i, n_o).

Anggaran 2 — Fluks-purata-global untuk HEH_E: Tenaga elektrik bergantung kepada nln_l pada setiap pautan. Dalam vakum gandingan lemah, nln_l adalah besar dan lebih kurang seragam. Gantikan nilai nln_l yang bergantung tapak dengan satu purata global tunggal nˉl\bar{n}_l, menjadikan HEH_E satu fasa pepenjuru berkadar dengan konfigurasi fermion pada setiap tapak:

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 ke atas tapak dalam konfigurasi fermionik (ni=0,no=1)(n_i=0, n_o=1), dan hE0h_E^0 ialah fasa global yang boleh anda abaikan.

Anggaran 3 — Trotterization: Operator evolusi masa untuk satu langkah bertempoh δτ\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). Penguraian Trotter tertib pertama ini memperkenalkan ralat yang hilang apabila δτ→0\delta_\tau \to 0. Kita menetapkan δτ=0.0015\delta_\tau = 0.0015 sepanjang.

Hasil ketiga-tiga anggaran ini ialah hanya kedua-dua Qubit fermionik setiap tapak (ni,no)(n_i, n_o) yang dinamik — darjah kebebasan bosonik nln_l telah diserap ke dalam parameter berkesan. Ini menghasilkan Circuit padat dengan 2N2N Qubit untuk NN tapak kekisi, di mana setiap langkah Trotter mempunyai kedalaman gate dua-qubit tetap (13 setiap langkah).

Apa yang disimulasikan oleh tutorial ini​

Tutorial ini mensimulasikan perambatan hadron: bermula daripada vakum gandingan kuat (keadaan hasil darab), letakkan satu meson di tengah kekisi dan berevolusi dalam masa. Protokol pengukuran pembezaan — menjalankan Circuit dengan dan tanpa meson tengah, kemudian menolak — mengasingkan isyarat hadron koheren daripada kedua-dua hingar perkakasan dan kesan sempadan. Hasilnya ialah corak kon-cahaya bagi ayunan ketumpatan fermion yang menjadi ciri mod pernafasan meson terkurung.

Requirements​

Sebelum memulakan tutorial ini, pasang yang berikut:

  • Qiskit SDK v2.0 atau lebih baharu, dengan sokongan visualisasi

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

  • Pakej Pauli Propagation (pip install pauli-prop)

  • NumPy (pip install numpy)

  • Matplotlib (pip install matplotlib)

Setup​

Mulakan dengan mengimport pustaka yang diperlukan dan mentakrifkan fungsi pembantu yang membina Circuit kuantum untuk evolusi masa LSH. Terdapat tiga fungsi pembinaan Circuit teras:

  1. pair_hamiltonian_circuit: Mengimplementasikan unitari dua-qubit UIU_I untuk Hamiltonian interaksi anggaran antara tapak bersebelahan. Penguraian gate ialah: 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 unitari dua-qubit UEU_E untuk tenaga medan elektrik anggaran pada setiap tapak. Penguraian gate ialah: 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: Menghimpunkan Circuit Trotterized penuh, melapiskan sebutan interaksi, elektrik, dan jisim dengan gate SWAP untuk menguruskan konektiviti 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 berskala kecil​

Pertama, tunjukkan aliran kerja pada skala kecil menggunakan kekisi enam-tapak (12 Qubit), supaya anda boleh mengesahkan pembinaan Circuit dan memahami cerapan fizikal sebelum menjalankan pada perkakasan.

Langkah 1: Petakan input klasik kepada masalah kuantum​

Takrifkan parameter fizikal yang sepadan dengan rejim gandingan lemah yang dikaji dalam kertas kerja (x=100x = 100, m/g=1m/g = 1). Parameter Circuit terbitan ialah:

  • 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 (fasa medan elektrik)

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

Untuk setiap kiraan langkah Trotter, bina dua Circuit: satu memulakan meson di tengah (inverse_mid=True) dan satu menyediakan vakum gandingan kuat (inverse_mid=False). Protokol pengukuran pembezaan menolak evolusi vakum untuk mengasingkan isyarat 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: Optimumkan masalah untuk pelaksanaan pada perkakasan kuantum​

Takrifkan cerapan: pengukuran ZZ satu-qubit pada setiap Qubit. Daripada ⟨Z⟩\langle Z \rangle anda boleh mengekstrak kebarangkalian pekerjaan dan kemudian nombor fermion tersusun-selang-seli nf(r)n_f(r) pada setiap tapak kekisi 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

Step 3: Execute using Qiskit primitives​

Gunakan StatevectorEstimator untuk simulasi tanpa hingar yang tepat 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: Pasca-proses dan kembalikan hasil dalam format klasik yang dikehendaki​

Tukarkan nilai jangkaan kepada nombor fermion tersusun-selang-seli nf(r,t)n_f(r, t) dan aplikasikan protokol pengukuran pembezaan (meson −- vakum) untuk menghasilkan peta haba perambatan hadron. Ini menghasilkan semula struktur Rajah 3 daripada kertas rujukan: tapak kekisi rr pada paksi-x, langkah Trotter (masa) tt pada paksi-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 perkakasan berskala besar​

Kita kini menskalakan kepada kekisi 30-tapak (60 Qubit) pada perkakasan IBM Quantum. Pada skala ini, Circuit pada 10 langkah Trotter mengandungi lebih daripada 3400 gate dua-qubit dan 14,000 gate satu-qubit.

Langkah 1-4 (dimampatkan ke dalam satu blok kod)​

Aspek utama aliran kerja perkakasan:

  • 10 langkah Trotter untuk Circuit meson dan vakum (diselang-seli untuk hanyutan minimum)

  • Transpilation dengan optimization_level=1 — susun atur Circuit sudah isomorfik kepada topologi peranti (rantai linear), jadi tiada SWAP penghalaan diperlukan. Transpiler digunakan semata-mata untuk memilih rantai Qubit fizikal hingar-rendah dan menguraikan gate ke dalam set gate asli.

  • EstimatorV2 dengan mitigasi ralat pembacaan TREX dan Pauli twirling

  • Sesi Batch untuk menghantar semua jobs bersama

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

Penandaarasan klasik melalui Perambatan Pauli​

Kaedah Pauli Propagation (PPM) menyediakan simulasi klasik tanpa hingar bagi Circuit kuantum dengan merambatkan semula cerapan yang diukur melalui Circuit dalam gambaran Heisenberg. Di bawah lapisan Clifford (gate CNOT, H, S, X), operator Pauli dipetakan kepada operator Pauli lain tanpa meningkatkan bilangan sebutan. Lapisan bukan-Clifford (gate RzR_z dalam Circuit) boleh menyebabkan percabangan — dalam kes terburuk, menggandakan bilangan sebutan — tetapi banyak cabang mempunyai pekali kecil dan boleh dipotong.

Aliran kerja dengan pauli-prop ialah:

  1. Pisahkan Circuit kepada bahagian Clifford dan bukan-Cliffordnya menggunakan evolve_through_cliffords.

  2. Rambatkan setiap cerapan melalui bahagian bukan-Clifford menggunakan propagate_through_circuit, mengekalkan sehingga max_terms sebutan Pauli dan membuang sebutan dengan pekali di bawah ambang pemotongan atol.

  3. Evolusikan hasil melalui bahagian Clifford menggunakan sokongan Clifford terbina dalam Qiskit.

  4. Ekstrak nilai jangkaan dengan menjumlahkan pekali sebutan Pauli pepenjuru (yang hanya mengandungi II dan ZZ).

Ambang pemotongan​

Parameter atol dalam propagate_through_circuit mengawal betapa agresifnya cabang Pauli kecil dipangkas. Ambang yang sangat ketat (contohnya, 1e-12) mengekalkan hampir semua cabang dan memberikan hasil yang tepat, tetapi masa simulasi berkembang dengan curam mengikut kedalaman Circuit; simulasi 120-qubit dalam kertas kerja mengambil masa lebih kurang 8.5 jam dengan tetapan lalai. Menaikkan ambang (contohnya, kepada 1e-6 atau 1e-3) membuang sebutan yang pekalinya jatuh di bawah nilai tersebut, mengurangkan bilangan sebutan yang dijejak secara drastik dan mempercepatkan pengiraan. Pertukaran nilainya ialah ralat anggaran yang kecil dan terkawal yang boleh anda sahkan dengan membandingkan hasil pada ambang yang berbeza.

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 seterusnya​

Jika kamu rasa kerja ini menarik, pertimbangkan untuk meneroka bahan berikut:

Cadangan

Rujukan​

[1] Kertas asal: 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 — dikemas kini 9 Sep 2026
Versi bahasa Inggeris di doQumentation — dikemas kini 9 Sep 2026
Masalah laman web atau terjemahan?Masalah kandungan? Sunting di IBM Quantum docs