Vai al contenuto principale

Osservazione di dinamica adronica non-Abeliana robusta e coerente su processori quantistici rumorosi

Utilizzo stimato: 6 minuti su un processore Heron (ibm_boston o equivalente) (NOTA: Questa è solo una stima. Il tempo di esecuzione potrebbe variare.)

Risultati di apprendimento​

  • Come le teorie di gauge non-Abeliane (nello specifico SU(2)) possano essere riformulate usando il framework Loop-String-Hadron (LSH) per una simulazione quantistica efficiente

  • Come costruire circuiti di evoluzione temporale Trotterizzati per un'Hamiltoniana approssimata di teoria di gauge SU(2) e mapparli su qubit

  • Come eseguire questi circuiti su hardware IBM Quantum® usando la primitiva Qiskit Estimator con mitigazione degli errori di lettura

Prerequisiti​

Contesto​

Motivazione​

La Cromodinamica Quantistica (QCD), la teoria di gauge SU(3) della forza forte, lega i quark in adroni e governa il confinamento e la rottura delle stringhe. I metodi classici di QCD reticolare eccellono nelle proprietà statiche ma non possono simulare la dinamica in tempo reale a causa del problema del segno. I computer quantistici offrono una via per aggirare questa barriera codificando i gradi di libertà del campo di gauge direttamente sui qubit.

Questo tutorial dimostra una simulazione di questo tipo: usa hardware IBM Quantum per simulare la propagazione adronica in tempo reale in una teoria di gauge reticolare SU(2) in (1+1) dimensioni — la più semplice teoria di gauge non-Abeliana e un passo intermedio verso la QCD completa.

L'Hamiltoniana di Kogut-Susskind​

La teoria è formulata su un reticolo spaziale 1D con fermioni staggered (materia) sui siti e campi di gauge SU(2) sui link. Dopo il riscalamento in forma adimensionale, l'Hamiltoniana è:

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

dove HEH_E è l'energia del campo cromoelettrico, HMH_M è il termine di massa staggered, HIH_I è il termine di interazione materia-gauge (hopping), μ=2mgx\mu = 2\frac{m}{g}\sqrt{x} codifica la massa del fermione, e x=1g2a2x = \frac{1}{g^2 a^2} è la forza di interazione. Il limite del continuo della teoria si trova a N→∞N \to \infty e x→∞x \to \infty.

Il framework Loop-String-Hadron (LSH)​

Una sfida chiave è che lo spazio di Hilbert del campo di gauge su ciascun link è di dimensione infinita. Il framework Loop-String-Hadron (LSH) affronta questo problema riformulando la teoria in termini di variabili gauge-invarianti — anelli di flusso, stringhe che collegano cariche separate, e adroni (coppie di fermioni gauge-singoletto su un sito). Nella base LSH, la legge di Gauss è soddisfatta automaticamente per costruzione, quindi ogni stato di base è fisico. Ogni sito del reticolo è caratterizzato da tre numeri quantici (nl,ni,no)(n_l, n_i, n_o) che rappresentano il numero di anelli, la stringa entrante e la stringa uscente, dove ni,no∈{0,1}n_i, n_o \in \{0,1\} sono fermionici e nl≥0n_l \geq 0 è bosonico. Il numero fermionico locale è definito da questi come nf(r)=ni(r)+no(r)n_f(r) = n_i(r) + n_o(r) per i siti pari e nf(r)=2−[ni(r)+no(r)]n_f(r) = 2 - [n_i(r) + n_o(r)] per i siti dispari.

Dall'Hamiltoniana completa al circuito quantistico: tre approssimazioni chiave​

Il circuito quantistico non simula esattamente l'Hamiltoniana SU(2) completa. Invece, implementa una serie controllata di approssimazioni valide nel regime di accoppiamento debole (x≫1x \gg 1). Comprendere cosa è e cosa non è approssimato è essenziale:

Approssimazione 1 — Limite di accoppiamento debole per HIH_I: L'Hamiltoniana di interazione completa HI(LSH)H_I^{\text{(LSH)}} (Eq. 16 in [1]) contiene prefattori che dipendono dal numero quantico bosonico nln_l tramite termini come 1/nl+11/\sqrt{n_l+1}. Nel regime di accoppiamento debole (x≫1x \gg 1), la dinamica è dominata dal termine elettrico HEH_E, che favorisce stati con nln_l grande. Per nl≫1n_l \gg 1, il rapporto nl/(nl+1)→1n_l/(n_l+1) \to 1 e tutti questi prefattori si semplificano a unità. L'Hamiltoniana di interazione si riduce quindi a un hopping puramente locale tra primi vicini:

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

