{
"cells": [
{
"cell_type": "markdown",
"id": "title-cell",
"metadata": {},
"source": [
"---\n",
"title: \"Observation d'une dynamique hadronique non abélienne robuste et cohérente sur des processeurs quantiques sujets au bruit\"\n",
"description: \"Simuler la dynamique des hadrons dans la théorie de jauge sur réseau SU(2) à l'aide du cadre LSH sur du matériel quantique d' IBM.\"\n",
"---\n",
"\n",
"{/* cspell:ignore Kogut Susskind expvals Pstep Nstep vmax vmin imshow fontsize cbar Ilčić */}\n",
"\n",
"\n",
"\n",
"# Observation d'une dynamique hadronique non abélienne robuste et cohérente sur des processeurs quantiques sujets au bruit\n",
"\n",
"*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.)*\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "learning-outcomes",
"metadata": {},
"source": [
"\n",
"\n",
"## Acquis d'apprentissage\n",
"\n",
"À l'issue de ce tutoriel, vous aurez acquis les connaissances suivantes :\n",
"\n",
"* 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\n",
"* 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\n",
"* 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\n",
"\n",
"\n",
"\n",
"## Prérequis\n",
"\n",
"Nous vous recommandons de vous familiariser avec les sujets suivants :\n",
"\n",
"* [Notions fondamentales sur les circuits et les portes quantiques](/learning/courses/basics-of-quantum-information)\n",
"* [Introduction à la primitive « Estimator » de Qiskit](/docs/guides/get-started-with-estimator)\n",
"* Une connaissance de base des concepts de la théorie quantique des champs (utile mais non obligatoire; la section « Contexte » aborde les éléments essentiels)\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "background",
"metadata": {},
"source": [
"\n",
"\n",
"## Arrière-plan\n",
"\n",
"\n",
"\n",
"### Motivation\n",
"\n",
"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.\n",
"\n",
"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.\n",
"\n",
"\n",
"\n",
"### L'hamiltonien de Kogut-Susskind\n",
"\n",
"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 :\n",
"\n",
"$W = H_E^{\\text{(KS)}} + \\mu H_M + x H_I^{\\text{(KS)}},$\n",
"\n",
"où $H_E$ représente l'énergie du champ chromoélectrique, $H_M$ le terme de masse décalée, $H_I$ le terme d'interaction matière-jauge (saut), $\\mu = 2\\frac{m}{g}\\sqrt{x}$ la masse du fermion, et $x = \\frac{1}{g^2 a^2}$ l'intensité d'interaction. La limite continue de la théorie est décrite sur $N \\to \\infty$ et $x \\to \\infty$.\n",
"\n",
"\n",
"\n",
"### Le cadre « Loop-String-Hadron » (LSH)\n",
"\n",
"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 $(n_l, n_i, n_o)$ représentant le nombre de boucles, la corde entrante et la corde sortante, où $n_i, n_o \\in \\{0,1\\}$ sont fermioniques et $n_l \\geq 0$ est bosonique. Le nombre de fermions local est défini à partir de ces valeurs comme suit : $n_f(r) = n_i(r) + n_o(r)$ pour les sites pairs et $n_f(r) = 2 - [n_i(r) + n_o(r)]$ pour les sites impairs.\n",
"\n",
"\n",
"\n",
"### De l'hamiltonien complet au circuit quantique : trois approximations clés\n",
"\n",
"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 \\gg 1$ ). Il est essentiel de bien comprendre ce qui est approximé et ce qui ne l'est pas :\n",
"\n",
"**Approximation 1 — Limite de couplage faible pour l’ $H_I$ :** L’hamiltonien à interactions complètes $H_I^{\\text{(LSH)}}$ (équation La formule (16) [\\[1\\]](#references) contient des préfacteurs qui dépendent du nombre quantique bosonique $n_l$ via des termes tels que $1/\\sqrt{n_l+1}$. Dans le régime de couplage faible ( $x \\gg 1$ ), la dynamique est dominée par le terme électrique $H_E$, qui favorise les états présentant une valeur élevée de $n_l$. Pour $n_l \\gg 1$, le rapport $n_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 :\n",
"\n",
"$H_I^{\\text{approx}} = -\\sum_r \\left[\\sigma^-(r)\\sigma^+(r+1) + \\sigma^+(r)\\sigma^-(r+1)\\right],$\n",
"\n",
"qui est indépendante de l' $n_l$ e et n'agit que sur les qubits fermioniques $(n_i, n_o)$.\n",
"\n",
"**Approximation 2 — Flux moyen global pour l’ $H_E$ :** l’énergie électrique dépend de $n_l$ à chaque maillon. Dans le vide à couplage faible, l' $n_l$ e est importante et approximativement uniforme. Remplacer les valeurs de l' $n_l$, qui dépendent du site, par une seule $\\bar{n}_l$ moyenne globale, ce qui fait de l' $H_E$ une phase diagonale proportionnelle à la configuration des fermions à chaque site :\n",
"\n",
"$H_E^{\\text{approx}} = N h_E^0 + \\sum_{\\{r'\\}} \\left(\\frac{\\bar{n}_l}{2} + \\frac{3}{4}\\right)$\n",
"\n",
"où $\\{r'\\}$ correspond à la somme sur les sites dans la configuration fermionique $(n_i=0, n_o=1)$, et $h_E^0$ est une phase globale que vous pouvez ignorer.\n",
"\n",
"**Approximation 3 — Trotterisation :** L'opérateur d'évolution temporelle pour un pas de durée $\\delta_\\tau$ se décompose comme suit :\n",
"\n",
"$e^{-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}}}$\n",
"\n",
"où $c = \\delta_\\tau x$, $\\tilde{m} = \\delta_\\tau \\mu$ et $\\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 $\\delta_\\tau \\to 0$. Nous fixons $\\delta_\\tau = 0.0015$ tout au long du calcul.\n",
"\n",
"**Il résulte** de ces trois approximations que seuls les deux qubits fermioniques par site $(n_i, n_o)$ sont dynamiques — le degré de liberté bosonique $n_l$ a été intégré dans les paramètres effectifs. On obtient ainsi un circuit compact comportant $2N$ qubits pour $N$ sites du réseau, chaque étape de Trotter présentant une profondeur de porte constante de deux qubits (13 par étape).\n",
"\n",
"\n",
"\n",
"### Ce que simule ce tutoriel\n",
"\n",
"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é.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "requirements",
"metadata": {},
"source": [
"\n",
"\n",
"## Exigences\n",
"\n",
"Avant de commencer ce tutoriel, veuillez installer les éléments suivants :\n",
"\n",
"* Qiskit SDK v2.0 ou version ultérieure, avec prise en charge [de la visualisation](/docs/api/qiskit/visualization)\n",
"* Qiskit Runtime v0.22 ou version ultérieure (`pip install qiskit-ibm-runtime`)\n",
"* Bibliothèque de propagation de Pauli (`pip install pauli-prop`)\n",
"* NumPy (`pip install numpy`)\n",
"* Matplotlib (`pip install matplotlib`)\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "setup-header",
"metadata": {},
"source": [
"\n",
"\n",
"## Configuration\n",
"\n",
"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 :\n",
"\n",
"1. **`pair_hamiltonian_circuit`**: Met en œuvre l' $U_I$ unitaire à deux qubits pour l'hamiltonien d'interaction approximatif entre des sites voisins. La décomposition en portes est la suivante : $\\text{CNOT} \\to H \\to R_z(-c) \\to \\text{CNOT} \\to R_z(c) \\to \\text{CNOT} \\to H \\to \\text{CNOT}$.\n",
"\n",
"2. **`electric_hamiltonian_circuit`**: Met en œuvre l' $U_E$ unitaire à deux qubits pour l'énergie approximative du champ électrique à chaque site. La décomposition en portes est la suivante : $X \\to R_z(\\theta/2) \\to \\text{CNOT} \\to R_z(-\\theta/2) \\to \\text{CNOT} \\to R_z(\\theta/2) \\to X$.\n",
"\n",
"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.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "setup-imports",
"metadata": {},
"outputs": [],
"source": [
"# Import libraries\n",
"\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from matplotlib.colors import TwoSlopeNorm\n",
"from qiskit.circuit import QuantumCircuit\n",
"from qiskit.quantum_info import SparsePauliOp\n",
"from typing import Optional\n",
"\n",
"import warnings\n",
"\n",
"warnings.filterwarnings(\"ignore\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "setup-functions",
"metadata": {},
"outputs": [],
"source": [
"def pair_hamiltonian_circuit(c: float) -> QuantumCircuit:\n",
" \"\"\"Two-qubit unitary for the approximate interaction Hamiltonian H_I.\n",
"\n",
" Implements exp(-i * c * H_I^approx) for one pair of neighboring sites,\n",
" where c = delta_tau * x.\n",
" \"\"\"\n",
" qc_temp = QuantumCircuit(2)\n",
" qc_temp.cx(1, 0)\n",
" qc_temp.h(1)\n",
" qc_temp.rz(-c, 1)\n",
" qc_temp.cx(0, 1)\n",
" qc_temp.rz(c, 1)\n",
" qc_temp.cx(0, 1)\n",
" qc_temp.h(1)\n",
" qc_temp.cx(1, 0)\n",
" return qc_temp\n",
"\n",
"\n",
"def electric_hamiltonian_circuit(theta: float) -> QuantumCircuit:\n",
" \"\"\"Two-qubit unitary for the approximate electric field Hamiltonian H_E.\n",
"\n",
" Implements exp(-i * theta * H_E^approx) for one lattice site,\n",
" where theta = -delta_tau * (n_bar_l / 2 + 3/4).\n",
" \"\"\"\n",
" qc_temp = QuantumCircuit(2)\n",
" qc_temp.x(0)\n",
" qc_temp.rz(theta / 2, 0)\n",
" qc_temp.cx(0, 1)\n",
" qc_temp.rz(-theta / 2, 1)\n",
" qc_temp.cx(0, 1)\n",
" qc_temp.rz(theta / 2, 1)\n",
" qc_temp.x(0)\n",
" return qc_temp\n",
"\n",
"\n",
"def construct_circuit(\n",
" num_lattice_point: int,\n",
" num_trotter_steps: int,\n",
" c: float,\n",
" theta: float,\n",
" m: float,\n",
" theory: Optional[int] = 2,\n",
" barriers: Optional[bool] = False,\n",
" measurement: Optional[bool] = False,\n",
" add_init_state: Optional[bool] = True,\n",
" inverse_mid: Optional[bool] = False,\n",
") -> QuantumCircuit:\n",
" \"\"\"Construct the full Trotterized time-evolution circuit.\n",
"\n",
" Builds a circuit implementing n Trotter steps of the approximate SU(2)\n",
" LSH Hamiltonian evolution. The qubit layout uses a zigzag ordering:\n",
" n_i(0), n_i(1), n_o(0), n_o(1), n_i(2), n_i(3), n_o(2), n_o(3), ...\n",
" which minimizes the number of SWAP layers needed.\n",
"\n",
" Args:\n",
" num_lattice_point: Number of lattice sites\n",
" (num_qubits = 2 * num_lattice_point).\n",
" num_trotter_steps: Number of Trotter steps.\n",
" c: Interaction parameter (delta_tau * x).\n",
" theta: Electric field phase parameter.\n",
" m: Mass parameter (m_tilde = delta_tau * mu).\n",
" theory: 1 for single chain, 2 for SU(2). Default 2.\n",
" barriers: Insert barriers between Trotter layers for\n",
" visualization.\n",
" measurement: Append measurements at the end.\n",
" add_init_state: Prepare the half-filled (strong-coupling vacuum)\n",
" initial state.\n",
" inverse_mid: Swap the central sites\n",
" (for differential measurement protocol).\n",
" \"\"\"\n",
" num_qubits = theory * num_lattice_point\n",
" qc = QuantumCircuit(num_qubits)\n",
"\n",
" if num_trotter_steps <= 0:\n",
" return qc\n",
"\n",
" # --- Initial state preparation ---\n",
" if add_init_state:\n",
" i = 1\n",
" while i < num_lattice_point:\n",
" for j in range(theory):\n",
" qc.x(i + j * num_lattice_point)\n",
" i = i + 2\n",
" if inverse_mid:\n",
" mid_lattice_qubits = [num_qubits // 2 - 1, num_qubits // 2]\n",
" qc.x(mid_lattice_qubits)\n",
" else:\n",
" i = 1\n",
" while i < num_qubits - 1:\n",
" qc.swap(i, i + 1)\n",
" i = i + 4\n",
"\n",
" # --- Trotter steps ---\n",
" for step in range(num_trotter_steps):\n",
" if barriers:\n",
" qc.barrier()\n",
"\n",
" # First SWAP layer (skipped at step 0 — absorbed into initial state mapping)\n",
" if step > 0:\n",
" i = 1\n",
" while i < num_qubits - 1:\n",
" qc.swap(i, i + 1)\n",
" i = i + 4\n",
"\n",
" # First layer of pair interactions\n",
" j = 0\n",
" while j < num_qubits - 2:\n",
" circ = pair_hamiltonian_circuit(c)\n",
" qc.compose(circ, [j, j + 1], inplace=True)\n",
" j = j + 2\n",
" if num_lattice_point % 2 == 0:\n",
" circ = pair_hamiltonian_circuit(c)\n",
" qc.compose(circ, [j, j + 1], inplace=True)\n",
"\n",
" # Second SWAP layer\n",
" i = 1\n",
" while i < num_qubits - 1:\n",
" qc.swap(i, i + 1)\n",
" i = i + theory\n",
"\n",
" # Second layer of pair interactions\n",
" j = 2\n",
" while j < num_qubits - 3:\n",
" circ = pair_hamiltonian_circuit(c)\n",
" qc.compose(circ, [j, j + 1], inplace=True)\n",
" j = j + 2\n",
" if num_lattice_point % 2 != 0:\n",
" circ = pair_hamiltonian_circuit(c)\n",
" qc.compose(circ, [j, j + 1], inplace=True)\n",
"\n",
" # Third SWAP layer\n",
" i = 3\n",
" while i < num_qubits - 1:\n",
" qc.swap(i, i + 1)\n",
" i = i + 2 * theory\n",
"\n",
" # Electric field term\n",
" if theta != 0:\n",
" e_circ = electric_hamiltonian_circuit(theta)\n",
" for j in range(num_lattice_point):\n",
" qc.compose(e_circ, [2 * j, 2 * j + 1], inplace=True)\n",
"\n",
" # Mass term: Rz(-m_tilde) for even sites, Rz(m_tilde) for odd sites\n",
" for q in range(num_qubits):\n",
" if q % 2 == 0:\n",
" qc.rz(-1 * m, q)\n",
" else:\n",
" qc.rz(m, q)\n",
"\n",
" if measurement:\n",
" qc.measure_all()\n",
"\n",
" return qc"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "setup-postprocess",
"metadata": {},
"outputs": [],
"source": [
"def get_probabilities(expval: float):\n",
" \"\"\"Convert a Z-expectation value to site occupation probability.\n",
"\n",
" Since = p(0) - p(1), the occupation probability is p(1) = (1 - ) / 2.\n",
" \"\"\"\n",
" p1 = round((1 - expval) / 2, 3)\n",
" return p1\n",
"\n",
"\n",
"def get_number(expval_data, num_lattice_point):\n",
" \"\"\"Convert raw Z-expectation values to staggered fermion number n_f at each site.\n",
"\n",
" n_f(r) = n_i(r) + n_o(r) for even r\n",
" n_f(r) = 2 - [n_i(r) + n_o(r)] for odd r\n",
"\n",
" The two qubits per site encode (n_i, n_o), and occupation probabilities\n",
" give us and .\n",
" \"\"\"\n",
" N = []\n",
" for expvals in expval_data:\n",
" Pstep = [get_probabilities(expval) for expval in expvals]\n",
" Nstep = []\n",
" for k in range(num_lattice_point):\n",
" val = Pstep[2 * k] + Pstep[2 * k + 1]\n",
" a = 2 * (k % 2) + (1 - 2 * (k % 2)) * val\n",
" Nstep.append(float(a))\n",
" N.append(Nstep)\n",
" return N\n",
"\n",
"\n",
"def calculate_difference(N, N_mid, num_lattice_point):\n",
" \"\"\"Differential measurement protocol: |n_f(meson) - n_f(vacuum)|.\n",
"\n",
" Subtracting the vacuum (SCV) evolution from the meson evolution\n",
" isolates the coherent hadron signal from symmetric noise and boundary effects.\n",
" \"\"\"\n",
" N_diff = []\n",
" for i in range(len(N)):\n",
" Nstep_diff = []\n",
" for j in range(num_lattice_point):\n",
" Nstep_diff.append(abs(N[i][j] - N_mid[i][j]))\n",
" N_diff.append(Nstep_diff)\n",
" return N_diff"
]
},