Skip to main content
IBM Quantum Platform

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

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

Acquis d'apprentissage

À l'issue de ce tutoriel, vous aurez acquis les connaissances suivantes :

  • Comment les théories de jauge sur treillis 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 type Trotter pour un hamiltonien approximatif d'une théorie de jauge SU(2) et les transposer en qubits
  • Comment exécuter ces circuits sur du matériel d’ IBM Quantum® s à l’aide de la primitive « Qiskit Estimator » avec atténuation des erreurs de lecture

Prérequis

Nous vous recommandons de vous familiariser avec les sujets suivants :

Arrière-plan

Motivation

La chromodynamique quantique (QCD), théorie de jauge SU(3) de la force forte, lie les quarks pour former des hadrons et régit le confinement et la rupture des cordes. Les méthodes classiques de la QCD sur réseau permettent d'étudier avec brio les propriétés statiques, mais ne permettent pas de simuler la dynamique en temps réel en raison du problème du signe. Les ordinateurs quantiques permettent de contourner cet obstacle en codant directement les degrés de liberté des champs de jauge sur les qubits.

Ce tutoriel présente une simulation de ce type : il s'agit d'utiliser le matériel d' IBM Quantum pour simuler la propagation en temps réel des hadrons dans une théorie de jauge sur réseau SU(2) de dimension (1+1) — la théorie de jauge non abélienne la plus simple et un tremplin vers la QCD complète.

L'hamiltonien de Kogut-Susskind

La théorie est formulée sur un réseau spatial de type « 1D » comportant des fermions (matière) disposés en quinconce sur les sites et des champs de jauge SU(2) sur les liaisons. Après mise à l'échelle pour obtenir une forme adimensionnelle, l'hamiltonien s'écrit :

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

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

Le cadre « Loop-String-Hadron » (LSH)

L'un des principaux défis réside dans le fait que l'espace de Hilbert du champ de jauge sur chaque liaison 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 sous la jauge : des boucles de flux, des cordes reliant des charges séparées et des hadrons (paires de fermions singulets de jauge en un site). Dans la base LSH, la loi de Gauss est automatiquement respectée par construction; chaque état de base est donc 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 de fermions local est défini à partir de ces valeurs comme suit : 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.

De l'hamiltonien complet au circuit quantique : trois approximations clés

Le circuit quantique ne simule pas exactement l'hamiltonien SU(2) complet. Elle 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 bien comprendre ce qui est approximé et ce qui ne l'est pas :

Approximation 1 — Limite de couplage faible pour l’ HIH_I : L’hamiltonien à interactions complètes HI(LSH)H_I^{\text{(LSH)}} (équation La formule (16) [1] contient des préfacteurs qui dépendent du nombre quantique bosonique nln_l via des termes tels que 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 présentant une valeur élevée de 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é. L'hamiltonien d'interaction se réduit alors à un saut entre voisins les plus proches, de nature purement locale :

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épendante de l' nln_l e et n'agit que sur les qubits fermioniques (ni,no)(n_i, n_o).

Approximation 2 — Flux moyen global pour l’ HEH_E : l’énergie électrique dépend de nln_l à chaque maillon. Dans le vide à couplage faible, l' nln_l e est importante et approximativement uniforme. Remplacer les valeurs de l' nln_l, qui dépendent du site, par une seule nˉl\bar{n}_l moyenne globale, ce qui fait de l' HEH_E une phase diagonale proportionnelle à la configuration des fermions à 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'\} correspond à la 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 vous pouvez ignorer.

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

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 tend vers zéro lorsque δτ→0\delta_\tau \to 0. Nous fixons δτ=0.0015\delta_\tau = 0.0015 tout au long du calcul.

Il résulte de ces trois approximations 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é intégré dans les paramètres effectifs. On obtient ainsi un circuit compact comportant 2N2N qubits pour NN sites du réseau, chaque étape de Trotter présentant une profondeur de porte constante de deux qubits (13 par étape).

Ce que simule ce tutoriel

Ce tutoriel simule la propagation des hadrons : à partir d'un vide à couplage fort (un état de produit), placez un méson au centre du réseau et suivez son évolution dans le temps. Le protocole de mesure différentielle — consistant à faire fonctionner le circuit avec et sans le méson central, puis à soustraire les résultats — permet d'isoler le signal hadronique cohérent à la fois du bruit matériel et des effets de frontière. Il en résulte un motif en cône de lumière caractérisé par des oscillations de la densité des fermions, propres à un mode de respiration d'un méson confiné.

Exigences

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

  • Qiskit SDK v2.0 ou version ultérieure, avec prise en charge de la visualisation
  • Qiskit Runtime v0.22 ou version ultérieure (pip install qiskit-ibm-runtime)
  • Bibliothèque de propagation de Pauli (pip install pauli-prop)
  • NumPy (pip install numpy)
  • Matplotlib (pip install matplotlib)

Configuration

Commencez par importer les bibliothèques nécessaires et par définir les fonctions d'aide qui permettent de construire les circuits quantiques pour l'évolution temporelle LSH. Il existe trois fonctions essentielles pour la conception de circuits :

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

