{
"cells": [
{
"cell_type": "markdown",
"id": "title-cell",
"metadata": {},
"source": [
"---\n",
"title: \"Osservazione di una dinamica adronica non abeliana robusta e coerente su processori quantistici soggetti a rumore\"\n",
"description: \"Simulare la dinamica degli adroni nella teoria di gauge su reticolo SU(2) utilizzando il framework LSH su un hardware quant IBM.\"\n",
"---\n",
"\n",
"{/* cspell:ignore Kogut Susskind expvals Pstep Nstep vmax vmin imshow fontsize cbar Ilčić */}\n",
"\n",
"\n",
"\n",
"# Osservazione di una dinamica adronica non abeliana robusta e coerente su processori quantistici soggetti a rumore\n",
"\n",
"*Stima del tempo di esecuzione: 6 minuti su un processore Heron (ibm\\_boston o equivalente) (NOTA: Si tratta solo di una stima. (La durata potrebbe variare.)*\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "learning-outcomes",
"metadata": {},
"source": [
"\n",
"\n",
"## Risultati di apprendimento\n",
"\n",
"Al termine di questo tutorial, avrai appreso quanto segue:\n",
"\n",
"* Come le teorie di gauge su reticoli non abeliani (in particolare SU(2)) possano essere riformulate utilizzando il quadro Loop-String-Hadron (LSH) per una simulazione quantistica efficiente\n",
"* Come costruire circuiti di evoluzione temporale di tipo Trotter per un hamiltoniano approssimativo di una teoria di gauge SU(2) e mapparli su qubit\n",
"* Come eseguire questi circuiti su un hardware d IBM Quantum® e utilizzando la primitiva Qiskit Estimator con mitigazione degli errori di lettura\n",
"\n",
"\n",
"\n",
"## Prerequisiti\n",
"\n",
"Si consiglia di approfondire i seguenti argomenti:\n",
"\n",
"* [Nozioni di base sui circuiti e sulle porte quantistiche](/learning/courses/basics-of-quantum-information)\n",
"* [Introduzione alla primitiva Estimator di Qiskit](/docs/guides/get-started-with-estimator)\n",
"* Conoscenza di base dei concetti della teoria quantistica dei campi (utile ma non indispensabile; la sezione dedicata alle nozioni di base ne illustra gli elementi essenziali)\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "background",
"metadata": {},
"source": [
"\n",
"\n",
"## Sfondo\n",
"\n",
"\n",
"\n",
"### Motivazione\n",
"\n",
"La cromodinamica quantistica (QCD), la teoria di gauge SU(3) della forza forte, lega i quark in adroni e regola il confinamento e la rottura delle stringhe. I metodi classici della QCD su reticolo eccellono nell'analisi delle proprietà statiche, ma non sono in grado di simulare le dinamiche in tempo reale a causa del problema del segno. I computer quantistici offrono una via per aggirare questa barriera codificando i gradi di libertà dei campi di gauge direttamente sui qubit.\n",
"\n",
"Questo tutorial illustra una simulazione di questo tipo: utilizza l'hardware \" IBM Quantum \" per simulare la propagazione degli adroni in tempo reale in una teoria di gauge a reticolo SU(2) a (1+1) dimensioni — la teoria di gauge non abeliana più semplice e un primo passo verso la QCD completa.\n",
"\n",
"\n",
"\n",
"### L'hamiltoniano di Kogut-Susskind\n",
"\n",
"La teoria è formulata su un reticolo spaziale di tipo “ 1D ”, con fermioni sfalsati (materia) sui siti e campi di gauge SU(2) sui legami. Dopo averlo riportato in forma adimensionale, l'hamiltoniano è:\n",
"\n",
"$W = H_E^{\\text{(KS)}} + \\mu H_M + x H_I^{\\text{(KS)}},$\n",
"\n",
"dove $H_E$ rappresenta l'energia del campo cromoelettrico, $H_M$ è il termine di massa sfalsata, $H_I$ è il termine di interazione materia-calibro (hopping), $\\mu = 2\\frac{m}{g}\\sqrt{x}$ codifica la massa del fermione e $x = \\frac{1}{g^2 a^2}$ è l'intensità di interazione. Il limite al continuo della teoria è descritto all'indirizzo $N \\to \\infty$ e $x \\to \\infty$.\n",
"\n",
"\n",
"\n",
"### Il framework Loop-String-Hadron (LSH)\n",
"\n",
"Una delle principali difficoltà risiede nel fatto che lo spazio di Hilbert del campo di gauge su ciascun collegamento è di dimensione infinita. Il quadro teorico **Loop-String-Hadron (LSH)** affronta questo problema riformulando la teoria in termini di variabili invarianti di gauge: loop di flusso, stringhe che collegano cariche separate e adroni (coppie di fermioni singolette di gauge in un sito). Nella base LSH, la legge di Gauss è soddisfatta automaticamente per costruzione, quindi ogni stato di base è fisico. Ogni sito del reticolo è caratterizzato da tre numeri quantici $(n_l, n_i, n_o)$ che rappresentano il numero di anello, la stringa in entrata e la stringa in uscita, dove $n_i, n_o \\in \\{0,1\\}$ sono fermionici e $n_l \\geq 0$ è bosonico. Da queste espressioni si definisce il numero locale di fermioni come $n_f(r) = n_i(r) + n_o(r)$ per i siti pari e $n_f(r) = 2 - [n_i(r) + n_o(r)]$ per i siti dispari.\n",
"\n",
"\n",
"\n",
"### Dall’Hamiltoniano completo al circuito quantistico: tre approssimazioni fondamentali\n",
"\n",
"Il circuito quantistico **non** simula esattamente l'Hamiltoniano SU(2) completo. Al contrario, implementa una serie controllata di approssimazioni valide nel **regime** di accoppiamento debole ( $x \\gg 1$ ). È fondamentale comprendere cosa viene approssimato e cosa no:\n",
"\n",
"**Approssimazione 1 — Limite di accoppiamento debole per l’ $H_I$ o:** l’Hamiltoniano a interazione completa $H_I^{\\text{(LSH)}}$ (Eq. 16 in [\\[1\\]](#references) ) contiene prefattori che dipendono dal numero quantico bosonico $n_l$ tramite termini del tipo $1/\\sqrt{n_l+1}$. Nel regime di accoppiamento debole ( $x \\gg 1$ ), la dinamica è dominata dal termine elettrico $H_E$, che favorisce stati con un valore elevato di $n_l$. Per $n_l \\gg 1$, il rapporto $n_l/(n_l+1) \\to 1$ e tutti questi prefattori si riducono all’unità. L'Hamiltoniano di interazione si riduce quindi a un salto tra vicini più prossimi di natura puramente locale:\n",
"\n",
"$H_I^{\\text{approx}} = -\\sum_r \\left[\\sigma^-(r)\\sigma^+(r+1) + \\sigma^+(r)\\sigma^-(r+1)\\right],$\n",
"\n",
"che è indipendente dall’ $n_l$ e e agisce solo sui qubit fermionici $(n_i, n_o)$.\n",
"\n",
"**Approssimazione 2 — Flusso medio globale per $H_E$ :** L'energia elettrica dipende da $n_l$ in ciascun collegamento. Nel vuoto a accoppiamento debole, l' $n_l$ e è elevata e approssimativamente uniforme. Sostituire i valori di $n_l$, che dipendono dal sito, con un unico valore medio globale $\\bar{n}_l$, rendendo $H_E$ una fase diagonale proporzionale alla configurazione dei fermioni in ciascun sito:\n",
"\n",
"$H_E^{\\text{approx}} = N h_E^0 + \\sum_{\\{r'\\}} \\left(\\frac{\\bar{n}_l}{2} + \\frac{3}{4}\\right)$\n",
"\n",
"dove $\\{r'\\}$ si somma su tutti i siti nella configurazione fermionica $(n_i=0, n_o=1)$, e $h_E^0$ è una fase globale che si può ignorare.\n",
"\n",
"**Approssimazione 3 — Trotterizzazione:** l'operatore di evoluzione temporale per un passo di durata $\\delta_\\tau$ si scompone come segue:\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",
"dove $c = \\delta_\\tau x$, $\\tilde{m} = \\delta_\\tau \\mu$ e $\\theta = -\\delta_\\tau(\\bar{n}_l/2 + 3/4)$. Questa decomposizione di Trotter di primo ordine introduce un errore che si annulla quando $\\delta_\\tau \\to 0$. Si fissa $\\delta_\\tau = 0.0015$ per tutto il calcolo.\n",
"\n",
"**Il risultato** di queste tre approssimazioni è che solo i due qubit fermionici per sito $(n_i, n_o)$ sono dinamici — il grado di libertà bosonico $n_l$ è stato assorbito nei parametri effettivi. Si ottiene così un circuito compatto con $2N$ qubit per $N$ siti del reticolo, in cui ogni passo di Trotter presenta una profondità costante di porte a due qubit (13 per passo).\n",
"\n",
"\n",
"\n",
"### Cosa simula questo tutorial\n",
"\n",
"Il tutorial simula **la propagazione degli adroni** : partendo dal vuoto a forte accoppiamento (uno stato di prodotto), si posiziona un mesone al centro del reticolo e si osserva l'evoluzione nel tempo. Il protocollo di misurazione differenziale — che prevede di far funzionare il circuito con e senza il mesone centrale, per poi effettuare la sottrazione — isola il segnale adronico coerente sia dal rumore dell'hardware che dagli effetti di confine. Il risultato è un motivo a cono di luce costituito da oscillazioni della densità dei fermioni, caratteristico di una modalità di respirazione mesonica confinata.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "requirements",
"metadata": {},
"source": [
"\n",
"\n",
"## Requisiti\n",
"\n",
"Prima di iniziare questo tutorial, installa quanto segue:\n",
"\n",
"* Qiskit SDK v2.0 o versioni successive, con supporto [alla visualizzazione](/docs/api/qiskit/visualization)\n",
"* Qiskit Runtime v0.22 o versioni successive (`pip install qiskit-ibm-runtime`)\n",
"* Pacchetto Pauli Propagation (`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",
"## Configura\n",
"\n",
"Inizia importando le librerie necessarie e definendo le funzioni di supporto che costruiscono i circuiti quantistici per l'evoluzione temporale LSH. Esistono tre funzioni fondamentali per la creazione di circuiti:\n",
"\n",
"1. **`pair_hamiltonian_circuit`**: Implementa l' $U_I$ e unitaria a due qubit per l'Hamiltoniano di interazione approssimativo tra siti adiacenti. La scomposizione in porte è la seguente: $\\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`**: Implementa l' $U_E$ e unitaria a due qubit per l'energia approssimativa del campo elettrico in ciascun sito. La scomposizione in porte è la seguente: $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`**: Realizza il circuito \"Trotterizzato\" completo, sovrapponendo i termini di interazione, elettrici e di massa con porte SWAP per gestire la connettività dei qubit.\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",
"## Esempio di simulatore su piccola scala\n",
"\n",
"In primo luogo, illustra il flusso di lavoro su piccola scala utilizzando un reticolo a sei siti (12 qubit), in modo da poter verificare la costruzione del circuito e comprendere le grandezze fisiche osservabili prima di eseguire il codice sull'hardware.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "step1-header",
"metadata": {},
"source": [
"\n",
"\n",
"### Fase 1: Mappare gli input classici su un problema quantistico\n",
"\n",
"Definire i parametri fisici corrispondenti al regime di accoppiamento debole studiato nell'articolo ( $x = 100$, $m/g = 1$ ). I parametri del circuito ricavati sono:\n",
"\n",
"* $c = \\delta_\\tau \\cdot x = 0.15$ (parametro di interazione)\n",
"* $\\theta = -\\delta_\\tau (\\bar{n}_l/2 + 3/4) = 0.01$ (fase del campo elettrico)\n",
"* $\\tilde{m} = \\delta_\\tau \\cdot \\mu = 0.03$ (parametro di massa)\n",
"\n",
"Per ogni conteggio di passi di Trotter, si costruiscono **due circuiti** : uno che inizializza un mesone al centro (`inverse_mid=True`) e uno che prepara il vuoto a forte accoppiamento (`inverse_mid=False`). Il protocollo di misurazione differenziale sottrae l'evoluzione del vuoto per isolare il segnale adronico.\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",
"### Fase 2: Ottimizzare il problema per l'esecuzione su hardware quantistico\n",
"\n",
"Definire le grandezze osservabili: misurazioni di tipo “ $Z$ ” su un singolo qubit per ciascun qubit. Da $\\langle Z \\rangle$ è possibile ricavare le probabilità di occupazione e quindi il numero di fermioni sfalsato $n_f(r)$ in ciascun sito del reticolo $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",
"### Passaggio 3: Eseguire il comando utilizzando Qiskit primitives\n",
"\n",
"Utilizzare `StatevectorEstimator` per una simulazione esatta e priva di rumore su piccola scala.\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",
"### Fase 4: Elaborazione finale e restituzione del risultato nel formato classico desiderato\n",
"\n",
"Convertire i valori attesi nel numero di fermioni sfalsato $n_f(r, t)$ e applicare il protocollo di misurazione differenziale (mesone $-$ vuoto) per generare la mappa termica della propagazione degli adroni. Questo grafico riproduce la struttura della Figura 3 dell'articolo di riferimento: il sito reticolare $r$ sull'asse x, il passo di Trotter (tempo) $t$ sull'asse y e $n_f(r,t)$ come scala cromatica.\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",
"## Esempio di hardware su larga scala\n",
"\n",
"Ora passiamo a un reticolo di 30 siti (60 qubit) su un hardware dell IBM Quantum. A questa scala, il circuito a 10 passi di Trotter comprende oltre 3.400 porte a due qubit e 14.000 porte a un qubit.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "hardware-steps",
"metadata": {},
"source": [
"\n",
"\n",
"### Passaggi da 1 a 4 (raggruppati in un unico blocco di codice)\n",
"\n",
"Aspetti chiave del flusso di lavoro hardware:\n",
"\n",
"* 10 passi di Trotter per i circuiti del mesone e del vuoto (intercalati per ridurre al minimo la deriva)\n",
"* Trasposizione con `optimization_level=1` — il layout del circuito è già isomorfo alla topologia del dispositivo (una catena lineare), quindi non sono necessari SWAP di instradamento. Il transpiler viene utilizzato esclusivamente per selezionare una catena di qubit fisici a basso rumore e per scomporre i gate nel set di gate nativo.\n",
"* `EstimatorV2` con la mitigazione degli errori di lettura TREX e il “Pauli twirling”\n",
"* `Batch` sessione per inviare tutti i lavori contemporaneamente\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",
"## Benchmarking classico tramite la propagazione di Pauli\n",
"\n",
"Il metodo di propagazione di Pauli (PPM) fornisce una simulazione classica priva di rumore del circuito quantistico, propagando a ritroso le grandezze osservabili misurate attraverso il circuito nella rappresentazione di Heisenberg. Negli strati di Clifford (porte CNOT, H, S, X), gli operatori di Pauli si mappano su altri operatori di Pauli senza aumentare il numero di termini. Gli strati non-Clifford (le porte \" $R_z$ \" presenti nel circuito) possono causare ramificazioni — nel peggiore dei casi, raddoppiando il numero di termini — ma molte ramificazioni hanno coefficienti piccoli e possono essere troncate.\n",
"\n",
"Il flusso di lavoro con [`pauli-prop`](https://github.com/Qiskit/pauli-prop) è il seguente:\n",
"\n",
"1. **Dividi** il circuito nelle sue parti Clifford e non Clifford utilizzando `evolve_through_cliffords`.\n",
"2. `atol`**Propagare** ciascun osservabile attraverso la parte non-Clifford utilizzando `propagate_through_circuit`, mantenendo fino a `max_terms` termini di Pauli ed eliminando i termini con coefficienti inferiori alla soglia di troncamento.\n",
"3. **Evolvi** il risultato attraverso la parte Clifford utilizzando il supporto integrato di Qiskit per Clifford.\n",
"4. **Si ricava** il valore atteso sommando i coefficienti dei termini di Pauli diagonali (che contengono solo $I$ e $Z$ ).\n",
"\n",
"\n",
"\n",
"### Soglia di troncamento\n",
"\n",
"Il `atol` parametro in `propagate_through_circuit` determina l'intensità con cui vengono eliminati i rami di Pauli di piccole dimensioni. Una soglia molto stretta (ad esempio, `1e-12`) conserva quasi tutti i rami e fornisce risultati esatti, ma il tempo di simulazione cresce rapidamente con la profondità del circuito; la simulazione a 120 qubit descritta [nell](https://arxiv.org/abs/2602.18080) ’articolo ha richiesto circa 8.5 ore con le impostazioni predefinite. Alzando la soglia (ad esempio, a `1e-6` o `1e-3`) si scartano i termini i cui coefficienti sono inferiori a tale valore, riducendo drasticamente il numero di termini monitorati e velocizzando il calcolo. Il compromesso consiste in un errore di approssimazione minimo e controllabile, che è possibile verificare confrontando i risultati ottenuti con soglie diverse.\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",
"## Passi successivi\n",
"\n",
"Se questo lavoro ti è sembrato interessante, ti invitiamo a dare un'occhiata al seguente materiale:\n",
"\n",
"\n",
" * [Documentazione sulle primitive di Qiskit Estimator](/docs/guides/get-started-with-estimator) — per ulteriori dettagli sulla configurazione delle opzioni di mitigazione degli errori\n",
" * [Tecniche di mitigazione e soppressione degli errori](/docs/guides/error-mitigation-and-suppression-techniques) — per saperne di più su TREX, ZNE e altri metodi di mitigazione\n",
" * [Qiskit Pauli Propagation (pauli-prop)](https://github.com/Qiskit/pauli-prop) — Simulazione classica accelerata da Rust tramite retropropagazione di Pauli\n",
"\n",
"\n",
"\n",
"\n",
"## Riferimenti\n",
"\n",
"\\[1] L'articolo originale: Ilčić, Majumdar, Mathew et al. \"Osservazione di dinamiche adroniche non abeliane robuste e coerenti su processori quantistici soggetti a rumore\" [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
}