
{
"cell_type": "markdown",
"id": "sim-header",
"metadata": {},
"source": [
"\n",
"\n",
"## Exemple de simulateur à petite échelle\n",
"\n",
"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.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "step1-header",
"metadata": {},
"source": [
"\n",
"\n",
"### Étape 1 : Mettre en correspondance les entrées classiques avec un problème quantique\n",
"\n",
"Définissez les paramètres physiques correspondant au régime de couplage faible étudié dans l'article ( $x = 100$, $m/g = 1$ ). Les paramètres de circuit qui en découlent sont les suivants :\n",
"\n",
"* $c = \\delta_\\tau \\cdot x = 0.15$ (paramètre d'interaction)\n",
"* $\\theta = -\\delta_\\tau (\\bar{n}_l/2 + 3/4) = 0.01$ (phase du champ électrique)\n",
"* $\\tilde{m} = \\delta_\\tau \\cdot \\mu = 0.03$ (paramètre de masse)\n",
"\n",
"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.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "step1-params",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Lattice sites: 6, Qubits: 12\n",
"Parameters: c=0.15, theta=0.01, m_tilde=0.03\n"
]
}
],
"source": [
"# Physical / circuit parameters\n",
"num_lattice_point = 6 # 6 lattice sites -> 12 qubits for SU(2)\n",
"num_qubits = 2 * num_lattice_point\n",
"c = 0.15 # delta_tau * x\n",
"theta = 0.01 # electric field phase\n",
"m = 0.03 # m_tilde = delta_tau * mu\n",
"trotter_steps = range(1, 11) # 10 Trotter steps\n",
"\n",
"print(f\"Lattice sites: {num_lattice_point}, Qubits: {num_qubits}\")\n",
"print(f\"Parameters: c={c}, theta={theta}, m_tilde={m}\")"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "step1-circuits",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Circuit for 1 Trotter step: 12 qubits, depth 26\n"
]
},
{
"data": {
"text/plain": [
""
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Build circuits: meson initial state and vacuum (SCV) initial state\n",
"circuits_mid = [\n",
" construct_circuit(\n",
" num_lattice_point,\n",
" d,\n",
" c,\n",
" theta,\n",
" m,\n",
" barriers=False,\n",
" measurement=False,\n",
" add_init_state=True,\n",
" inverse_mid=True,\n",
" )\n",
" for d in trotter_steps\n",
"]\n",
"\n",
"circuits = [\n",
" construct_circuit(\n",
" num_lattice_point,\n",
" d,\n",
" c,\n",
" theta,\n",
" m,\n",
" barriers=False,\n",
" measurement=False,\n",
" add_init_state=True,\n",
" inverse_mid=False,\n",
" )\n",
" for d in trotter_steps\n",
"]\n",
"\n",
"# Visualize a single Trotter step\n",
"print(\n",
" f\"Circuit for 1 Trotter step: {circuits[0].num_qubits} qubits, depth {circuits[0].depth()}\"\n",
")\n",
"circuits[0].draw(\"mpl\", fold=-1)"
]
},