Lumaktaw sa pangunahing nilalaman

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

Tinatayang paggamit: 6 na minuto sa isang Heron processor (ibm_boston o katumbas) (TANDAAN: Tantiya lamang ito. Maaaring mag-iba ang iyong runtime.)

Mga learning outcome​

  • Kung paano maaaring i-reformulate ang non-Abelian lattice gauge theories (partikular ang SU(2)) gamit ang Loop-String-Hadron (LSH) framework para sa episyenteng quantum simulation

  • Kung paano bumuo ng Trotterized time-evolution circuits para sa isang approximate SU(2) gauge theory Hamiltonian at i-map ang mga ito sa qubits

  • Kung paano patakbuhin ang mga circuit na ito sa IBM Quantum® hardware gamit ang Qiskit Estimator primitive na may readout error mitigation

Prerequisites​

Background​

Motivation​

Ang Quantum Chromodynamics (QCD), ang SU(3) gauge theory ng malakas na puwersa, ay nagbubuklod sa mga quark tungo sa hadrons at pinamamahalaan ang confinement at string breaking. Ang mga classical lattice QCD method ay mahusay sa static properties ngunit hindi maaaring gayahin ang real-time dynamics dahil sa sign problem. Nag-aalok ang mga quantum computer ng ruta upang malampasan ang hadlang na ito sa pamamagitan ng direktang pag-encode ng gauge field degrees of freedom sa qubits.

Ipinapakita ng tutorial na ito ang ganitong simulation: gamitin ang IBM Quantum hardware upang gayahin ang real-time hadron propagation sa isang (1+1)-dimensional SU(2) lattice gauge theory — ang pinakasimpleng non-Abelian gauge theory at isang hakbang patungo sa ganap na QCD.

Ang Kogut-Susskind Hamiltonian​

Ang teorya ay binubuo sa isang 1D spatial lattice na may staggered fermions (matter) sa mga site at SU(2) gauge fields sa mga link. Pagkatapos ng pag-rescale sa dimensionless form, ang Hamiltonian ay:

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

kung saan ang HEH_E ay ang chromoelectric field energy, ang HMH_M ay ang staggered mass term, ang HIH_I ay ang matter-gauge interaction (hopping) term, ang μ=2mgx\mu = 2\frac{m}{g}\sqrt{x} ay nag-encode ng fermion mass, at ang x=1g2a2x = \frac{1}{g^2 a^2} ay ang interaction strength. Ang continuum limit ng teorya ay nasa N→∞N \to \infty at x→∞x \to \infty.

Ang Loop-String-Hadron (LSH) na framework​

Isang pangunahing hamon ay ang gauge field Hilbert space sa bawat link ay infinite-dimensional. Tinutugunan ito ng Loop-String-Hadron (LSH) framework sa pamamagitan ng pag-reformulate ng teorya sa mga terminong gauge-invariant na variable — mga loop ng flux, mga string na nag-uugnay sa mga hiwalay na charge, at mga hadron (gauge-singlet fermion pairs sa isang site). Sa LSH basis, ang batas ni Gauss ay natutugunan nang awtomatiko sa pamamagitan ng konstruksyon, kaya bawat basis state ay physical. Ang bawat lattice site ay kinakarakter ng tatlong quantum number (nl,ni,no)(n_l, n_i, n_o) na kumakatawan sa loop number, incoming string, at outgoing string, kung saan ang ni,no∈{0,1}n_i, n_o \in \{0,1\} ay fermionic at ang nl≥0n_l \geq 0 ay bosonic. Ang local fermion number ay idinidepine mula rito bilang nf(r)=ni(r)+no(r)n_f(r) = n_i(r) + n_o(r) para sa even sites at nf(r)=2−[ni(r)+no(r)]n_f(r) = 2 - [n_i(r) + n_o(r)] para sa odd sites.

Mula sa buong Hamiltonian patungo sa quantum circuit: tatlong pangunahing approximation​

Ang quantum circuit ay hindi eksaktong ginagaya ang buong SU(2) Hamiltonian. Sa halip, ipinapatupad nito ang isang kontroladong serye ng mga approximation na balido sa weak-coupling regime (x≫1x \gg 1). Mahalagang maintindihan kung ano ang na-approximate at kung ano ang hindi:

Approximation 1 — Weak-coupling limit para sa HIH_I: Ang buong interaction Hamiltonian HI(LSH)H_I^{\text{(LSH)}} (Eq. 16 sa [1]) ay naglalaman ng mga prefactor na nakadepende sa bosonic quantum number nln_l sa pamamagitan ng mga terminong tulad ng 1/nl+11/\sqrt{n_l+1}. Sa weak-coupling regime (x≫1x \gg 1), ang dynamics ay dominado ng electric term HEH_E, na pumapabor sa mga state na may malaking nln_l. Para sa nl≫1n_l \gg 1, ang ratio nl/(nl+1)→1n_l/(n_l+1) \to 1 at lahat ng mga prefactor na ito ay nagsi-simplify sa unity. Ang interaction Hamiltonian ay babagsak sa isang purely local nearest-neighbor 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],

