Sari la conținutul principal

Observarea dinamicii hadronice non-abeliene robuste și coerente pe procesoare cuantice zgomotoase

Timp estimat de utilizare: 6 minute pe un procesor Heron (ibm_boston sau echivalent) (NOTĂ: Aceasta este doar o estimare. Timpul tău de execuție poate varia.)

Rezultate ale învățării​

  • Cum pot fi reformulate teoriile gauge de rețea non-abeliene (în mod specific SU(2)) folosind cadrul Loop-String-Hadron (LSH) pentru o simulare cuantică eficientă

  • Cum să construiești circuite de evoluție temporală trotterizate pentru un hamiltonian aproximativ al teoriei gauge SU(2) și cum să le mapezi pe qubiți

  • Cum să rulezi aceste circuite pe hardware IBM Quantum® folosind primitiva Estimator din Qiskit, cu atenuarea erorilor de citire

Cerințe prealabile​

Context​

Motivație​

Cromodinamica cuantică (QCD), teoria gauge SU(3) a forței tari, leagă quarcii în hadroni și guvernează confinarea și ruperea firului. Metodele clasice de QCD pe rețea excelează la proprietățile statice, dar nu pot simula dinamica în timp real din cauza problemei de semn. Calculatoarele cuantice oferă o cale de a ocoli acest obstacol, codificând gradele de libertate ale câmpului gauge direct pe qubiți.

Acest tutorial demonstrează o astfel de simulare: folosește hardware IBM Quantum pentru a simula propagarea în timp real a hadronilor într-o teorie gauge de rețea SU(2) (1+1)-dimensională — cea mai simplă teorie gauge non-abeliană și un pas important către QCD completă.

Hamiltonianul Kogut-Susskind​

Teoria este formulată pe o rețea spațială 1D cu fermioni staggered (materie) pe site-uri și câmpuri gauge SU(2) pe legături. După rescalarea la o formă adimensională, hamiltonianul este:

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

unde HEH_E este energia câmpului cromoelectric, HMH_M este termenul de masă staggered, HIH_I este termenul de interacțiune materie-gauge (hopping), μ=2mgx\mu = 2\frac{m}{g}\sqrt{x} codifică masa fermionului, iar x=1g2a2x = \frac{1}{g^2 a^2} este intensitatea interacțiunii. Limita de continuu a teoriei se află la N→∞N \to \infty și x→∞x \to \infty.

Cadrul Loop-String-Hadron (LSH)​

O provocare esențială este că spațiul Hilbert al câmpului gauge pe fiecare legătură este infinit-dimensional. Cadrul Loop-String-Hadron (LSH) rezolvă acest lucru reformulând teoria în termeni de variabile gauge-invariante — bucle de flux, fire care conectează sarcini separate și hadroni (perechi de fermioni gauge-singlet la un site). În baza LSH, legea lui Gauss este satisfăcută automat prin construcție, astfel încât fiecare stare de bază este fizică. Fiecare site al rețelei este caracterizat de trei numere cuantice (nl,ni,no)(n_l, n_i, n_o) ce reprezintă numărul de buclă, firul de intrare și firul de ieșire, unde ni,no∈{0,1}n_i, n_o \in \{0,1\} sunt fermionice, iar nl≥0n_l \geq 0 este bosonic. Numărul local de fermioni este definit din acestea ca nf(r)=ni(r)+no(r)n_f(r) = n_i(r) + n_o(r) pentru site-urile pare și nf(r)=2−[ni(r)+no(r)]n_f(r) = 2 - [n_i(r) + n_o(r)] pentru site-urile impare.

De la hamiltonianul complet la circuitul cuantic: trei aproximări-cheie​

Circuitul cuantic nu simulează exact hamiltonianul SU(2) complet. În schimb, implementează o serie controlată de aproximări valabile în regimul de cuplaj slab (x≫1x \gg 1). Este esențial să înțelegi ce este și ce nu este aproximat:

Aproximarea 1 — Limita de cuplaj slab pentru HIH_I: Hamiltonianul complet de interacțiune HI(LSH)H_I^{\text{(LSH)}} (Ec. 16 din [1]) conține prefactori care depind de numărul cuantic bosonic nln_l prin termeni de tipul 1/nl+11/\sqrt{n_l+1}. În regimul de cuplaj slab (x≫1x \gg 1), dinamica este dominată de termenul electric HEH_E, care favorizează stările cu nln_l mare. Pentru nl≫1n_l \gg 1, raportul nl/(nl+1)→1n_l/(n_l+1) \to 1 și toți acești prefactori se simplifică la unitate. Hamiltonianul de interacțiune se reduce atunci la un hopping pur local între vecini apropiați:

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

