{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "title-cell",
      "metadata": {},
      "source": [
        "---\n",
        "title: \"Observación de una dinámica hadrónica no abeliana robusta y coherente en procesadores cuánticos con ruido\"\n",
        "description: \"Simular la dinámica de hadrones de la teoría de gauge en red SU(2) utilizando el marco LSH en un hardware cuántico de tip IBM.\"\n",
        "---\n",
        "\n",
        "{/* cspell:ignore Kogut Susskind expvals Pstep Nstep vmax vmin imshow fontsize cbar Ilčić */}\n",
        "\n",
        "<span id=\"observation-of-robust-and-coherent-non-abelian-hadron-dynamics-on-noisy-quantum-processors\" />\n",
        "\n",
        "# Observación de una dinámica hadrónica no abeliana robusta y coherente en procesadores cuánticos con ruido\n",
        "\n",
        "*Estimación de tiempo de ejecución: 6 minutos en un procesador Heron (ibm\\_boston o equivalente) (NOTA: Se trata únicamente de una estimación. (El tiempo de ejecución puede variar.)*\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "learning-outcomes",
      "metadata": {},
      "source": [
        "<span id=\"learning-outcomes\" />\n",
        "\n",
        "## Resultados del aprendizaje\n",
        "\n",
        "Al finalizar este tutorial, habrás aprendido lo siguiente:\n",
        "\n",
        "* Cómo se pueden reformular las teorías de gauge en retículos no abelianos (concretamente SU(2)) utilizando el marco Loop-String-Hadron (LSH) para una simulación cuántica eficiente\n",
        "* Cómo construir circuitos de evolución temporal «trotterizados» para un hamiltoniano aproximado de una teoría de gauge SU(2) y mapearlos en qubits\n",
        "* Cómo ejecutar estos circuitos en un hardware de tip IBM Quantum® utilizando la primitiva «Qiskit Estimator» con mitigación de errores de lectura\n",
        "\n",
        "<span id=\"prerequisites\" />\n",
        "\n",
        "## Requisitos previos\n",
        "\n",
        "Se recomienda que te familiarices con los siguientes temas:\n",
        "\n",
        "* [Conceptos básicos sobre circuitos y puertas cuánticas](/learning/courses/basics-of-quantum-information)\n",
        "* [Introducción a la primitiva «Estimator» de Qiskit](/docs/guides/get-started-with-estimator)\n",
        "* Conocimientos básicos sobre los conceptos de la teoría cuántica de campos (resulta útil, pero no es imprescindible; en la sección de antecedentes se tratan los aspectos fundamentales)\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "background",
      "metadata": {},
      "source": [
        "<span id=\"background\" />\n",
        "\n",
        "## En segundo plano\n",
        "\n",
        "<span id=\"motivation\" />\n",
        "\n",
        "### Motivación\n",
        "\n",
        "La cromodinámica cuántica (QCD), la teoría de gauge SU(3) de la fuerza fuerte, une a los quarks en hadrones y rige el confinamiento y la ruptura de cuerdas. Los métodos clásicos de la QCD de red destacan en el estudio de las propiedades estáticas, pero no pueden simular la dinámica en tiempo real debido al problema del signo. Los ordenadores cuánticos ofrecen una vía para sortear esta barrera mediante la codificación de los grados de libertad de los campos de gauge directamente en los qubits.\n",
        "\n",
        "Este tutorial muestra una simulación de este tipo: utiliza el hardware d IBM Quantum para simular la propagación de hadrones en tiempo real en una teoría de gauge en red SU(2) de dimensión (1+1) —la teoría de gauge no abeliana más simple y un paso previo hacia la QCD completa—.\n",
        "\n",
        "<span id=\"the-kogut-susskind-hamiltonian\" />\n",
        "\n",
        "### El hamiltoniano de Kogut-Susskind\n",
        "\n",
        "La teoría se formula en una red espacial de tipo « 1D », con fermiones (materia) distribuidos de forma escalonada en los vértices y campos de gauge SU(2) en los enlaces. Tras expresar el hamiltoniano en forma adimensional, este queda así:\n",
        "\n",
        "$W = H_E^{\\text{(KS)}} + \\mu H_M + x H_I^{\\text{(KS)}},$\n",
        "\n",
        "donde $H_E$ es la energía del campo cromoeléctrico, $H_M$ es el término de masa escalonada, $H_I$ es el término de interacción materia-calibre (salto), $\\mu = 2\\frac{m}{g}\\sqrt{x}$ codifica la masa del fermión y $x = \\frac{1}{g^2 a^2}$ es la intensidad de la interacción. El límite continuo de la teoría se encuentra en $N \\to \\infty$ y $x \\to \\infty$.\n",
        "\n",
        "<span id=\"the-loop-string-hadron-lsh-framework\" />\n",
        "\n",
        "### El marco Loop-String-Hadron (LSH)\n",
        "\n",
        "Uno de los principales retos es que el espacio de Hilbert del campo de gauge en cada enlace es de dimensión infinita. El marco **Loop-String-Hadron (LSH)** aborda esta cuestión reformulando la teoría en términos de variables invariantes de gauge: bucles de flujo, cuerdas que conectan cargas separadas y hadrones (pares de fermiones singlet de gauge en un sitio). En la base LSH, la ley de Gauss se cumple automáticamente por definición, por lo que todos los estados de la base son físicos. Cada sitio de la red se caracteriza por tres números cuánticos $(n_l, n_i, n_o)$, que representan el número de bucle, la cuerda entrante y la cuerda saliente, donde $n_i, n_o \\in \\{0,1\\}$ son fermiónicos y $n_l \\geq 0$ es bosónico. A partir de estos, el número de fermiones local se define como « $n_f(r) = n_i(r) + n_o(r)$ » para los sitios pares y « $n_f(r) = 2 - [n_i(r) + n_o(r)]$ » para los sitios impares.\n",
        "\n",
        "<span id=\"from-full-hamiltonian-to-the-quantum-circuit-three-key-approximations\" />\n",
        "\n",
        "### Del hamiltoniano completo al circuito cuántico: tres aproximaciones clave\n",
        "\n",
        "El circuito cuántico **no** simula con exactitud el hamiltoniano SU(2) completo. En su lugar, aplica una serie controlada de aproximaciones que son válidas en el **régimen de acoplamiento débil** ( $x \\gg 1$ ). Es fundamental comprender qué se aproxima y qué no:\n",
        "\n",
        "**Aproximación 1 — Límite de acoplamiento débil para « $H_I$ »:** El hamiltoniano de interacción completa $H_I^{\\text{(LSH)}}$ (ecuación 16 en [\\[1\\]](#references) ) contiene prefactores que dependen del número cuántico bosónico $n_l$ a través de términos como $1/\\sqrt{n_l+1}$. En el régimen de acoplamiento débil ( $x \\gg 1$ ), la dinámica viene dominada por el término eléctrico $H_E$, que favorece los estados con un $n_l$ elevado. Para $n_l \\gg 1$, la relación $n_l/(n_l+1) \\to 1$ y todos estos prefactores se simplifican a la unidad. El hamiltoniano de interacción se reduce entonces a un salto entre vecinos más cercanos de carácter puramente local:\n",
        "\n",
        "$H_I^{\\text{approx}} = -\\sum_r \\left[\\sigma^-(r)\\sigma^+(r+1) + \\sigma^+(r)\\sigma^-(r+1)\\right],$\n",
        "\n",
        "que es independiente de $n_l$ y actúa únicamente sobre los qubits fermiónicos $(n_i, n_o)$.\n",
        "\n",
        "**Aproximación 2 — Flujo medio global para « $H_E$ »:** La energía eléctrica depende de « $n_l$ » en cada enlace. En el vacío de acoplamiento débil, $n_l$ es grande y aproximadamente uniforme. Sustituye los valores de $n_l$, que dependen de cada sitio, por un único valor medio global $\\bar{n}_l$, de modo que $H_E$ sea una fase diagonal proporcional a la configuración de fermiones en cada sitio:\n",
        "\n",
        "$H_E^{\\text{approx}} = N h_E^0 + \\sum_{\\{r'\\}} \\left(\\frac{\\bar{n}_l}{2} + \\frac{3}{4}\\right)$\n",
        "\n",
        "donde $\\{r'\\}$ se suma por todos los sitios en la configuración fermiónica $(n_i=0, n_o=1)$, y $h_E^0$ es una fase global que puedes ignorar.\n",
        "\n",
        "**Aproximación 3 — Trotterización:** El operador de evolución temporal para un paso de duración $\\delta_\\tau$ se descompone de la siguiente manera:\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",
        "donde $c = \\delta_\\tau x$, $\\tilde{m} = \\delta_\\tau \\mu$ y $\\theta = -\\delta_\\tau(\\bar{n}_l/2 + 3/4)$. Esta descomposición de Trotter de primer orden introduce un error que se anula cuando $\\delta_\\tau \\to 0$. Fijamos $\\delta_\\tau = 0.0015$ en todo momento.\n",
        "\n",
        "**El resultado** de estas tres aproximaciones es que solo los dos qubits fermiónicos por sitio $(n_i, n_o)$ son dinámicos; el grado de libertad bosónico $n_l$ se ha integrado en los parámetros efectivos. Esto da como resultado un circuito compacto con $2N$ qubits para $N$ nodos de la red, en el que cada paso de Trotter tiene una profundidad constante de puertas de dos qubits (13 por paso).\n",
        "\n",
        "<span id=\"what-this-tutorial-simulates\" />\n",
        "\n",
        "### Qué simula este tutorial\n",
        "\n",
        "El tutorial simula **la propagación de hadrones** : partiendo del vacío de acoplamiento fuerte (un estado de producto), se coloca un mesón en el centro de la red y se deja que evolucione en el tiempo. El protocolo de medición diferencial —que consiste en hacer funcionar el circuito con y sin el mesón central y, a continuación, restar los resultados— aísla la señal coherente de los hadrones tanto del ruido del hardware como de los efectos de contorno. El resultado es un patrón en forma de cono de luz de oscilaciones en la densidad de fermiones, característico de un modo de respiración de mesones confinados.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "requirements",
      "metadata": {},
      "source": [
        "<span id=\"requirements\" />\n",
        "\n",
        "## Requisitos\n",
        "\n",
        "Antes de empezar este tutorial, instala lo siguiente:\n",
        "\n",
        "* Qiskit SDK v2.0 o posterior, con soporte [para visualización](/docs/api/qiskit/visualization)\n",
        "* Qiskit Runtime v0.22 o posterior (`pip install qiskit-ibm-runtime`)\n",
        "* Paquete de propagación 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": [
        "<span id=\"setup\" />\n",
        "\n",
        "## Configuración\n",
        "\n",
        "Empieza importando las bibliotecas necesarias y definiendo las funciones auxiliares que construyen los circuitos cuánticos para la evolución temporal LSH. Hay tres funciones básicas para la construcción de circuitos:\n",
        "\n",
        "1. **`pair_hamiltonian_circuit`**: Implementa el operador unitario de dos qubits $U_I$ para el hamiltoniano de interacción aproximado entre sitios vecinos. La descomposición de la puerta es: $\\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 el operador unitario de dos qubits $U_E$ para la energía aproximada del campo eléctrico en cada sitio. La descomposición de la puerta es: $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`**: Monta el circuito «Trotterizado» completo, combinando los términos de interacción, eléctricos y de masa con puertas SWAP para gestionar la conectividad de los 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 <Z> = p(0) - p(1), the occupation probability is p(1) = (1 - <Z>) / 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 <n_i> and <n_o>.\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": [
        "<span id=\"small-scale-simulator-example\" />\n",
        "\n",
        "## Ejemplo de simulador a pequeña escala\n",
        "\n",
        "En primer lugar, muestra el flujo de trabajo a pequeña escala utilizando una red de seis sitios (12 qubits), de modo que puedas verificar la construcción del circuito y comprender los observables físicos antes de ejecutarlo en el hardware.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "step1-header",
      "metadata": {},
      "source": [
        "<span id=\"step-1-map-classical-inputs-to-a-quantum-problem\" />\n",
        "\n",
        "### Paso 1: Asignar entradas clásicas a un problema cuántico\n",
        "\n",
        "Defina los parámetros físicos que se ajustan al régimen de acoplamiento débil estudiado en el artículo ( $x = 100$, $m/g = 1$ ). Los parámetros del circuito obtenidos son:\n",
        "\n",
        "* $c = \\delta_\\tau \\cdot x = 0.15$ (parámetro de interacción)\n",
        "* $\\theta = -\\delta_\\tau (\\bar{n}_l/2 + 3/4) = 0.01$ (fase del campo eléctrico)\n",
        "* $\\tilde{m} = \\delta_\\tau \\cdot \\mu = 0.03$ (parámetro de masa)\n",
        "\n",
        "Para cada recuento de pasos de Trotter, se construyen **dos circuitos** : uno que inicializa un mesón en el centro (`inverse_mid=True`) y otro que prepara el vacío de acoplamiento fuerte (`inverse_mid=False`). El protocolo de medición diferencial resta la evolución del vacío para aislar la señal del hadrón.\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": [
              "<Image src=\"/docs/images/tutorials/loop-string-hadron-dynamics/extracted-outputs/step1-circuits-1.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "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": [
        "<span id=\"step-2-optimize-problem-for-quantum-hardware-execution\" />\n",
        "\n",
        "### Paso 2: Optimizar el problema para su ejecución en hardware cuántico\n",
        "\n",
        "Definir las magnitudes observables: mediciones de « $Z$ » de un solo qubit en cada qubit. En $\\langle Z \\rangle$ se pueden extraer las probabilidades de ocupación y, a continuación, el número de fermiones escalonado $n_f(r)$ en cada sitio de la red $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": [
        "<span id=\"step-3-execute-using-qiskit-primitives\" />\n",
        "\n",
        "### Paso 3: Ejecutar utilizando Qiskit primitives\n",
        "\n",
        "Utilízalo `StatevectorEstimator` para realizar simulaciones exactas y sin ruido a pequeña escala.\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": [
        "<span id=\"step-4-post-process-and-return-result-in-desired-classical-format\" />\n",
        "\n",
        "### Paso 4: Realizar el posprocesamiento y obtener el resultado en el formato clásico deseado\n",
        "\n",
        "Convertir los valores esperados al número de fermiones escalonado $n_f(r, t)$ y aplicar el protocolo de medición diferencial (mesón $-$ vacío) para generar el mapa de calor de propagación de los hadrones. Esto reproduce la estructura de la figura 3 del artículo de referencia: el sitio de la red $r$ en el eje x, el paso de Trotter (tiempo) $t$ en el eje y, y $n_f(r,t)$ como escala de colores.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "id": "step4-postprocess",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/loop-string-hadron-dynamics/extracted-outputs/step4-postprocess-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "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": [
        "<span id=\"large-scale-hardware-example\" />\n",
        "\n",
        "## Ejemplo de hardware a gran escala\n",
        "\n",
        "Ahora ampliamos a una red de 30 sitios (60 qubits) en un hardware de tip IBM Quantum. A esta escala, el circuito de 10 pasos de Trotter consta de más de 3.400 puertas de dos qubits y 14.000 puertas de un solo qubit.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "hardware-steps",
      "metadata": {},
      "source": [
        "<span id=\"steps-1-4-compressed-into-a-single-code-block\" />\n",
        "\n",
        "### Pasos 1-4 (agrupados en un único bloque de código)\n",
        "\n",
        "Aspectos clave del flujo de trabajo de hardware:\n",
        "\n",
        "* 10 pasos de Trotter para los circuitos del mesón y del vacío (intercalados para minimizar la deriva)\n",
        "* Transpilación con `optimization_level=1` — el diseño del circuito ya es isomórfico a la topología del dispositivo (una cadena lineal), por lo que no se necesitan SWAP de enrutamiento. El transpilador se utiliza exclusivamente para seleccionar una cadena de qubits físicos con bajo nivel de ruido y descomponer las puertas en el conjunto de puertas nativas.\n",
        "* `EstimatorV2` con mitigación de errores de lectura de TREX y giro de Pauli\n",
        "* `Batch` sesión para enviar todos los trabajos a la vez\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": [
              "<Image src=\"/docs/images/tutorials/loop-string-hadron-dynamics/extracted-outputs/19ee420d-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "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": [
        "<span id=\"classical-benchmarking-via-pauli-propagation\" />\n",
        "\n",
        "## Evaluación comparativa clásica mediante la propagación de Pauli\n",
        "\n",
        "El método de propagación de Pauli (PPM) ofrece una simulación clásica sin ruido del circuito cuántico mediante la retropropagación de los observables medidos a lo largo del circuito en el marco de Heisenberg. En las capas de Clifford (puertas CNOT, H, S y X), los operadores de Pauli se mapean a otros operadores de Pauli sin aumentar el número de términos. Las capas que no son de Clifford (las puertas « $R_z$ » del circuito) pueden provocar ramificaciones —en el peor de los casos, duplicando el número de términos—, pero muchas ramificaciones tienen coeficientes pequeños y pueden truncarse.\n",
        "\n",
        "El proceso con [`pauli-prop`](https://github.com/Qiskit/pauli-prop) es el siguiente:\n",
        "\n",
        "1. **Divide** el circuito en sus partes de Clifford y no Clifford utilizando `evolve_through_cliffords`.\n",
        "2. `atol`**Propaga** cada observable a través de la parte no-Clifford utilizando `propagate_through_circuit`, conservando hasta `max_terms` términos de Pauli y descartando los términos con coeficientes inferiores al umbral de truncamiento.\n",
        "3. **Calcula** el resultado mediante la operación de Clifford utilizando la compatibilidad integrada de Qiskit con las operaciones de Clifford.\n",
        "4. **Se obtiene** el valor esperado sumando los coeficientes de los términos diagonales de Pauli (que contienen únicamente $I$ y $Z$ ).\n",
        "\n",
        "<span id=\"truncation-threshold\" />\n",
        "\n",
        "### Umbral de truncamiento\n",
        "\n",
        "El `atol` parámetro determina la intensidad con la que se podan las ramas pequeñas de `propagate_through_circuit` Pauli. Un umbral muy ajustado (por ejemplo, `1e-12`) conserva casi todas las ramificaciones y ofrece resultados exactos, pero el tiempo de simulación aumenta considerablemente con la profundidad del circuito; la simulación de 120 qubits que se presenta en el [artículo](https://arxiv.org/abs/2602.18080) tardó aproximadamente 8.5 horas con la configuración predeterminada. Al elevar el umbral (por ejemplo, a `1e-6` o `1e-3`), se descartan los términos cuyos coeficientes son inferiores a ese valor, lo que reduce drásticamente el número de términos analizados y agiliza el cálculo. La contrapartida es un pequeño error de aproximación controlable que puedes comprobar comparando los resultados con distintos umbrales.\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": [
              "<Image src=\"/docs/images/tutorials/loop-string-hadron-dynamics/extracted-outputs/pauli-prop-timing-plot-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "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": [
              "<Image src=\"/docs/images/tutorials/loop-string-hadron-dynamics/extracted-outputs/pauli-prop-heatmap-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "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": [
        "<span id=\"next-steps\" />\n",
        "\n",
        "## Próximos pasos\n",
        "\n",
        "Si este trabajo te ha parecido interesante, te recomendamos que eches un vistazo al siguiente material:\n",
        "\n",
        "<Admonition type=\"tip\" title=\"Recomendaciones\">\n",
        "  * [Documentación](/docs/guides/get-started-with-estimator) sobre la primitiva «Estimator» de Qiskit: para obtener más información sobre cómo configurar las opciones de mitigación de errores\n",
        "  * [Técnicas de mitigación y supresión de errores](/docs/guides/error-mitigation-and-suppression-techniques) : para conocer TREX, ZNE y otros métodos de mitigación\n",
        "  * [Qiskit Pauli Propagation (pauli-prop)](https://github.com/Qiskit/pauli-prop) : simulación clásica acelerada con Rust mediante la retropropagación de Pauli\n",
        "</Admonition>\n",
        "\n",
        "<span id=\"references\" />\n",
        "\n",
        "## Referencias\n",
        "\n",
        "\\[1] El artículo original: Ilčić, Majumdar, Mathew et al. «Observación de una dinámica hadrónica no abeliana robusta y coherente en procesadores cuánticos con ruido» [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
}