che è indipendente da nln_l e agisce solo sui qubit fermionici (ni,no)(n_i, n_o).

Approssimazione 2 — Flusso medio globale per HEH_E: L'energia elettrica dipende da nln_l a ciascun link. Nel vuoto di accoppiamento debole, nln_l è grande e approssimativamente uniforme. Sostituiamo i valori di nln_l dipendenti dal sito con una singola media globale nˉl\bar{n}_l, rendendo HEH_E una fase diagonale proporzionale alla configurazione fermionica a ciascun sito:

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)

dove {r′}\{r'\} somma sui siti nella configurazione fermionica (ni=0,no=1)(n_i=0, n_o=1), e hE0h_E^0 è una fase globale che puoi ignorare.

Approssimazione 3 — Trotterizzazione: L'operatore di evoluzione temporale per un passo di durata δτ\delta_\tau è decomposto come:

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

dove c=δτxc = \delta_\tau x, m~=δτμ\tilde{m} = \delta_\tau \mu, e θ=−δτ(nˉl/2+3/4)\theta = -\delta_\tau(\bar{n}_l/2 + 3/4). Questa decomposizione di Trotter del primo ordine introduce un errore che svanisce quando δτ→0\delta_\tau \to 0. Fissiamo δτ=0.0015\delta_\tau = 0.0015 per tutto il tutorial.

Il risultato di queste tre approssimazioni è che solo i due qubit fermionici per sito (ni,no)(n_i, n_o) sono dinamici — il grado di libertà bosonico nln_l è stato assorbito in parametri effettivi. Questo produce un circuito compatto con 2N2N qubit per NN siti del reticolo, dove ogni passo di Trotter ha una profondità costante di porte a due qubit (13 per passo).

Cosa simula questo tutorial​

Il tutorial simula la propagazione adronica: partendo dal vuoto di accoppiamento forte (uno stato prodotto), viene posizionato un mesone al centro del reticolo ed evoluto nel tempo. Il protocollo di misurazione differenziale — eseguendo il circuito con e senza il mesone centrale, e poi sottraendo — isola il segnale coerente dell'adrone sia dal rumore hardware sia dagli effetti di bordo. Il risultato è un pattern a cono di luce di oscillazioni della densità fermionica caratteristico di un modo di respiro di un mesone confinato.

Requisiti​

Prima di iniziare questo tutorial, installa quanto segue:

  • Qiskit SDK v2.0 o successivo, con supporto per la visualizzazione

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

  • Pacchetto Pauli Propagation (pip install pauli-prop)

  • NumPy (pip install numpy)

  • Matplotlib (pip install matplotlib)

Configurazione​

Inizia importando le librerie necessarie e definendo le funzioni helper che costruiscono i circuiti quantistici per l'evoluzione temporale LSH. Ci sono tre funzioni principali per la costruzione dei circuiti:

  1. pair_hamiltonian_circuit: Implementa l'unitaria a due qubit UIU_I per l'Hamiltoniana di interazione approssimata tra siti vicini. La decomposizione della porta è: 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: Implementa l'unitaria a due qubit UEU_E per l'energia del campo elettrico approssimata a ciascun sito. La decomposizione della porta è: 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: Assembla il circuito Trotterizzato completo, stratificando i termini di interazione, elettrico e di massa con porte SWAP per gestire la connettività dei 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

Esempio di simulatore su piccola scala​

Innanzitutto, dimostra il flusso di lavoro su piccola scala usando un reticolo a sei siti (12 qubit), in modo da poter verificare la costruzione del circuito e comprendere gli osservabili fisici prima di eseguire su hardware.

Passo 1: Mappa gli input classici in un problema quantistico​

Definisci i parametri fisici corrispondenti al regime di accoppiamento debole studiato nell'articolo (x=100x = 100, m/g=1m/g = 1). I parametri del circuito derivati sono:

  • c=δτ⋅x=0.15c = \delta_\tau \cdot x = 0.15 (parametro di interazione)

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

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

Per ogni conteggio di passi di Trotter, costruisci due circuiti: uno che inizializza un mesone al centro (inverse_mid=True) e uno che prepara il vuoto di accoppiamento forte (inverse_mid=False). Il protocollo di misurazione differenziale sottrae l'evoluzione del vuoto per isolare il segnale dell'adrone.

# 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

Passo 2: Ottimizza il problema per l'esecuzione su hardware quantistico​

Definisci gli osservabili: misurazioni ZZ a singolo qubit su ogni qubit. Da ⟨Z⟩\langle Z \rangle puoi estrarre le probabilità di occupazione e quindi il numero fermionico staggered nf(r)n_f(r) a ciascun sito del reticolo 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

Passo 3: Esegui usando le primitive Qiskit​

Usa StatevectorEstimator per una simulazione esatta senza rumore su piccola scala.

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

Passo 4: Post-elabora e restituisci il risultato nel formato classico desiderato​

Converti i valori di aspettazione nel numero fermionico staggered nf(r,t)n_f(r, t) e applica il protocollo di misurazione differenziale (mesone −- vuoto) per produrre la mappa di calore della propagazione adronica. Questo riproduce la struttura della Figura 3 dall'articolo di riferimento: sito del reticolo rr sull'asse x, passo di Trotter (tempo) tt sull'asse y, e nf(r,t)n_f(r,t) come scala di colore.

# 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

Esempio su hardware su larga scala​

Ora scaliamo a un reticolo a 30 siti (60 qubit) su hardware IBM Quantum. A questa scala, il circuito a 10 passi di Trotter comprende oltre 3400 porte a due qubit e 14.000 porte a singolo qubit.

Passi 1-4 (compressi in un unico blocco di codice)​

Aspetti chiave del flusso di lavoro hardware:

  • 10 passi di Trotter per i circuiti del mesone e del vuoto (interfogliati per una deriva minima)

  • Transpilazione con optimization_level=1 — il layout del circuito è già isomorfo alla topologia del dispositivo (una catena lineare), quindi non sono necessari SWAP di instradamento. Il transpiler viene usato esclusivamente per selezionare una catena a basso rumore di qubit fisici e decomporre le porte nel set di porte nativo.

  • EstimatorV2 con mitigazione degli errori di lettura TREX e Pauli twirling

  • Sessione Batch per inviare tutti i job insieme

# -------------------------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 classico tramite Pauli Propagation​

Il Pauli Propagation Method (PPM) fornisce una simulazione classica senza rumore del circuito quantistico retro-propagando gli osservabili misurati attraverso il circuito nella rappresentazione di Heisenberg. Sotto gli strati di Clifford (porte CNOT, H, S, X), gli operatori di Pauli si mappano su altri operatori di Pauli senza aumentare il numero di termini. Gli strati non-Clifford (le porte RzR_z nel circuito) possono causare ramificazioni — nel caso peggiore, raddoppiando il numero di termini — ma molti rami hanno coefficienti piccoli e possono essere troncati.

Il flusso di lavoro con pauli-prop è:

  1. Dividi il circuito nelle sue parti Clifford e non-Clifford usando evolve_through_cliffords.

  2. Propaga ciascun osservabile attraverso la parte non-Clifford usando propagate_through_circuit, mantenendo fino a max_terms termini di Pauli e scartando i termini con coefficienti inferiori alla soglia di troncamento atol.

  3. Evolvi il risultato attraverso la parte Clifford usando il supporto Clifford integrato di Qiskit.

  4. Estrai il valore di aspettazione sommando i coefficienti dei termini di Pauli diagonali (contenenti solo II e ZZ).

Soglia di troncamento​

Il parametro atol in propagate_through_circuit controlla quanto aggressivamente vengono eliminati i piccoli rami di Pauli. Una soglia molto stretta (per esempio, 1e-12) mantiene quasi tutti i rami e fornisce risultati esatti, ma il tempo di simulazione cresce ripidamente con la profondità del circuito; la simulazione a 120 qubit nell'articolo ha richiesto circa 8,5 ore con le impostazioni predefinite. Aumentare la soglia (per esempio, a 1e-6 o 1e-3) scarta i termini i cui coefficienti scendono sotto quel valore, riducendo drasticamente il numero di termini tracciati e velocizzando il calcolo. Il compromesso è un piccolo errore di approssimazione controllabile che puoi convalidare confrontando i risultati a soglie diverse.

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

Prossimi passi​

Se hai trovato interessante questo lavoro, prendi in considerazione di esplorare il seguente materiale:

Consigli

Riferimenti​

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

Fonte: IBM Quantum docs — aggiornato il 9 set 2026
Versione inglese su doQumentation — aggiornata il 9 set 2026
Problema con il sito o la traduzione?Problema di contenuto? Modifica su IBM Quantum docs