care este independent de nln_l și acționează doar asupra qubiților fermionici (ni,no)(n_i, n_o).

Aproximarea 2 — Fluxul mediat global pentru HEH_E: Energia electrică depinde de nln_l la fiecare legătură. În vidul de cuplaj slab, nln_l este mare și aproximativ uniform. Înlocuiește valorile nln_l dependente de site cu o singură medie globală nˉl\bar{n}_l, transformând HEH_E într-o fază diagonală proporțională cu configurația fermionică la fiecare 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)

unde {r′}\{r'\} însumează peste site-urile din configurația fermionică (ni=0,no=1)(n_i=0, n_o=1), iar hE0h_E^0 este o fază globală pe care o poți ignora.

Aproximarea 3 — Trotterizare: Operatorul de evoluție temporală pentru un pas de durată δτ\delta_\tau este descompus ca:

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

unde c=δτxc = \delta_\tau x, m~=δτμ\tilde{m} = \delta_\tau \mu și θ=−δτ(nˉl/2+3/4)\theta = -\delta_\tau(\bar{n}_l/2 + 3/4). Această descompunere Trotter de ordinul întâi introduce o eroare care se anulează pe măsură ce δτ→0\delta_\tau \to 0. Fixăm δτ=0.0015\delta_\tau = 0.0015 pe tot parcursul.

Rezultatul acestor trei aproximări este că doar cei doi qubiți fermionici per site (ni,no)(n_i, n_o) sunt dinamici — gradul de libertate bosonic nln_l a fost absorbit în parametri efectivi. Aceasta produce un circuit compact cu 2N2N qubiți pentru NN site-uri de rețea, unde fiecare pas Trotter are o adâncime constantă de porți cu doi qubiți (13 per pas).

Ce simulează acest tutorial​

Tutorialul simulează propagarea hadronilor: pornind de la vidul de cuplaj puternic (o stare produs), plasează un meson în centrul rețelei și evoluează în timp. Protocolul de măsurare diferențială — rularea circuitului cu și fără mesonul central, urmată de scădere — izolează semnalul hadronic coerent atât de zgomotul hardware-ului, cât și de efectele de frontieră. Rezultatul este un tipar de con de lumină al oscilațiilor densității de fermioni caracteristic unui mod de respirație al unui meson confinat.

Cerințe​

Înainte de a începe acest tutorial, instalează următoarele:

  • Qiskit SDK v2.0 sau o versiune ulterioară, cu suport pentru vizualizare

  • Qiskit Runtime v0.22 sau o versiune ulterioară (pip install qiskit-ibm-runtime)

  • Pachetul Pauli Propagation (pip install pauli-prop)

  • NumPy (pip install numpy)

  • Matplotlib (pip install matplotlib)

Configurare​

Începe prin a importa bibliotecile necesare și a defini funcțiile ajutătoare care construiesc circuitele cuantice pentru evoluția temporală LSH. Există trei funcții de bază pentru construirea circuitelor:

  1. pair_hamiltonian_circuit: Implementează unitara cu doi qubiți UIU_I pentru hamiltonianul de interacțiune aproximativ dintre site-uri vecine. Descompunerea în porți este: 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: Implementează unitara cu doi qubiți UEU_E pentru energia câmpului electric aproximativ la fiecare site. Descompunerea în porți este: 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: Asamblează circuitul trotterizat complet, stratificând termenii de interacțiune, electric și de masă cu porți SWAP pentru a gestiona conectivitatea qubiților.

# 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

Exemplu de simulator la scară mică​

Mai întâi, demonstrează fluxul de lucru la scară mică folosind o rețea cu șase site-uri (12 qubiți), astfel încât să poți verifica construcția circuitului și să înțelegi observabilele fizice înainte de a rula pe hardware.

Pasul 1: Maparea intrărilor clasice la o problemă cuantică​

Definește parametrii fizici corespunzători regimului de cuplaj slab studiat în articol (x=100x = 100, m/g=1m/g = 1). Parametrii circuitului derivați sunt:

  • c=δτ⋅x=0.15c = \delta_\tau \cdot x = 0.15 (parametrul de interacțiune)

  • θ=−δτ(nˉl/2+3/4)=0.01\theta = -\delta_\tau (\bar{n}_l/2 + 3/4) = 0.01 (faza câmpului electric)

  • m~=δτ⋅μ=0.03\tilde{m} = \delta_\tau \cdot \mu = 0.03 (parametrul de masă)

Pentru fiecare număr de pași Trotter, construiește două circuite: unul care inițializează un meson în centru (inverse_mid=True) și unul care pregătește vidul de cuplaj puternic (inverse_mid=False). Protocolul de măsurare diferențială scade evoluția vidului pentru a izola semnalul hadronic.

# 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

Pasul 2: Optimizarea problemei pentru execuție pe hardware cuantic​

Definește observabilele: măsurători ZZ pe un singur qubit pentru fiecare qubit. Din ⟨Z⟩\langle Z \rangle poți extrage probabilitățile de ocupare și apoi numărul de fermioni staggered nf(r)n_f(r) la fiecare site al rețelei 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

Pasul 3: Execuția folosind primitivele Qiskit​

Folosește StatevectorEstimator pentru o simulare exactă fără zgomot la scară mică.

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

Pasul 4: Post-procesarea și returnarea rezultatului în formatul clasic dorit​

Convertește valorile de așteptare în numărul de fermioni staggered nf(r,t)n_f(r, t) și aplică protocolul de măsurare diferențială (meson −- vid) pentru a produce harta termică a propagării hadronilor. Aceasta reproduce structura Figurii 3 din articolul de referință: site-ul rețelei rr pe axa x, pasul Trotter (timp) tt pe axa y, iar nf(r,t)n_f(r,t) ca scală de culoare.

# 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

Exemplu de hardware la scară mare​

Acum scalăm la o rețea cu 30 de site-uri (60 de qubiți) pe hardware IBM Quantum. La această scară, circuitul cu 10 pași Trotter cuprinde peste 3400 de porți cu doi qubiți și 14.000 de porți cu un singur qubit.

Pașii 1-4 (comprimați într-un singur bloc de cod)​

Aspecte-cheie ale fluxului de lucru pe hardware:

  • 10 pași Trotter pentru circuitele de meson și vid (intercalate pentru derivă minimă)

  • Transpilare cu optimization_level=1 — configurația circuitului este deja izomorfă cu topologia dispozitivului (un lanț liniar), deci nu sunt necesare SWAP-uri de rutare. Transpilerul este folosit doar pentru a selecta un lanț de qubiți fizici cu zgomot redus și pentru a descompune porțile în setul de porți native.

  • EstimatorV2 cu atenuarea erorilor de citire TREX și twirling Pauli

  • Sesiune Batch pentru a trimite toate joburile împreună

# -------------------------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 clasic prin Pauli Propagation​

Metoda Pauli Propagation (PPM) oferă o simulare clasică fără zgomot a circuitului cuantic, propagând înapoi observabilele măsurate prin circuit în imaginea Heisenberg. Sub straturi Clifford (porți CNOT, H, S, X), operatorii Pauli se mapează în alți operatori Pauli fără a crește numărul de termeni. Straturile non-Clifford (porțile RzR_z din circuit) pot cauza ramificare — în cel mai rău caz, dublând numărul de termeni — dar multe ramuri au coeficienți mici și pot fi trunchiate.

Fluxul de lucru cu pauli-prop este:

  1. Împarte circuitul în părțile sale Clifford și non-Clifford folosind evolve_through_cliffords.

  2. Propagă fiecare observabilă prin partea non-Clifford folosind propagate_through_circuit, păstrând până la max_terms termeni Pauli și eliminând termenii cu coeficienți sub pragul de trunchiere atol.

  3. Evoluează rezultatul prin partea Clifford folosind suportul Clifford integrat al Qiskit.

  4. Extrage valoarea de așteptare însumând coeficienții termenilor Pauli diagonali (care conțin doar II și ZZ).

Pragul de trunchiere​

Parametrul atol din propagate_through_circuit controlează cât de agresiv sunt eliminate ramurile Pauli mici. Un prag foarte strict (de exemplu, 1e-12) păstrează aproape toate ramurile și oferă rezultate exacte, dar timpul de simulare crește abrupt odată cu adâncimea circuitului; simularea cu 120 de qubiți din articol a durat aproximativ 8,5 ore cu setările implicite. Ridicarea pragului (de exemplu, la 1e-6 sau 1e-3) elimină termenii ai căror coeficienți sunt sub acea valoare, reducând drastic numărul de termeni urmăriți și accelerând calculul. Compromisul este o eroare de aproximare mică, controlabilă, pe care o poți valida comparând rezultatele la praguri diferite.

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

Pași următori​

Dacă acest material ți s-a părut interesant, ia în considerare explorarea următoarelor resurse:

Recomandări

Referințe​

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

Sursă: IBM Quantum docs — actualizat 9 sept. 2026
Versiunea în engleză pe doQumentation — actualizată 9 sept. 2026
Problemă cu site-ul sau cu traducerea?Problemă de conținut? Editează pe IBM Quantum docs