# 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

Commencez par illustrer le déroulement du processus à petite échelle à l'aide d'un réseau à six sites (12 qubits), afin de pouvoir vérifier la construction du circuit et comprendre les observables physiques avant de lancer l'exécution sur le matériel.

Étape 1 : Mettre en correspondance les entrées classiques avec un problème quantique

Définissez 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 qui en découlent sont les suivants :

  • 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 pas de Trotter, on construit deux circuits : l'un initialisant un méson au centre (inverse_mid=True) et l'autre préparant le vide à couplage fort (inverse_mid=False). Le protocole de mesure différentielle soustrait l'évolution du vide afin d'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}")

Output:

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)

Output:

Circuit for 1 Trotter step: 12 qubits, depth 26
Output of the previous code cell

Étape 2 : Optimiser le problème en vue de son exécution sur du matériel quantique

Définir les grandeurs observables : mesures d’ ZZ s d’un seul qubit sur chaque qubit. À partir de ⟨Z⟩\langle Z \rangle, vous pouvez extraire les probabilités d'occupation, puis le nombre de fermions échelonné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)}")

Output:

Number of observables: 12

Étape 3 : Exécutez la commande à l'aide d' Qiskit primitives

Utilisez StatevectorEstimator pour une simulation exacte et 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")

Output:

Computed expectation values for 10 Trotter steps

Étape 4 : Traitement ultérieur et restitution du résultat dans le format classique souhaité

Convertir les valeurs attendues en nombre de fermions décalés nf(r,t)n_f(r, t) et appliquer le protocole de mesure différentielle ( −- -méson-vide) afin de générer la carte thermique de propagation des hadrons. Ce graphique reproduit la structure de la figure 3 de l'article de référence : les sites 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 couleurs.

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

Output of the previous code cell

Exemple de matériel à grande échelle

Nous passons désormais à un réseau de 30 sites (60 qubits) sur le matériel d' IBM Quantum. À cette échelle, le circuit de 10 étapes de Trotter comprend plus de 3 400 portes à deux qubits et 14 000 portes à un seul qubit.

Étapes 1 à 4 (regroupées dans un seul bloc de code)

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

  • 10 pas de Trotter pour les circuits mésoniques et de vide (entrelacés pour minimiser la dérive)
  • Transpilation avec optimization_level=1 — la configuration du circuit est déjà isomorphe à la topologie du dispositif (une chaîne linéaire); aucun SWAP de routage n'est donc nécessaire. Le transpileur sert uniquement à sélectionner une chaîne de qubits physiques à faible bruit et à décomposer les portes en un ensemble de portes natives.
  • EstimatorV2 grâce à l'atténuation des erreurs de lecture TREX et à la rotation de Pauli
  • Batch session permettant de soumettre tous les travaux en une seule fois
# -------------------------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:

Output of the previous code cell

Analyse comparative classique par propagation de Pauli

La méthode de propagation de Pauli (PPM) permet une simulation classique sans bruit du circuit quantique en propageant en retour les observables mesurées à travers le circuit dans le cadre de Heisenberg. Dans 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 » du circuit) peuvent entraîner des ramifications — qui, dans le pire des cas, doublent le nombre de termes — mais de nombreuses ramifications ont des coefficients faibles et peuvent être tronquées.

Le déroulement des opérations avec pauli-prop est le suivant :

  1. Divisez le circuit en ses parties « Clifford » et « non-Clifford » à l'aide de evolve_through_cliffords.
  2. atolPropager chaque observable à travers la partie non-Clifford à l'aide de propagate_through_circuit, en conservant jusqu'à max_terms termes de Pauli et en écartant les termes dont les coefficients sont inférieurs au seuil de troncature.
  3. Faites évoluer le résultat via la partie Clifford en utilisant la prise en charge intégrée de Clifford par Qiskit.
  4. On obtient la valeur attendue en additionnant les coefficients des termes de Pauli diagonaux (qui ne contiennent que II et ZZ ).

Seuil de troncature

Le atol paramètre détermine l'intensité avec laquelle les petites branches de propagate_through_circuit Pauli sont éliminées. Un seuil très serré (par exemple, 1e-12) conserve la quasi-totalité des branches et fournit des résultats exacts, mais la durée de simulation augmente fortement avec la profondeur du circuit; la simulation de 120 qubits présentée dans l 'article a pris environ 8.5 heures avec les paramètres par défaut. Le fait de relever le seuil (par exemple, à 1e-6 ou 1e-3) permet d'écarter les termes dont les coefficients sont inférieurs à cette valeur, ce qui réduit considérablement le nombre de termes pris en compte et accélère le calcul. En contrepartie, on obtient une petite erreur d'approximation maîtrisable, que vous pouvez vérifier en comparant les résultats obtenus avec 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}")

Output:

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:

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:

Output of the previous code cell

Etapes suivantes

Si ce travail vous a intéressé, n'hésitez pas à consulter les ressources suivantes :

Recommandations

Références

[1] 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)

Cette page a-t-elle été utile ?
Signaler un bogue, une coquille ou proposer du contenu sur GitHub.