Aller au contenu principal

Observation d'une dynamique hadronique non abélienne robuste et cohérente sur des processeurs quantiques bruités

Utilisation estimée : 6 minutes sur un processeur Heron (ibm_boston ou équivalent) (REMARQUE : il s'agit d'une estimation uniquement. Votre temps d'exécution peut varier.)

Résultats d'apprentissage​

  • Comment les théories de jauge sur réseau non abéliennes (en particulier SU(2)) peuvent être reformulées à l'aide du cadre Loop-String-Hadron (LSH) pour une simulation quantique efficace

  • Comment construire des circuits d'évolution temporelle de Trotter pour un hamiltonien de théorie de jauge SU(2) approximatif et les faire correspondre à des qubits

  • Comment exécuter ces circuits sur du matériel IBM Quantum® à l'aide de la primitive Estimator de Qiskit avec atténuation des erreurs de lecture

Prérequis​

Contexte​

Motivation​

La chromodynamique quantique (QCD), la théorie de jauge SU(3) de la force forte, lie les quarks en hadrons et régit le confinement et la rupture de corde. Les méthodes classiques de QCD sur réseau excellent pour les propriétés statiques mais ne peuvent pas simuler la dynamique en temps réel en raison du problème de signe. Les ordinateurs quantiques offrent un moyen de contourner cette barrière en encodant directement les degrés de liberté des champs de jauge sur des qubits.

Ce tutoriel présente une telle simulation : utiliser du matériel IBM Quantum pour simuler la propagation d'un hadron en temps réel dans une théorie de jauge sur réseau SU(2) à (1+1) dimensions — la théorie de jauge non abélienne la plus simple et un tremplin vers la QCD complète.

Le hamiltonien de Kogut-Susskind​

La théorie est formulée sur un réseau spatial 1D avec des fermions décalés (matière) sur les sites et des champs de jauge SU(2) sur les liens. Après remise à l'échelle sous forme adimensionnelle, le hamiltonien est :

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

où HEH_E est l'énergie du champ chromoélectrique, HMH_M est le terme de masse décalé, HIH_I est le terme d'interaction matière-jauge (saut), μ=2mgx\mu = 2\frac{m}{g}\sqrt{x} encode la masse du fermion, et x=1g2a2x = \frac{1}{g^2 a^2} est l'intensité de l'interaction. La limite continue de la théorie se situe à N→∞N \to \infty et x→∞x \to \infty.

Le cadre Loop-String-Hadron (LSH)​

Un défi clé est que l'espace de Hilbert du champ de jauge sur chaque lien est de dimension infinie. Le cadre Loop-String-Hadron (LSH) résout ce problème en reformulant la théorie en termes de variables invariantes de jauge — des boucles de flux, des cordes reliant des charges séparées, et des hadrons (paires de fermions singulets de jauge sur un site). Dans la base LSH, la loi de Gauss est satisfaite automatiquement par construction, de sorte que chaque état de base est physique. Chaque site du réseau est caractérisé par trois nombres quantiques (nl,ni,no)(n_l, n_i, n_o) représentant le nombre de boucles, la corde entrante et la corde sortante, où ni,no∈{0,1}n_i, n_o \in \{0,1\} sont fermioniques et nl≥0n_l \geq 0 est bosonique. Le nombre fermionique local est défini à partir de ceux-ci comme nf(r)=ni(r)+no(r)n_f(r) = n_i(r) + n_o(r) pour les sites pairs et nf(r)=2−[ni(r)+no(r)]n_f(r) = 2 - [n_i(r) + n_o(r)] pour les sites impairs.

Du hamiltonien complet au circuit quantique : trois approximations clés​

Le circuit quantique ne simule pas exactement le hamiltonien SU(2) complet. Il met plutôt en œuvre une série contrôlée d'approximations valables dans le régime de couplage faible (x≫1x \gg 1). Il est essentiel de comprendre ce qui est approximé et ce qui ne l'est pas :

Approximation 1 — Limite de couplage faible pour HIH_I : Le hamiltonien d'interaction complet HI(LSH)H_I^{\text{(LSH)}} (Éq. 16 dans [1]) contient des préfacteurs qui dépendent du nombre quantique bosonique nln_l via des termes comme 1/nl+11/\sqrt{n_l+1}. Dans le régime de couplage faible (x≫1x \gg 1), la dynamique est dominée par le terme électrique HEH_E, qui favorise les états avec un grand nln_l. Pour nl≫1n_l \gg 1, le rapport nl/(nl+1)→1n_l/(n_l+1) \to 1 et tous ces préfacteurs se simplifient à l'unité. Le hamiltonien d'interaction se réduit alors à un saut purement local entre plus proches voisins :

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

qui est indépendant de nln_l et n'agit que sur les qubits fermioniques (ni,no)(n_i, n_o).

Approximation 2 — Flux moyen global pour HEH_E : L'énergie électrique dépend de nln_l sur chaque lien. Dans le vide de couplage faible, nln_l est grand et approximativement uniforme. On remplace les valeurs de nln_l dépendantes du site par une seule moyenne globale nˉl\bar{n}_l, faisant de HEH_E une phase diagonale proportionnelle à la configuration fermionique de chaque 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)

où {r′}\{r'\} somme sur les sites dans la configuration fermionique (ni=0,no=1)(n_i=0, n_o=1), et hE0h_E^0 est une phase globale que tu peux ignorer.

Approximation 3 — Trotterisation : L'opérateur d'évolution temporelle pour un pas de durée δτ\delta_\tau est décomposé comme :

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

où c=δτxc = \delta_\tau x, m~=δτμ\tilde{m} = \delta_\tau \mu, et θ=−δτ(nˉl/2+3/4)\theta = -\delta_\tau(\bar{n}_l/2 + 3/4). Cette décomposition de Trotter du premier ordre introduit une erreur qui s'annule lorsque δτ→0\delta_\tau \to 0. Nous fixons δτ=0.0015\delta_\tau = 0.0015 tout au long.

Le résultat de ces trois approximations est que seuls les deux qubits fermioniques par site (ni,no)(n_i, n_o) sont dynamiques — le degré de liberté bosonique nln_l a été absorbé dans des paramètres effectifs. Cela donne un circuit compact avec 2N2N qubits pour NN sites du réseau, où chaque pas de Trotter a une profondeur de porte à deux qubits constante (13 par pas).

Ce que ce tutoriel simule​

Le tutoriel simule la propagation d'un hadron : en partant du vide de couplage fort (un état produit), on place un méson au centre du réseau et on fait évoluer dans le temps. Le protocole de mesure différentielle — exécuter le circuit avec et sans le méson central, puis soustraire — isole le signal hadronique cohérent du bruit matériel et des effets de bord. Le résultat est un motif en cône de lumière d'oscillations de densité fermionique caractéristique d'un mode de respiration de méson confiné.

Prérequis techniques​

Avant de commencer ce tutoriel, installe les éléments suivants :

  • Qiskit SDK v2.0 ou une version ultérieure, avec la prise en charge de la visualisation

  • Qiskit Runtime v0.22 ou une version ultérieure (pip install qiskit-ibm-runtime)

  • Le package Pauli Propagation (pip install pauli-prop)

  • NumPy (pip install numpy)

  • Matplotlib (pip install matplotlib)

Configuration​

Commence par importer les bibliothèques nécessaires et définir les fonctions d'aide qui construisent les circuits quantiques pour l'évolution temporelle LSH. Il existe trois fonctions principales de construction de circuits :

  1. pair_hamiltonian_circuit : Met en œuvre l'unitaire à deux qubits UIU_I pour le hamiltonien d'interaction approximatif entre sites voisins. La décomposition en portes est : 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 : Met en œuvre l'unitaire à deux qubits UEU_E pour l'énergie du champ électrique approximatif à chaque site. La décomposition en portes est : 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 : Assemble le circuit de Trotter complet, en superposant les termes d'interaction, électrique et de masse avec des portes SWAP pour gérer la connectivité des qubits.

# 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

Exemple de simulateur à petite échelle​

Tout d'abord, montre le flux de travail à petite échelle en utilisant un réseau à six sites (12 qubits), afin de pouvoir vérifier la construction du circuit et comprendre les observables physiques avant d'exécuter sur du matériel.

Étape 1 : Faire correspondre les entrées classiques à un problème quantique​

Définis les paramètres physiques correspondant au régime de couplage faible étudié dans l'article (x=100x = 100, m/g=1m/g = 1). Les paramètres de circuit dérivés sont :

  • c=δτ⋅x=0.15c = \delta_\tau \cdot x = 0.15 (paramètre d'interaction)

  • θ=−δτ(nˉl/2+3/4)=0.01\theta = -\delta_\tau (\bar{n}_l/2 + 3/4) = 0.01 (phase du champ électrique)

  • m~=δτ⋅μ=0.03\tilde{m} = \delta_\tau \cdot \mu = 0.03 (paramètre de masse)

Pour chaque nombre de pas de Trotter, construis deux circuits : l'un initialisant un méson au centre (inverse_mid=True) et l'autre préparant le vide de couplage fort (inverse_mid=False). Le protocole de mesure différentielle soustrait l'évolution du vide pour isoler le signal hadronique.

# 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

Étape 2 : Optimiser le problème pour l'exécution sur du matériel quantique​

Définis les observables : des mesures ZZ à un seul qubit sur chaque qubit. À partir de ⟨Z⟩\langle Z \rangle, tu peux extraire les probabilités d'occupation puis le nombre de fermions décalés nf(r)n_f(r) à chaque site du réseau 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

Étape 3 : Exécuter en utilisant les primitives Qiskit​

Utilise StatevectorEstimator pour une simulation exacte sans bruit à petite échelle.

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

Étape 4 : Post-traiter et retourner le résultat dans le format classique souhaité​

Convertis les valeurs moyennes en nombre de fermions décalés nf(r,t)n_f(r, t) et applique le protocole de mesure différentielle (méson −- vide) pour produire la carte thermique de propagation du hadron. Cela reproduit la structure de la Figure 3 de l'article de référence : le site du réseau rr sur l'axe des x, le pas de Trotter (temps) tt sur l'axe des y, et nf(r,t)n_f(r,t) comme échelle de couleur.

# 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

Exemple à grande échelle sur du matériel​

Nous passons maintenant à un réseau de 30 sites (60 qubits) sur du matériel IBM Quantum. À cette échelle, le circuit à 10 pas de Trotter comprend plus de 3400 portes à deux qubits et 14 000 portes à un qubit.

Étapes 1 à 4 (compressées en un seul bloc de code)​

Aspects clés du flux de travail matériel :

  • 10 pas de Trotter pour les circuits du méson et du vide (entrelacés pour une dérive minimale)

  • Transpilation avec optimization_level=1 — la disposition du circuit est déjà isomorphe à la topologie de l'appareil (une chaîne linéaire), donc aucun routage par SWAP n'est nécessaire. Le transpileur est utilisé uniquement pour sélectionner une chaîne de qubits physiques à faible bruit et décomposer les portes dans l'ensemble de portes natif.

  • EstimatorV2 avec atténuation des erreurs de lecture TREX et twirling de Pauli

  • Session Batch pour soumettre tous les jobs ensemble

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

Analyse comparative classique via Pauli Propagation​

La méthode de propagation de Pauli (PPM) fournit une simulation classique sans bruit du circuit quantique en rétropropageant les observables mesurées à travers le circuit dans l'image de Heisenberg. Sous les couches de Clifford (portes CNOT, H, S, X), les opérateurs de Pauli se transforment en d'autres opérateurs de Pauli sans augmenter le nombre de termes. Les couches non-Clifford (les portes RzR_z dans le circuit) peuvent provoquer une ramification — dans le pire des cas, doublant le nombre de termes — mais de nombreuses branches ont de petits coefficients et peuvent être tronquées.

Le flux de travail avec pauli-prop est le suivant :

  1. Diviser le circuit en ses parties Clifford et non-Clifford à l'aide de evolve_through_cliffords.

  2. Propager chaque observable à travers la partie non-Clifford à l'aide de propagate_through_circuit, en conservant jusqu'à max_terms termes de Pauli et en supprimant les termes dont les coefficients sont inférieurs au seuil de troncature atol.

  3. Faire évoluer le résultat à travers la partie Clifford en utilisant la prise en charge intégrée de Clifford de Qiskit.

  4. Extraire la valeur moyenne en additionnant les coefficients des termes de Pauli diagonaux (contenant uniquement II et ZZ).

Seuil de troncature​

Le paramètre atol dans propagate_through_circuit contrôle avec quelle agressivité les petites branches de Pauli sont élaguées. Un seuil très strict (par exemple, 1e-12) conserve presque toutes les branches et donne des résultats exacts, mais le temps de simulation augmente fortement avec la profondeur du circuit ; la simulation à 120 qubits de l'article a pris environ 8,5 heures avec les paramètres par défaut. Augmenter le seuil (par exemple, à 1e-6 ou 1e-3) élimine les termes dont les coefficients sont inférieurs à cette valeur, réduisant considérablement le nombre de termes suivis et accélérant le calcul. Le compromis est une erreur d'approximation faible et contrôlable que tu peux valider en comparant les résultats à différents seuils.

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

Étapes suivantes​

Si tu as trouvé ce travail intéressant, envisage d'explorer les documents suivants :

Recommandations

Références​

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

Source : IBM Quantum docs — mis à jour le 9 sept. 2026
Version anglaise sur doQumentation — mise à jour le 9 sept. 2026
Problème avec le site ou la traduction ?Problème de contenu ? Modifier sur IBM Quantum docs