
{
"cell_type": "markdown",
"id": "step2-header",
"metadata": {},
"source": [
"\n",
"\n",
"### Étape 2 : Optimiser le problème en vue de son exécution sur du matériel quantique\n",
"\n",
"Définir les grandeurs observables : mesures d’ $Z$ s d’un seul qubit sur chaque qubit. À partir de $\\langle Z \\rangle$, vous pouvez extraire les probabilités d'occupation, puis le nombre de fermions échelonnés $n_f(r)$ à chaque site du réseau $r$.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "step2-observables",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Number of observables: 12\n"
]
}
],
"source": [
"# Z observable on each qubit\n",
"observables = [\n",
" SparsePauliOp(\"I\" * i + \"Z\" + \"I\" * (num_qubits - i - 1))\n",
" for i in range(num_qubits)\n",
"]\n",
"\n",
"print(f\"Number of observables: {len(observables)}\")"
]
},
{
"cell_type": "markdown",
"id": "step3-header",
"metadata": {},
"source": [
"\n",
"\n",
"### Étape 3 : Exécutez la commande à l'aide d' Qiskit primitives\n",
"\n",
"Utilisez `StatevectorEstimator` pour une simulation exacte et sans bruit à petite échelle.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "step3-simulate",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Computed expectation values for 10 Trotter steps\n"
]
}
],
"source": [
"from qiskit.primitives import StatevectorEstimator\n",
"\n",
"estimator = StatevectorEstimator()\n",
"\n",
"# Run meson circuits\n",
"pubs_mid = [(circuit, observables) for circuit in circuits_mid]\n",
"result_mid = estimator.run(pubs_mid).result()\n",
"\n",
"# Run vacuum (SCV) circuits\n",
"pubs = [(circuit, observables) for circuit in circuits]\n",
"result = estimator.run(pubs).result()\n",
"\n",
"# Extract expectation values\n",
"raw_expvals_mid = [\n",
" result_mid[i].data.evs[::-1] for i in range(len(circuits_mid))\n",
"]\n",
"raw_expvals = [result[i].data.evs[::-1] for i in range(len(circuits))]\n",
"\n",
"print(f\"Computed expectation values for {len(raw_expvals)} Trotter steps\")"
]
},