na independiyente sa nln_l at kumikilos lamang sa fermionic (ni,no)(n_i, n_o) qubits.

Approximation 2 — Global-average-flux para sa HEH_E: Ang electric energy ay nakadepende sa nln_l sa bawat link. Sa weak-coupling vacuum, ang nln_l ay malaki at halos uniporme. Palitan ang site-dependent na mga value ng nln_l ng isang solong global average nˉl\bar{n}_l, na ginagawang diagonal phase ang HEH_E na proporsyonal sa fermion configuration sa bawat site:

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)

kung saan ang {r′}\{r'\} ay nagsu-sum sa mga site sa fermionic configuration (ni=0,no=1)(n_i=0, n_o=1), at ang hE0h_E^0 ay isang global phase na maaari mong balewalain.

Approximation 3 — Trotterization: Ang time-evolution operator para sa isang hakbang na may tagal na δτ\delta_\tau ay na-decompose bilang:

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

kung saan ang c=δτxc = \delta_\tau x, m~=δτμ\tilde{m} = \delta_\tau \mu, at θ=−δτ(nˉl/2+3/4)\theta = -\delta_\tau(\bar{n}_l/2 + 3/4). Ang first-order Trotter decomposition na ito ay nagpapakilala ng error na nawawala habang δτ→0\delta_\tau \to 0. Itinatakda natin ang δτ=0.0015\delta_\tau = 0.0015 sa buong kabuuan.

Ang resulta ng tatlong approximation na ito ay tanging ang dalawang fermionic qubit bawat site (ni,no)(n_i, n_o) ang dynamical — ang bosonic nln_l degree of freedom ay na-absorb na sa effective parameters. Nagbubunga ito ng compact circuit na may 2N2N qubits para sa NN lattice sites, kung saan ang bawat Trotter step ay may constant two-qubit gate depth (13 bawat hakbang).

Ang isinisimulate ng tutorial na ito​

Ginagaya ng tutorial na ito ang hadron propagation: simula sa strong-coupling vacuum (isang product state), maglagay ng meson sa gitna ng lattice at ebolusyunin sa oras. Ang differential measurement protocol — pagpapatakbo ng circuit na may at walang central meson, pagkatapos ay pagbabawas — ay ihinihiwalay ang coherent hadron signal mula sa parehong hardware noise at boundary effects. Ang resulta ay isang light-cone pattern ng fermion density oscillations na katangian ng confined meson breathing mode.

Requirements​

Bago simulan ang tutorial na ito, i-install ang mga sumusunod:

  • Qiskit SDK v2.0 o mas bago, na may visualization support

  • Qiskit Runtime v0.22 o mas bago (pip install qiskit-ibm-runtime)

  • Pauli Propagation package (pip install pauli-prop)

  • NumPy (pip install numpy)

  • Matplotlib (pip install matplotlib)

Setup​

Simulan sa pag-import ng mga kailangang library at pagdedefine ng mga helper function na bumubuo ng quantum circuits para sa LSH time evolution. May tatlong core circuit-building function:

  1. pair_hamiltonian_circuit: Ipinapatupad ang two-qubit unitary UIU_I para sa approximate interaction Hamiltonian sa pagitan ng mga karatig na site. Ang gate decomposition ay: 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: Ipinapatupad ang two-qubit unitary UEU_E para sa approximate electric field energy sa bawat site. Ang gate decomposition ay: 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: Isinasama ang buong Trotterized circuit, sa pamamagitan ng pag-layer ng interaction, electric, at mass terms na may SWAP gates upang pamahalaan ang qubit connectivity.

# 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

Halimbawa ng maliitang simulator​

Una, ipakita ang workflow sa maliit na sukat gamit ang isang six-site lattice (12 qubits), upang ma-verify mo ang circuit construction at maintindihan ang physical observables bago patakbuhin sa hardware.

Hakbang 1: I-map ang mga classical input sa isang quantum na problema​

Idepine ang mga physical parameter na tumutugma sa weak-coupling regime na pinag-aralan sa papel (x=100x = 100, m/g=1m/g = 1). Ang mga derived circuit parameter ay:

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

  • θ=−δτ(nˉl/2+3/4)=0.01\theta = -\delta_\tau (\bar{n}_l/2 + 3/4) = 0.01 (phase ng electric field)

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

Para sa bawat bilang ng Trotter step, bumuo ng dalawang circuit: isang nag-iinitialize ng meson sa gitna (inverse_mid=True) at isa na naghahanda ng strong-coupling vacuum (inverse_mid=False). Ang differential measurement protocol ay ibinabawas ang vacuum evolution upang ihiwalay ang hadron signal.

# 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

Hakbang 2: I-optimize ang problema para sa pagpapatupad sa quantum hardware​

Idepine ang mga observable: single-qubit ZZ measurements sa bawat qubit. Mula sa ⟨Z⟩\langle Z \rangle maaari mong i-extract ang occupation probabilities at pagkatapos ang staggered fermion number nf(r)n_f(r) sa bawat lattice site 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

Hakbang 3: Isagawa gamit ang Qiskit primitives​

Gamitin ang StatevectorEstimator para sa exact noiseless simulation sa maliit na sukat.

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

Hakbang 4: I-post-process at ibalik ang resulta sa nais na classical format​

I-convert ang mga expectation value sa staggered fermion number nf(r,t)n_f(r, t) at ilapat ang differential measurement protocol (meson −- vacuum) upang makabuo ng hadron propagation heatmap. Ito ay muling gumagawa ng structure ng Figure 3 mula sa reference paper: lattice site rr sa x-axis, Trotter step (oras) tt sa y-axis, at nf(r,t)n_f(r,t) bilang color scale.

# 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

Halimbawa ng malakihang hardware​

Ngayon ay maglalaki tayo sa isang 30-site lattice (60 qubits) sa IBM Quantum hardware. Sa sukat na ito, ang circuit sa 10 Trotter steps ay binubuo ng higit sa 3400 two-qubit gates at 14,000 single-qubit gates.

Mga Hakbang 1-4 (na pinagsama sa iisang code block)​

Mga pangunahing aspeto ng hardware workflow:

  • 10 Trotter steps para sa meson at vacuum circuits (interleaved para sa minimal drift)

  • Transpilation na may optimization_level=1 — ang circuit layout ay isomorphic na sa device topology (isang linear chain), kaya walang kailangang routing SWAPs. Ginagamit lamang ang transpiler upang pumili ng low-noise chain ng physical qubits at i-decompose ang mga gate sa native gate set.

  • EstimatorV2 na may TREX readout error mitigation at Pauli twirling

  • Batch session upang isumite ang lahat ng jobs nang magkasama

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

Classical benchmarking gamit ang Pauli Propagation​

Ang Pauli Propagation Method (PPM) ay nagbibigay ng noiseless classical simulation ng quantum circuit sa pamamagitan ng back-propagating ng mga nasukat na observable sa buong circuit sa Heisenberg picture. Sa ilalim ng Clifford layers (CNOT, H, S, X gates), ang Pauli operators ay nagma-map sa iba pang Pauli operators nang hindi dinaragdagan ang bilang ng mga terms. Ang Non-Clifford layers (ang mga RzR_z gate sa circuit) ay maaaring magdulot ng branching — sa pinakamasamang kaso, dinodoble ang bilang ng mga terms — ngunit maraming branches ang may maliliit na coefficient at maaaring i-truncate.

Ang workflow gamit ang pauli-prop ay:

  1. Hatiin ang circuit sa Clifford at non-Clifford na bahagi gamit ang evolve_through_cliffords.

  2. I-propagate ang bawat observable sa non-Clifford na bahagi gamit ang propagate_through_circuit, na iniingatan ang hanggang max_terms na Pauli terms at binabalewala ang mga terms na may coefficient na mas mababa sa truncation threshold na atol.

  3. Ebolusyunin ang resulta sa Clifford na bahagi gamit ang built-in Clifford support ng Qiskit.

  4. I-extract ang expectation value sa pamamagitan ng pagsuma ng mga coefficient ng diagonal Pauli terms (naglalaman lamang ng II at ZZ).

Threshold ng truncation​

Kinokontrol ng parameter na atol sa propagate_through_circuit kung gaano ka-agresibo ang pagpuputol ng maliliit na Pauli branches. Ang isang napakahigpit na threshold (halimbawa, 1e-12) ay iniingatan halos lahat ng branches at nagbibigay ng eksaktong resulta, ngunit ang simulation time ay lumalaki nang malaki kasabay ng circuit depth; ang 120-qubit simulation sa papel ay umabot ng humigit-kumulang 8.5 oras gamit ang default settings. Ang pagtaas ng threshold (halimbawa, sa 1e-6 o 1e-3) ay tinatanggal ang mga terms na ang coefficient ay bumagsak sa ibaba ng value na iyon, dramatikong binabawasan ang bilang ng mga na-track na terms at pinapabilis ang computation. Ang trade-off ay isang maliit, makokontrol na approximation error na maaari mong i-validate sa pamamagitan ng paghahambing ng mga resulta sa iba't ibang threshold.

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

Mga susunod na hakbang​

Kung nakita mong interesante ang gawaing ito, isaalang-alang ang paggalugad sa mga sumusunod na materyal:

Mga Rekomendasyon

References​

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

Pinagmulan: IBM Quantum docs — na-update noong Set 9, 2026
Bersyong Ingles sa doQumentation — na-update noong Set 9, 2026
Isyu sa site o pagsasalin?Isyu sa nilalaman? I-edit sa IBM Quantum docs