
{
"cell_type": "markdown",
"id": "step4-header",
"metadata": {},
"source": [
"\n",
"\n",
"### Étape 4 : Traitement ultérieur et restitution du résultat dans le format classique souhaité\n",
"\n",
"Convertir les valeurs attendues en nombre de fermions décalés $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 $r$ sur l'axe des x, le pas de Trotter (temps) $t$ sur l'axe des y, et $n_f(r,t)$ comme échelle de couleurs.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "step4-postprocess",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Compute fermion numbers\n",
"N_mid_sim = get_number(raw_expvals_mid, num_lattice_point)\n",
"N_sim = get_number(raw_expvals, num_lattice_point)\n",
"N_diff_sim = calculate_difference(N_mid_sim, N_sim, num_lattice_point)\n",
"\n",
"# --- Reproduce Figure 3 style: Staggered Fermionic Occupation Number Dynamics ---\n",
"fig, axes = plt.subplots(1, 3, figsize=(18, 5))\n",
"\n",
"# Convert to numpy arrays for plotting\n",
"N_mid_arr = np.array(N_mid_sim)\n",
"N_arr = np.array(N_sim)\n",
"N_diff_arr = np.array(N_diff_sim)\n",
"\n",
"# Color scheme\n",
"vmax = max(max(sublist) for sublist in N_arr)\n",
"vmin = -vmax\n",
"\n",
"# Meson evolution\n",
"norm1 = TwoSlopeNorm(vmin=vmin, vcenter=0, vmax=vmax)\n",
"im0 = axes[0].imshow(\n",
" N_mid_arr,\n",
" aspect=\"auto\",\n",
" origin=\"lower\",\n",
" cmap=\"RdBu_r\",\n",
" norm=norm1,\n",
" extent=[0, num_lattice_point, 0.5, len(trotter_steps) + 0.5],\n",
")\n",
"axes[0].set_xlabel(\"Lattice site $r$\", fontsize=12)\n",
"axes[0].set_ylabel(\"Trotter step $t$\", fontsize=12)\n",
"axes[0].set_title(\"$n_f(r,t)$ — Meson initial state\", fontsize=12)\n",
"plt.colorbar(im0, ax=axes[0], label=\"$n_f(r,t)$\")\n",
"\n",
"# Vacuum (SCV) evolution\n",
"im1 = axes[1].imshow(\n",
" N_arr,\n",
" aspect=\"auto\",\n",
" origin=\"lower\",\n",
" cmap=\"RdBu_r\",\n",
" norm=norm1,\n",
" extent=[0, num_lattice_point, 0.5, len(trotter_steps) + 0.5],\n",
")\n",
"axes[1].set_xlabel(\"Lattice site $r$\", fontsize=12)\n",
"axes[1].set_ylabel(\"Trotter step $t$\", fontsize=12)\n",
"axes[1].set_title(\"$n_f(r,t)$ — Vacuum (SCV)\", fontsize=12)\n",
"plt.colorbar(im1, ax=axes[1], label=\"$n_f(r,t)$\")\n",
"\n",
"# Differential: meson - vacuum\n",
"norm2 = TwoSlopeNorm(vmin=vmin, vcenter=0, vmax=vmax)\n",
"im2 = axes[2].imshow(\n",
" N_diff_arr,\n",
" aspect=\"auto\",\n",
" origin=\"lower\",\n",
" cmap=\"RdBu_r\",\n",
" norm=norm2,\n",
" extent=[0, num_lattice_point, 0.5, len(trotter_steps) + 0.5],\n",
")\n",
"axes[2].set_xlabel(\"Lattice site $r$\", fontsize=12)\n",
"axes[2].set_ylabel(\"Trotter step $t$\", fontsize=12)\n",
"axes[2].set_title(\n",
" \"Staggered Fermionic Occupation Number Dynamics\\n$|n_f^{\\\\mathrm{meson}} - n_f^{\\\\mathrm{vacuum}}|$\",\n",
" fontsize=12,\n",
")\n",
"plt.colorbar(im2, ax=axes[2], label=\"$n_f(r,t)$\")\n",
"\n",
"plt.suptitle(\n",
" f\"Hadron propagation — {num_lattice_point}-site lattice (StatevectorEstimator)\",\n",
" fontsize=14,\n",
" y=1.02,\n",
")\n",
"plt.tight_layout()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "hardware-header",
"metadata": {},
"source": [
"\n",
"\n",
"## Exemple de matériel à grande échelle\n",
"\n",
"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.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "hardware-steps",
"metadata": {},
"source": [
"\n",
"\n",
"### Étapes 1 à 4 (regroupées dans un seul bloc de code)\n",
"\n",
"Aspects clés du flux de travail matériel :\n",
"\n",
"* 10 pas de Trotter pour les circuits mésoniques et de vide (entrelacés pour minimiser la dérive)\n",
"* 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.\n",
"* `EstimatorV2` grâce à l'atténuation des erreurs de lecture TREX et à la rotation de Pauli\n",
"* `Batch` session permettant de soumettre tous les travaux en une seule fois\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "hardware-code",
"metadata": {},
"outputs": [],
"source": [
"# -------------------------Step 1: Define parameters & build circuits-------------------------\n",
"\n",
"from qiskit_ibm_runtime import QiskitRuntimeService\n",
"from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager\n",
"from qiskit_ibm_runtime import EstimatorV2, Batch\n",
"from qiskit_ibm_runtime.options import (\n",
" EstimatorOptions,\n",
" ResilienceOptionsV2,\n",
" TwirlingOptions,\n",
" DynamicalDecouplingOptions,\n",
")\n",
"\n",
"service = QiskitRuntimeService()\n",
"\n",
"num_lattice_point_hw = 30\n",
"num_qubits_hw = 2 * num_lattice_point_hw # 60 qubits\n",
"c_hw = 0.15\n",
"theta_hw = 0.01\n",
"m_hw = 0.03\n",
"trotter_steps_hw = range(1, 11) # 10 Trotter steps\n",
"\n",
"# Build meson and vacuum circuits\n",
"circuits_mid_hw = [\n",
" construct_circuit(\n",
" num_lattice_point_hw,\n",
" d,\n",
" c_hw,\n",
" theta_hw,\n",
" m_hw,\n",
" barriers=False,\n",
" measurement=False,\n",
" add_init_state=True,\n",
" inverse_mid=True,\n",
" )\n",
" for d in trotter_steps_hw\n",
"]\n",
"\n",
"circuits_hw = [\n",
" construct_circuit(\n",
" num_lattice_point_hw,\n",
" d,\n",
" c_hw,\n",
" theta_hw,\n",
" m_hw,\n",
" barriers=False,\n",
" measurement=False,\n",
" add_init_state=True,\n",
" inverse_mid=False,\n",
" )\n",
" for d in trotter_steps_hw\n",
"]\n",
"\n",
"print(f\"Built {len(circuits_hw)} circuit pairs for {num_qubits_hw} qubits\")\n",
"\n",
"# -------------------------Step 2: Transpile for hardware-------------------------\n",
"# The circuit topology is a linear chain, isomorphic to the device topology.\n",
"# We use optimization_level=1 since no routing SWAPs are needed — the transpiler\n",
"# only needs to select a low-noise qubit chain and decompose to native gates.\n",
"\n",
"backend = service.backend(\"ibm_boston\")\n",
"\n",
"layout = [\n",
" 140,\n",
" 141,\n",
" 142,\n",
" 143,\n",
" 136,\n",
" 123,\n",
" 122,\n",
" 121,\n",
" 116,\n",
" 101,\n",
" 102,\n",
" 103,\n",
" 96,\n",
" 83,\n",
" 82,\n",
" 81,\n",
" 76,\n",
" 61,\n",
" 62,\n",
" 63,\n",
" 64,\n",
" 65,\n",
" 66,\n",
" 67,\n",
" 68,\n",
" 69,\n",
" 78,\n",
" 89,\n",
" 88,\n",
" 87,\n",
" 97,\n",
" 107,\n",
" 106,\n",
" 105,\n",
" 117,\n",
" 125,\n",
" 126,\n",
" 127,\n",
" 137,\n",
" 147,\n",
" 148,\n",
" 149,\n",
" 150,\n",
" 151,\n",
" 152,\n",
" 153,\n",
" 154,\n",
" 155,\n",
" 139,\n",
" 135,\n",
" 134,\n",
" 133,\n",
" 132,\n",
" 131,\n",
" 130,\n",
" 129,\n",
" 118,\n",
" 109,\n",
" 110,\n",
" 111,\n",
"]\n",
"\n",
"\n",
"pm = generate_preset_pass_manager(\n",
" optimization_level=1, backend=backend, initial_layout=layout\n",
")\n",
"\n",
"isa_circuits_mid = pm.run(circuits_mid_hw)\n",
"isa_circuits = pm.run(circuits_hw)\n",
"\n",
"print(f\"Transpiled circuits. Example depth: {isa_circuits[0].depth()}\")\n",
"\n",
"# Define and layout-map observables\n",
"observables_hw = [\n",
" SparsePauliOp(\"I\" * i + \"Z\" + \"I\" * (num_qubits_hw - i - 1))\n",
" for i in range(num_qubits_hw)\n",
"]\n",
"\n",
"isa_observables_mid = [\n",
" [obs.apply_layout(isa_circuits_mid[i].layout) for obs in observables_hw]\n",
" for i in range(len(isa_circuits_mid))\n",
"]\n",
"isa_observables = [\n",
" [obs.apply_layout(isa_circuits[i].layout) for obs in observables_hw]\n",
" for i in range(len(isa_circuits))\n",
"]\n",
"\n",
"# Build PUBs — interleave meson and vacuum for each Trotter step\n",
"isa_pubs_mid = [\n",
" (circ, obs) for circ, obs in zip(isa_circuits_mid, isa_observables_mid)\n",
"]\n",
"isa_pubs = [(circ, obs) for circ, obs in zip(isa_circuits, isa_observables)]\n",
"\n",
"pubs_to_execute = [\n",
" [isa_pubs_mid[i], isa_pubs[i]] for i in range(len(isa_pubs))\n",
"]\n",
"\n",
"# -------------------------Step 3: Execute on hardware-------------------------\n",
"\n",
"twirling_options = TwirlingOptions(\n",
" enable_gates=True,\n",
" enable_measure=True,\n",
" shots_per_randomization=\"auto\",\n",
" strategy=\"active-circuit\",\n",
")\n",
"\n",
"resilience_options = ResilienceOptionsV2(\n",
" measure_mitigation=True, # TREX readout error mitigation\n",
" zne_mitigation=False, # ZNE turned off\n",
")\n",
"\n",
"dd_options = DynamicalDecouplingOptions(\n",
" enable=False # Circuit is sufficiently dense\n",
")\n",
"\n",
"options = EstimatorOptions(\n",
" resilience=resilience_options,\n",
" twirling=twirling_options,\n",
" dynamical_decoupling=dd_options,\n",
" default_shots=10_000,\n",
")\n",
"\n",
"ids = []\n",
"with Batch(backend=backend) as batch:\n",
" for idx, pub in enumerate(pubs_to_execute):\n",
" print(f\"Submitting job for Trotter step {idx + 1}\")\n",
" estimator = EstimatorV2(mode=batch, options=options)\n",
" estimator.skip_transpilation = True\n",
" job = estimator.run(pub)\n",
" ids.append(job.job_id())\n",
" batch_id = batch.session_id\n",
"\n",
"job_info = {\"ids\": ids, \"batch_id\": batch_id}\n",
"print(f\"Submitted {len(ids)} jobs. Batch ID: {batch_id}\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f03fb6e6-2ba7-49bb-b1f5-eb9b6eda993b",
"metadata": {},
"outputs": [],
"source": [
"print(ids)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "72d09009-e0d2-4bb0-9157-2e88a9d973ea",
"metadata": {},
"outputs": [],
"source": [
"# -------------------------Step 4: Post-process results-------------------------\n",
"\n",
"jobs = [service.job(job_id) for job_id in ids]\n",
"results = [job.result() for job in jobs]\n",
"\n",
"# Extract expectation values (index 0 = meson, index 1 = vacuum)\n",
"raw_expvals_mid_hw = [result[0].data.evs[::-1] for result in results]\n",
"raw_expvals_hw = [result[1].data.evs[::-1] for result in results]\n",
"\n",
"# Compute fermion numbers and differential\n",
"N_mid_hw = get_number(raw_expvals_mid_hw, num_lattice_point_hw)\n",
"N_hw = get_number(raw_expvals_hw, num_lattice_point_hw)\n",
"N_diff_hw = calculate_difference(N_mid_hw, N_hw, num_lattice_point_hw)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "19ee420d",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"N_diff_hw_arr = np.array(N_diff_hw)\n",
"\n",
"fig, ax = plt.subplots(figsize=(10, 6))\n",
"vmax = np.abs(N_hw).max()\n",
"vmin = -vmax\n",
"norm = TwoSlopeNorm(vmin=vmin, vcenter=0, vmax=vmax)\n",
"im = ax.imshow(\n",
" N_diff_hw_arr,\n",
" aspect=\"auto\",\n",
" origin=\"lower\",\n",
" cmap=\"RdBu_r\",\n",
" norm=norm,\n",
" extent=[0, 8, 0.5, len(trotter_steps_hw) + 0.5],\n",
")\n",
"ax.set_xlabel(\"Lattice site $r$\", fontsize=13)\n",
"ax.set_ylabel(\"Trotter step $t$\", fontsize=13)\n",
"ax.set_title(\n",
" \"Staggered Fermionic Occupation Number Dynamics\\nQuantum Simulation on IBM Hardware — 30-site lattice (60 qubits)\",\n",
" fontsize=13,\n",
")\n",
"cbar = plt.colorbar(im, ax=ax)\n",
"cbar.set_label(\"$n_f(r,t)$\", fontsize=12)\n",
"plt.tight_layout()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "de8e2aa6",
"metadata": {},
"source": [
"\n",
"\n",
"## Analyse comparative classique par propagation de Pauli\n",
"\n",
"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 « $R_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.\n",
"\n",
"Le déroulement des opérations avec [`pauli-prop`](https://github.com/Qiskit/pauli-prop) est le suivant :\n",
"\n",
"1. **Divisez** le circuit en ses parties « Clifford » et « non-Clifford » à l'aide de `evolve_through_cliffords`.\n",
"2. `atol`**Propager** 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.\n",
"3. **Faites évoluer** le résultat via la partie Clifford en utilisant la prise en charge intégrée de Clifford par Qiskit.\n",
"4. **On obtient** la valeur attendue en additionnant les coefficients des termes de Pauli diagonaux (qui ne contiennent que $I$ et $Z$ ).\n",
"\n",
"\n",
"\n",
"### Seuil de troncature\n",
"\n",
"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](https://arxiv.org/abs/2602.18080) '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.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 24,
"id": "0ed2dd40",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"PPM settings: atol=0.001, max_terms=66000\n",
"Trotter step 1: 5.0 s\n",
"Trotter step 2: 7.5 s\n",
"Trotter step 3: 11.2 s\n",
"Trotter step 4: 14.7 s\n",
"Trotter step 5: 18.3 s\n",
"Trotter step 6: 22.1 s\n",
"Trotter step 7: 25.6 s\n",
"Trotter step 8: 29.4 s\n",
"Trotter step 9: 33.2 s\n",
"Trotter step 10: 36.6 s\n",
"\n",
"Total PPM simulation time: 203.6 s\n",
"Truncation threshold used: 0.001\n"
]
}
],
"source": [
"import time\n",
"from pauli_prop import evolve_through_cliffords, propagate_through_circuit\n",
"\n",
"# ── PPM Configuration ──\n",
"# Truncation threshold: controls the speed/accuracy trade-off.\n",
"PPM_THRESHOLD = 1e-3\n",
"\n",
"# Maximum Pauli terms to track per observable (hard cap on memory/time)\n",
"PPM_MAX_TERMS = 66_000\n",
"\n",
"print(f\"PPM settings: atol={PPM_THRESHOLD}, max_terms={PPM_MAX_TERMS}\")\n",
"\n",
"# We propagate each single-qubit Z observable through each circuit.\n",
"# For PPM, we work with the un-transpiled circuits (ideal noiseless simulation).\n",
"\n",
"observables_pp = [\n",
" SparsePauliOp(\"I\" * i + \"Z\" + \"I\" * (num_qubits_hw - i - 1))\n",
" for i in range(num_qubits_hw)\n",
"]\n",
"\n",
"\n",
"def ppm_expectation_values(\n",
" circuit, observables, max_terms=PPM_MAX_TERMS, atol=PPM_THRESHOLD\n",
"):\n",
" \"\"\"Compute expectation values of single-qubit Z observables\n",
" via Pauli propagation.\n",
"\n",
" Args:\n",
" circuit: The quantum circuit to simulate.\n",
" observables: List of single-qubit Z observables.\n",
" max_terms: Maximum number of Pauli terms to retain (hard cap).\n",
" atol: Absolute tolerance — Pauli terms with coefficients below this\n",
" value are discarded during propagation. Larger values give\n",
" faster simulation at the cost of approximation accuracy.\n",
" \"\"\"\n",
" circuit = circuit.decompose([\"swap\"]) # decompose SWAPs into 3 CX gates\n",
" cliff, non_cliff = evolve_through_cliffords(circuit)\n",
"\n",
" evs = []\n",
" for obs in observables:\n",
" evolved_obs = propagate_through_circuit(\n",
" obs, non_cliff, max_terms=max_terms, atol=atol, frame=\"h\"\n",
" )[0]\n",
" evolved_obs.paulis = evolved_obs.paulis.evolve(cliff, frame=\"h\")\n",
" diagonal_mask = ~evolved_obs.paulis.x.any(axis=1)\n",
" ev = float(evolved_obs.coeffs[diagonal_mask].sum().real)\n",
" evs.append(ev)\n",
" return np.array(evs)\n",
"\n",
"\n",
"# Run PPM for each Trotter step and record wall-clock time\n",
"pp_expvals_mid = []\n",
"pp_expvals = []\n",
"pp_times = []\n",
"\n",
"for idx, d in enumerate(trotter_steps_hw):\n",
" t_start = time.perf_counter()\n",
"\n",
" # Meson circuit\n",
" evs_mid = ppm_expectation_values(circuits_mid_hw[idx], observables_pp)\n",
"\n",
" # Vacuum circuit\n",
" evs_vac = ppm_expectation_values(circuits_hw[idx], observables_pp)\n",
"\n",
" elapsed = time.perf_counter() - t_start\n",
" pp_times.append(elapsed)\n",
"\n",
" pp_expvals_mid.append(evs_mid[::-1])\n",
" pp_expvals.append(evs_vac[::-1])\n",
"\n",
" print(f\"Trotter step {d:2d}: {elapsed:.1f} s\")\n",
"\n",
"print(f\"\\nTotal PPM simulation time: {sum(pp_times):.1f} s\")\n",
"print(f\"Truncation threshold used: {PPM_THRESHOLD}\")"
]
},
{
"cell_type": "code",
"execution_count": 25,
"id": "pauli-prop-timing-plot",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# --- PPM simulation time vs. Trotter steps ---\n",
"fig, ax = plt.subplots(figsize=(8, 5))\n",
"ax.plot(\n",
" list(trotter_steps_hw),\n",
" pp_times,\n",
" \"o-\",\n",
" color=\"tab:blue\",\n",
" linewidth=2,\n",
" markersize=6,\n",
")\n",
"ax.set_xlabel(\"Trotter step\", fontsize=13)\n",
"ax.set_ylabel(\"Wall-clock time (s)\", fontsize=13)\n",
"ax.set_title(\n",
" \"Pauli Propagation simulation time vs. Trotter steps\\n(30-site lattice, 60 qubits)\",\n",
" fontsize=13,\n",
")\n",
"ax.grid(True, alpha=0.3)\n",
"plt.tight_layout()\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 37,
"id": "pauli-prop-heatmap",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# --- PPM heatmap and comparison with hardware ---\n",
"N_mid_pp = get_number(pp_expvals_mid, num_lattice_point_hw)\n",
"N_pp = get_number(pp_expvals, num_lattice_point_hw)\n",
"N_diff_pp = calculate_difference(N_mid_pp, N_pp, num_lattice_point_hw)\n",
"\n",
"N_diff_pp_arr = np.array(N_diff_pp)\n",
"\n",
"fig, axes = plt.subplots(1, 2, figsize=(18, 6))\n",
"vmax = np.abs(N_hw).max()\n",
"vmin = -vmax\n",
"norm = TwoSlopeNorm(vmin=vmin, vcenter=0, vmax=vmax)\n",
"\n",
"# PPM result\n",
"im0 = axes[0].imshow(\n",
" N_diff_pp_arr,\n",
" aspect=\"auto\",\n",
" origin=\"lower\",\n",
" cmap=\"RdBu_r\",\n",
" norm=norm,\n",
" extent=[0, num_lattice_point_hw, 0.5, len(trotter_steps_hw) + 0.5],\n",
")\n",
"axes[0].set_xlabel(\"Lattice site $r$\", fontsize=12)\n",
"axes[0].set_ylabel(\"Trotter step $t$\", fontsize=12)\n",
"axes[0].set_title(\n",
" \"Pauli Propagation\\n(classical noiseless simulation)\", fontsize=12\n",
")\n",
"plt.colorbar(im0, ax=axes[0], label=\"$n_f(r,t)$\")\n",
"\n",
"# Hardware result\n",
"im1 = axes[1].imshow(\n",
" N_diff_hw_arr,\n",
" aspect=\"auto\",\n",
" origin=\"lower\",\n",
" cmap=\"RdBu_r\",\n",
" norm=norm,\n",
" extent=[0, num_lattice_point_hw, 0.5, len(trotter_steps_hw) + 0.5],\n",
")\n",
"axes[1].set_xlabel(\"Lattice site $r$\", fontsize=12)\n",
"axes[1].set_ylabel(\"Trotter step $t$\", fontsize=12)\n",
"axes[1].set_title(\n",
" \"Quantum Simulation\\n(IBM Hardware, readout error mitigation only)\",\n",
" fontsize=12,\n",
")\n",
"plt.colorbar(im1, ax=axes[1], label=\"$n_f(r,t)$\")\n",
"\n",
"plt.suptitle(\n",
" \"Staggered Fermionic Occupation Number Dynamics — 30-site lattice\",\n",
" fontsize=14,\n",
" y=1.02,\n",
")\n",
"plt.tight_layout()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "next-steps",
"metadata": {},
"source": [
"\n",
"\n",
"## Etapes suivantes\n",
"\n",
"Si ce travail vous a intéressé, n'hésitez pas à consulter les ressources suivantes :\n",
"\n",
"\n",
" * [Documentation sur la primitive « Estimator » de Qiskit](/docs/guides/get-started-with-estimator) — pour plus de détails sur la configuration des options d'atténuation des erreurs\n",
" * [Techniques d'atténuation et de suppression des erreurs](/docs/guides/error-mitigation-and-suppression-techniques) — pour en savoir plus sur TREX, ZNE et d'autres méthodes d'atténuation\n",
" * [Qiskit Pauli Propagation (pauli-prop)](https://github.com/Qiskit/pauli-prop) — Simulation classique accélérée par Rust via la rétropropagation de Pauli\n",
"\n",
"\n",
"\n",
"\n",
"## Références\n",
"\n",
"\\[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](https://arxiv.org/abs/2602.18080) (2026)\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"id": "a1b8767d",
"source": "© IBM Corp., 2017-2026"
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3"
},
"hours": 1,
"qpuSeconds": 360
},
"nbformat": 4,
"nbformat_minor": 5
}