Saltar al contenido principal

Observation of robust and coherent non-Abelian hadron dynamics on noisy quantum processors

Uso estimado: 6 minutos en un procesador Heron (ibm_boston o equivalente) (NOTA: Esto es solo una estimación. Tu tiempo de ejecución puede variar.)

Resultados de aprendizaje​

  • Cómo se pueden reformular las teorías de calibre de red no abelianas (específicamente SU(2)) usando el marco Loop-String-Hadron (LSH) para una simulación cuántica eficiente

  • Cómo construir circuitos de evolución temporal trotterizados para un hamiltoniano aproximado de teoría de calibre SU(2) y mapearlos en qubits

  • Cómo ejecutar estos circuitos en hardware de IBM Quantum® usando el primitivo Estimator de Qiskit con mitigación de errores de lectura

Requisitos previos​

Contexto​

Motivación​

La cromodinámica cuántica (QCD), la teoría de calibre 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 QCD en red destacan en propiedades estáticas, pero no pueden simular la dinámica en tiempo real debido al problema del signo. Las computadoras cuánticas ofrecen una vía para sortear esta barrera al codificar los grados de libertad del campo de calibre directamente en qubits.

Este tutorial demuestra una simulación de este tipo: usar hardware de IBM Quantum para simular la propagación de hadrones en tiempo real en una teoría de calibre de red SU(2) en (1+1) dimensiones — la teoría de calibre no abeliana más simple y un paso hacia la QCD completa.

El hamiltoniano de Kogut-Susskind​

La teoría se formula en una red espacial 1D con fermiones escalonados (materia) en los sitios y campos de calibre SU(2) en los enlaces. Después de reescalar a forma adimensional, el hamiltoniano es:

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

donde HEH_E es la energía del campo cromoeléctrico, HMH_M es el término de masa escalonada, HIH_I es el término de interacción materia-calibre (hopping), μ=2mgx\mu = 2\frac{m}{g}\sqrt{x} codifica la masa del fermión, y x=1g2a2x = \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→∞N \to \infty y x→∞x \to \infty.

El marco Loop-String-Hadron (LSH)​

Un desafío clave es que el espacio de Hilbert del campo de calibre en cada enlace es de dimensión infinita. El marco Loop-String-Hadron (LSH) aborda esto reformulando la teoría en términos de variables invariantes de calibre — bucles de flujo, cuerdas que conectan cargas separadas, y hadrones (pares de fermiones singlete de calibre en un sitio). En la base LSH, la ley de Gauss se satisface automáticamente por construcción, por lo que todo estado base es físico. Cada sitio de la red se caracteriza por tres números cuánticos (nl,ni,no)(n_l, n_i, n_o) que representan el número de bucle, la cuerda entrante y la cuerda saliente, donde ni,no∈{0,1}n_i, n_o \in \{0,1\} son fermiónicos y nl≥0n_l \geq 0 es bosónico. El número de fermiones local se define a partir de estos como nf(r)=ni(r)+no(r)n_f(r) = n_i(r) + n_o(r) para sitios pares y nf(r)=2−[ni(r)+no(r)]n_f(r) = 2 - [n_i(r) + n_o(r)] para sitios impares.

Del hamiltoniano completo al circuito cuántico: tres aproximaciones clave​

El circuito cuántico no simula exactamente el hamiltoniano SU(2) completo. En su lugar, implementa una serie controlada de aproximaciones que son válidas en el régimen de acoplamiento débil (x≫1x \gg 1). Es esencial entender qué se aproxima y qué no:

Aproximación 1 — Límite de acoplamiento débil para HIH_I: El hamiltoniano de interacción completo HI(LSH)H_I^{\text{(LSH)}} (Ec. 16 en [1]) contiene prefactores que dependen del número cuántico bosónico nln_l a través de términos como 1/nl+11/\sqrt{n_l+1}. En el régimen de acoplamiento débil (x≫1x \gg 1), la dinámica está dominada por el término eléctrico HEH_E, que favorece estados con nln_l grande. Para nl≫1n_l \gg 1, la razón nl/(nl+1)→1n_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 hopping puramente local entre vecinos más cercanos:

HIapprox=−∑r[σ−(r)σ+(r+1)+σ+(r)σ−(r+1)],H_I^{\text{approx}} = -\sum_r \left[\sigma^-(r)\sigma^+(r+1) + \sigma^+(r)\sigma^-(r+1)\right],

que es independiente de nln_l y actúa solo sobre los qubits fermiónicos (ni,no)(n_i, n_o).

Aproximación 2 — Flujo promedio global para HEH_E: La energía eléctrica depende de nln_l en cada enlace. En el vacío de acoplamiento débil, nln_l es grande y aproximadamente uniforme. Reemplaza los valores de nln_l dependientes del sitio con un único promedio global nˉl\bar{n}_l, haciendo que HEH_E sea una fase diagonal proporcional a la configuración de fermiones en cada sitio:

HEapprox=NhE0+∑{r′}(nˉl2+34)H_E^{\text{approx}} = N h_E^0 + \sum_{\{r'\}} \left(\frac{\bar{n}_l}{2} + \frac{3}{4}\right)

donde {r′}\{r'\} suma sobre los sitios en la configuración fermiónica (ni=0,no=1)(n_i=0, n_o=1), y hE0h_E^0 es una fase global que puedes ignorar.

Aproximación 3 — Trotterización: El operador de evolución temporal para un paso de duración δτ\delta_\tau se descompone como:

e−iδτW≈e−im~HM e−iδτHEapprox e−icHIapproxe^{-i\delta_\tau W} \approx e^{-i\tilde{m} H_M} \, e^{-i\delta_\tau H_E^{\text{approx}}} \, e^{-ic H_I^{\text{approx}}}

donde c=δτxc = \delta_\tau x, m~=δτμ\tilde{m} = \delta_\tau \mu, y θ=−δτ(nˉl/2+3/4)\theta = -\delta_\tau(\bar{n}_l/2 + 3/4). Esta descomposición de Trotter de primer orden introduce un error que se anula cuando δτ→0\delta_\tau \to 0. Fijamos δτ=0.0015\delta_\tau = 0.0015 en todo momento.

El resultado de estas tres aproximaciones es que solo los dos qubits fermiónicos por sitio (ni,no)(n_i, n_o) son dinámicos — el grado de libertad bosónico nln_l se ha absorbido en parámetros efectivos. Esto produce un circuito compacto con 2N2N qubits para NN sitios de la red, donde cada paso de Trotter tiene una profundidad constante de compuertas de dos qubits (13 por paso).

Qué simula este tutorial​

El tutorial simula la propagación de hadrones: comenzando desde el vacío de acoplamiento fuerte (un estado producto), coloca un mesón en el centro de la red y evoluciona en el tiempo. El protocolo de medición diferencial — ejecutar el circuito con y sin el mesón central, y luego restar — aísla la señal coherente del hadrón tanto del ruido del hardware como de los efectos de frontera. El resultado es un patrón de cono de luz de oscilaciones de densidad de fermiones característico de un modo de respiración de un mesón confinado.

Requisitos​

Antes de comenzar este tutorial, instala lo siguiente:

  • Qiskit SDK v2.0 o posterior, con soporte de visualización

  • Qiskit Runtime v0.22 o posterior (pip install qiskit-ibm-runtime)

  • Paquete Pauli Propagation (pip install pauli-prop)

  • NumPy (pip install numpy)

  • Matplotlib (pip install matplotlib)

Configuración​

Comienza 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 de construcción de circuitos:

  1. pair_hamiltonian_circuit: Implementa la unitaria de dos qubits UIU_I para el hamiltoniano de interacción aproximado entre sitios vecinos. La descomposición de compuertas es: CNOT→H→Rz(−c)→CNOT→Rz(c)→CNOT→H→CNOT\text{CNOT} \to H \to R_z(-c) \to \text{CNOT} \to R_z(c) \to \text{CNOT} \to H \to \text{CNOT}.

  2. electric_hamiltonian_circuit: Implementa la unitaria de dos qubits UEU_E para la energía aproximada del campo eléctrico en cada sitio. La descomposición de compuertas es: X→Rz(θ/2)→CNOT→Rz(−θ/2)→CNOT→Rz(θ/2)→XX \to R_z(\theta/2) \to \text{CNOT} \to R_z(-\theta/2) \to \text{CNOT} \to R_z(\theta/2) \to X.

  3. construct_circuit: Ensambla el circuito trotterizado completo, disponiendo en capas los términos de interacción, eléctrico y de masa con compuertas SWAP para gestionar la conectividad de los qubits.

# Added by doQumentation — required packages for this notebook
!pip install -q matplotlib numpy pauli-prop qiskit qiskit-ibm-runtime
# Import libraries

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import TwoSlopeNorm
from qiskit.circuit import QuantumCircuit
from qiskit.quantum_info import SparsePauliOp
from typing import Optional

import warnings

warnings.filterwarnings("ignore")
def pair_hamiltonian_circuit(c: float) -> QuantumCircuit:
    """Two-qubit unitary for the approximate interaction Hamiltonian H_I.

    Implements exp(-i * c * H_I^approx) for one pair of neighboring sites,
    where c = delta_tau * x.
    """
    qc_temp = QuantumCircuit(2)
    qc_temp.cx(1, 0)
    qc_temp.h(1)
    qc_temp.rz(-c, 1)
    qc_temp.cx(0, 1)
    qc_temp.rz(c, 1)
    qc_temp.cx(0, 1)
    qc_temp.h(1)
    qc_temp.cx(1, 0)
    return qc_temp

def electric_hamiltonian_circuit(theta: float) -> QuantumCircuit:
    """Two-qubit unitary for the approximate electric field Hamiltonian H_E.

    Implements exp(-i * theta * H_E^approx) for one lattice site,
    where theta = -delta_tau * (n_bar_l / 2 + 3/4).
    """
    qc_temp = QuantumCircuit(2)
    qc_temp.x(0)
    qc_temp.rz(theta / 2, 0)
    qc_temp.cx(0, 1)
    qc_temp.rz(-theta / 2, 1)
    qc_temp.cx(0, 1)
    qc_temp.rz(theta / 2, 1)
    qc_temp.x(0)
    return qc_temp

def construct_circuit(
    num_lattice_point: int,
    num_trotter_steps: int,
    c: float,
    theta: float,
    m: float,
    theory: Optional[int] = 2,
    barriers: Optional[bool] = False,
    measurement: Optional[bool] = False,
    add_init_state: Optional[bool] = True,
    inverse_mid: Optional[bool] = False,
) -> QuantumCircuit:
    """Construct the full Trotterized time-evolution circuit.

    Builds a circuit implementing n Trotter steps of the approximate SU(2)
    LSH Hamiltonian evolution. The qubit layout uses a zigzag ordering:
    n_i(0), n_i(1), n_o(0), n_o(1), n_i(2), n_i(3), n_o(2), n_o(3), ...
    which minimizes the number of SWAP layers needed.

    Args:
        num_lattice_point: Number of lattice sites
        (num_qubits = 2 * num_lattice_point).
        num_trotter_steps: Number of Trotter steps.
        c: Interaction parameter (delta_tau * x).
        theta: Electric field phase parameter.
        m: Mass parameter (m_tilde = delta_tau * mu).
        theory: 1 for single chain, 2 for SU(2). Default 2.
        barriers: Insert barriers between Trotter layers for
        visualization.
        measurement: Append measurements at the end.
        add_init_state: Prepare the half-filled (strong-coupling vacuum)
        initial state.
        inverse_mid: Swap the central sites
        (for differential measurement protocol).
    """
    num_qubits = theory * num_lattice_point
    qc = QuantumCircuit(num_qubits)

    if num_trotter_steps <= 0:
        return qc

    # --- Initial state preparation ---
    if add_init_state:
        i = 1
        while i < num_lattice_point:
            for j in range(theory):
                qc.x(i + j * num_lattice_point)
            i = i + 2
        if inverse_mid:
            mid_lattice_qubits = [num_qubits // 2 - 1, num_qubits // 2]
            qc.x(mid_lattice_qubits)
    else:
        i = 1
        while i < num_qubits - 1:
            qc.swap(i, i + 1)
            i = i + 4

    # --- Trotter steps ---
    for step in range(num_trotter_steps):
        if barriers:
            qc.barrier()

        # First SWAP layer (skipped at step 0 — absorbed into initial state mapping)
        if step > 0:
            i = 1
            while i < num_qubits - 1:
                qc.swap(i, i + 1)
                i = i + 4

        # First layer of pair interactions
        j = 0
        while j < num_qubits - 2:
            circ = pair_hamiltonian_circuit(c)
            qc.compose(circ, [j, j + 1], inplace=True)
            j = j + 2
        if num_lattice_point % 2 == 0:
            circ = pair_hamiltonian_circuit(c)
            qc.compose(circ, [j, j + 1], inplace=True)

        # Second SWAP layer
        i = 1
        while i < num_qubits - 1:
            qc.swap(i, i + 1)
            i = i + theory

        # Second layer of pair interactions
        j = 2
        while j < num_qubits - 3:
            circ = pair_hamiltonian_circuit(c)
            qc.compose(circ, [j, j + 1], inplace=True)
            j = j + 2
        if num_lattice_point % 2 != 0:
            circ = pair_hamiltonian_circuit(c)
            qc.compose(circ, [j, j + 1], inplace=True)

        # Third SWAP layer
        i = 3
        while i < num_qubits - 1:
            qc.swap(i, i + 1)
            i = i + 2 * theory

        # Electric field term
        if theta != 0:
            e_circ = electric_hamiltonian_circuit(theta)
            for j in range(num_lattice_point):
                qc.compose(e_circ, [2 * j, 2 * j + 1], inplace=True)

        # Mass term: Rz(-m_tilde) for even sites, Rz(m_tilde) for odd sites
        for q in range(num_qubits):
            if q % 2 == 0:
                qc.rz(-1 * m, q)
            else:
                qc.rz(m, q)

    if measurement:
        qc.measure_all()

    return qc
def get_probabilities(expval: float):
    """Convert a Z-expectation value to site occupation probability.

    Since <Z> = p(0) - p(1), the occupation probability is p(1) = (1 - <Z>) / 2.
    """
    p1 = round((1 - expval) / 2, 3)
    return p1

def get_number(expval_data, num_lattice_point):
    """Convert raw Z-expectation values to staggered fermion number n_f at each site.

    n_f(r) = n_i(r) + n_o(r)           for even r
    n_f(r) = 2 - [n_i(r) + n_o(r)]     for odd r

    The two qubits per site encode (n_i, n_o), and occupation probabilities
    give us <n_i> and <n_o>.
    """
    N = []
    for expvals in expval_data:
        Pstep = [get_probabilities(expval) for expval in expvals]
        Nstep = []
        for k in range(num_lattice_point):
            val = Pstep[2 * k] + Pstep[2 * k + 1]
            a = 2 * (k % 2) + (1 - 2 * (k % 2)) * val
            Nstep.append(float(a))
        N.append(Nstep)
    return N

def calculate_difference(N, N_mid, num_lattice_point):
    """Differential measurement protocol: |n_f(meson) - n_f(vacuum)|.

    Subtracting the vacuum (SCV) evolution from the meson evolution
    isolates the coherent hadron signal from symmetric noise and boundary effects.
    """
    N_diff = []
    for i in range(len(N)):
        Nstep_diff = []
        for j in range(num_lattice_point):
            Nstep_diff.append(abs(N[i][j] - N_mid[i][j]))
        N_diff.append(Nstep_diff)
    return N_diff

Ejemplo con simulador a pequeña escala​

Primero, demuestra el flujo de trabajo a pequeña escala usando una red de seis sitios (12 qubits), para que puedas verificar la construcción del circuito y comprender los observables físicos antes de ejecutar en hardware.

Paso 1: Mapear entradas clásicas a un problema cuántico​

Define los parámetros físicos que coinciden con el régimen de acoplamiento débil estudiado en el artículo (x=100x = 100, m/g=1m/g = 1). Los parámetros del circuito derivados son:

  • c=δτ⋅x=0.15c = \delta_\tau \cdot x = 0.15 (parámetro de interacción)

  • θ=−δτ(nˉl/2+3/4)=0.01\theta = -\delta_\tau (\bar{n}_l/2 + 3/4) = 0.01 (fase del campo eléctrico)

  • m~=δτ⋅μ=0.03\tilde{m} = \delta_\tau \cdot \mu = 0.03 (parámetro de masa)

Para cada número de pasos de Trotter, construye 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.

# Physical / circuit parameters
num_lattice_point = 6  # 6 lattice sites -> 12 qubits for SU(2)
num_qubits = 2 * num_lattice_point
c = 0.15  # delta_tau * x
theta = 0.01  # electric field phase
m = 0.03  # m_tilde = delta_tau * mu
trotter_steps = range(1, 11)  # 10 Trotter steps

print(f"Lattice sites: {num_lattice_point}, Qubits: {num_qubits}")
print(f"Parameters: c={c}, theta={theta}, m_tilde={m}")
Lattice sites: 6, Qubits: 12
Parameters: c=0.15, theta=0.01, m_tilde=0.03
# Build circuits: meson initial state and vacuum (SCV) initial state
circuits_mid = [
    construct_circuit(
        num_lattice_point,
        d,
        c,
        theta,
        m,
        barriers=False,
        measurement=False,
        add_init_state=True,
        inverse_mid=True,
    )
    for d in trotter_steps
]

circuits = [
    construct_circuit(
        num_lattice_point,
        d,
        c,
        theta,
        m,
        barriers=False,
        measurement=False,
        add_init_state=True,
        inverse_mid=False,
    )
    for d in trotter_steps
]

# Visualize a single Trotter step
print(
    f"Circuit for 1 Trotter step: {circuits[0].num_qubits} qubits, depth {circuits[0].depth()}"
)
circuits[0].draw("mpl", fold=-1)
Circuit for 1 Trotter step: 12 qubits, depth 26

Output of the previous code cell

Paso 2: Optimizar el problema para la ejecución en hardware cuántico​

Define los observables: mediciones ZZ de un solo qubit en cada qubit. A partir de ⟨Z⟩\langle Z \rangle puedes extraer las probabilidades de ocupación y luego el número de fermiones escalonado nf(r)n_f(r) en cada sitio de la red rr.

# Z observable on each qubit
observables = [
    SparsePauliOp("I" * i + "Z" + "I" * (num_qubits - i - 1))
    for i in range(num_qubits)
]

print(f"Number of observables: {len(observables)}")
Number of observables: 12

Paso 3: Ejecutar usando primitivos de Qiskit​

Usa StatevectorEstimator para una simulación exacta sin ruido a pequeña escala.

from qiskit.primitives import StatevectorEstimator

estimator = StatevectorEstimator()

# Run meson circuits
pubs_mid = [(circuit, observables) for circuit in circuits_mid]
result_mid = estimator.run(pubs_mid).result()

# Run vacuum (SCV) circuits
pubs = [(circuit, observables) for circuit in circuits]
result = estimator.run(pubs).result()

# Extract expectation values
raw_expvals_mid = [
    result_mid[i].data.evs[::-1] for i in range(len(circuits_mid))
]
raw_expvals = [result[i].data.evs[::-1] for i in range(len(circuits))]

print(f"Computed expectation values for {len(raw_expvals)} Trotter steps")
Computed expectation values for 10 Trotter steps

Paso 4: Postprocesar y devolver el resultado en el formato clásico deseado​

Convierte los valores esperados al número de fermiones escalonado nf(r,t)n_f(r, t) y aplica el protocolo de medición diferencial (mesón −- vacío) para producir el mapa de calor de propagación del hadrón. Esto reproduce la estructura de la Figura 3 del artículo de referencia: sitio de la red rr en el eje x, paso de Trotter (tiempo) tt en el eje y, y nf(r,t)n_f(r,t) como escala de color.

# Compute fermion numbers
N_mid_sim = get_number(raw_expvals_mid, num_lattice_point)
N_sim = get_number(raw_expvals, num_lattice_point)
N_diff_sim = calculate_difference(N_mid_sim, N_sim, num_lattice_point)

# --- Reproduce Figure 3 style: Staggered Fermionic Occupation Number Dynamics ---
fig, axes = plt.subplots(1, 3, figsize=(18, 5))

# Convert to numpy arrays for plotting
N_mid_arr = np.array(N_mid_sim)
N_arr = np.array(N_sim)
N_diff_arr = np.array(N_diff_sim)

# Color scheme
vmax = max(max(sublist) for sublist in N_arr)
vmin = -vmax

# Meson evolution
norm1 = TwoSlopeNorm(vmin=vmin, vcenter=0, vmax=vmax)
im0 = axes[0].imshow(
    N_mid_arr,
    aspect="auto",
    origin="lower",
    cmap="RdBu_r",
    norm=norm1,
    extent=[0, num_lattice_point, 0.5, len(trotter_steps) + 0.5],
)
axes[0].set_xlabel("Lattice site $r$", fontsize=12)
axes[0].set_ylabel("Trotter step $t$", fontsize=12)
axes[0].set_title("$n_f(r,t)$ — Meson initial state", fontsize=12)
plt.colorbar(im0, ax=axes[0], label="$n_f(r,t)$")

# Vacuum (SCV) evolution
im1 = axes[1].imshow(
    N_arr,
    aspect="auto",
    origin="lower",
    cmap="RdBu_r",
    norm=norm1,
    extent=[0, num_lattice_point, 0.5, len(trotter_steps) + 0.5],
)
axes[1].set_xlabel("Lattice site $r$", fontsize=12)
axes[1].set_ylabel("Trotter step $t$", fontsize=12)
axes[1].set_title("$n_f(r,t)$ — Vacuum (SCV)", fontsize=12)
plt.colorbar(im1, ax=axes[1], label="$n_f(r,t)$")

# Differential: meson - vacuum
norm2 = TwoSlopeNorm(vmin=vmin, vcenter=0, vmax=vmax)
im2 = axes[2].imshow(
    N_diff_arr,
    aspect="auto",
    origin="lower",
    cmap="RdBu_r",
    norm=norm2,
    extent=[0, num_lattice_point, 0.5, len(trotter_steps) + 0.5],
)
axes[2].set_xlabel("Lattice site $r$", fontsize=12)
axes[2].set_ylabel("Trotter step $t$", fontsize=12)
axes[2].set_title(
    "Staggered Fermionic Occupation Number Dynamics\n$|n_f^{\\mathrm{meson}} - n_f^{\\mathrm{vacuum}}|$",
    fontsize=12,
)
plt.colorbar(im2, ax=axes[2], label="$n_f(r,t)$")

plt.suptitle(
    f"Hadron propagation — {num_lattice_point}-site lattice (StatevectorEstimator)",
    fontsize=14,
    y=1.02,
)
plt.tight_layout()
plt.show()

Output of the previous code cell

Ejemplo de hardware a gran escala​

Ahora escalamos a una red de 30 sitios (60 qubits) en hardware de IBM Quantum. A esta escala, el circuito con 10 pasos de Trotter comprende más de 3400 compuertas de dos qubits y 14,000 compuertas de un solo qubit.

Pasos 1-4 (comprimidos en un solo bloque de código)​

Aspectos clave del flujo de trabajo de hardware:

  • 10 pasos de Trotter para los circuitos del mesón y del vacío (intercalados para una deriva mínima)

  • Transpilación con optimization_level=1 — el diseño del circuito ya es isomorfo a la topología del dispositivo (una cadena lineal), por lo que no se necesitan SWAP de enrutamiento. El transpiler se usa únicamente para seleccionar una cadena de qubits físicos de bajo ruido y descomponer las compuertas en el conjunto de compuertas nativo.

  • EstimatorV2 con mitigación de errores de lectura TREX y Pauli twirling

  • Sesión Batch para enviar todos los trabajos juntos

# -------------------------Step 1: Define parameters & build circuits-------------------------

from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
from qiskit_ibm_runtime import EstimatorV2, Batch
from qiskit_ibm_runtime.options import (
    EstimatorOptions,
    ResilienceOptionsV2,
    TwirlingOptions,
    DynamicalDecouplingOptions,
)

service = QiskitRuntimeService()

num_lattice_point_hw = 30
num_qubits_hw = 2 * num_lattice_point_hw  # 60 qubits
c_hw = 0.15
theta_hw = 0.01
m_hw = 0.03
trotter_steps_hw = range(1, 11)  # 10 Trotter steps

# Build meson and vacuum circuits
circuits_mid_hw = [
    construct_circuit(
        num_lattice_point_hw,
        d,
        c_hw,
        theta_hw,
        m_hw,
        barriers=False,
        measurement=False,
        add_init_state=True,
        inverse_mid=True,
    )
    for d in trotter_steps_hw
]

circuits_hw = [
    construct_circuit(
        num_lattice_point_hw,
        d,
        c_hw,
        theta_hw,
        m_hw,
        barriers=False,
        measurement=False,
        add_init_state=True,
        inverse_mid=False,
    )
    for d in trotter_steps_hw
]

print(f"Built {len(circuits_hw)} circuit pairs for {num_qubits_hw} qubits")

# -------------------------Step 2: Transpile for hardware-------------------------
# The circuit topology is a linear chain, isomorphic to the device topology.
# We use optimization_level=1 since no routing SWAPs are needed — the transpiler
# only needs to select a low-noise qubit chain and decompose to native gates.

backend = service.backend("ibm_boston")

layout = [
    140,
    141,
    142,
    143,
    136,
    123,
    122,
    121,
    116,
    101,
    102,
    103,
    96,
    83,
    82,
    81,
    76,
    61,
    62,
    63,
    64,
    65,
    66,
    67,
    68,
    69,
    78,
    89,
    88,
    87,
    97,
    107,
    106,
    105,
    117,
    125,
    126,
    127,
    137,
    147,
    148,
    149,
    150,
    151,
    152,
    153,
    154,
    155,
    139,
    135,
    134,
    133,
    132,
    131,
    130,
    129,
    118,
    109,
    110,
    111,
]

pm = generate_preset_pass_manager(
    optimization_level=1, backend=backend, initial_layout=layout
)

isa_circuits_mid = pm.run(circuits_mid_hw)
isa_circuits = pm.run(circuits_hw)

print(f"Transpiled circuits. Example depth: {isa_circuits[0].depth()}")

# Define and layout-map observables
observables_hw = [
    SparsePauliOp("I" * i + "Z" + "I" * (num_qubits_hw - i - 1))
    for i in range(num_qubits_hw)
]

isa_observables_mid = [
    [obs.apply_layout(isa_circuits_mid[i].layout) for obs in observables_hw]
    for i in range(len(isa_circuits_mid))
]
isa_observables = [
    [obs.apply_layout(isa_circuits[i].layout) for obs in observables_hw]
    for i in range(len(isa_circuits))
]

# Build PUBs — interleave meson and vacuum for each Trotter step
isa_pubs_mid = [
    (circ, obs) for circ, obs in zip(isa_circuits_mid, isa_observables_mid)
]
isa_pubs = [(circ, obs) for circ, obs in zip(isa_circuits, isa_observables)]

pubs_to_execute = [
    [isa_pubs_mid[i], isa_pubs[i]] for i in range(len(isa_pubs))
]

# -------------------------Step 3: Execute on hardware-------------------------

twirling_options = TwirlingOptions(
    enable_gates=True,
    enable_measure=True,
    shots_per_randomization="auto",
    strategy="active-circuit",
)

resilience_options = ResilienceOptionsV2(
    measure_mitigation=True,  # TREX readout error mitigation
    zne_mitigation=False,  # ZNE turned off
)

dd_options = DynamicalDecouplingOptions(
    enable=False  # Circuit is sufficiently dense
)

options = EstimatorOptions(
    resilience=resilience_options,
    twirling=twirling_options,
    dynamical_decoupling=dd_options,
    default_shots=10_000,
)

ids = []
with Batch(backend=backend) as batch:
    for idx, pub in enumerate(pubs_to_execute):
        print(f"Submitting job for Trotter step {idx + 1}")
        estimator = EstimatorV2(mode=batch, options=options)
        estimator.skip_transpilation = True
        job = estimator.run(pub)
        ids.append(job.job_id())
    batch_id = batch.session_id

job_info = {"ids": ids, "batch_id": batch_id}
print(f"Submitted {len(ids)} jobs. Batch ID: {batch_id}")
print(ids)
# -------------------------Step 4: Post-process results-------------------------

jobs = [service.job(job_id) for job_id in ids]
results = [job.result() for job in jobs]

# Extract expectation values (index 0 = meson, index 1 = vacuum)
raw_expvals_mid_hw = [result[0].data.evs[::-1] for result in results]
raw_expvals_hw = [result[1].data.evs[::-1] for result in results]

# Compute fermion numbers and differential
N_mid_hw = get_number(raw_expvals_mid_hw, num_lattice_point_hw)
N_hw = get_number(raw_expvals_hw, num_lattice_point_hw)
N_diff_hw = calculate_difference(N_mid_hw, N_hw, num_lattice_point_hw)
N_diff_hw_arr = np.array(N_diff_hw)

fig, ax = plt.subplots(figsize=(10, 6))
vmax = np.abs(N_hw).max()
vmin = -vmax
norm = TwoSlopeNorm(vmin=vmin, vcenter=0, vmax=vmax)
im = ax.imshow(
    N_diff_hw_arr,
    aspect="auto",
    origin="lower",
    cmap="RdBu_r",
    norm=norm,
    extent=[0, 8, 0.5, len(trotter_steps_hw) + 0.5],
)
ax.set_xlabel("Lattice site $r$", fontsize=13)
ax.set_ylabel("Trotter step $t$", fontsize=13)
ax.set_title(
    "Staggered Fermionic Occupation Number Dynamics\nQuantum Simulation on IBM Hardware — 30-site lattice (60 qubits)",
    fontsize=13,
)
cbar = plt.colorbar(im, ax=ax)
cbar.set_label("$n_f(r,t)$", fontsize=12)
plt.tight_layout()
plt.show()

Output of the previous code cell

Evaluación comparativa clásica mediante Pauli Propagation​

El método Pauli Propagation (PPM) proporciona una simulación clásica sin ruido del circuito cuántico al retropropagar los observables medidos a través del circuito en la representación de Heisenberg. Bajo capas de Clifford (compuertas CNOT, H, S, X), los operadores de Pauli se mapean a otros operadores de Pauli sin aumentar el número de términos. Las capas no Clifford (las compuertas RzR_z en el circuito) pueden causar ramificación — en el peor de los casos, duplicando el número de términos — pero muchas ramas tienen coeficientes pequeños y se pueden truncar.

El flujo de trabajo con pauli-prop es:

  1. Divide el circuito en sus partes Clifford y no Clifford usando evolve_through_cliffords.

  2. Propaga cada observable a través de la parte no Clifford usando propagate_through_circuit, manteniendo hasta max_terms términos de Pauli y descartando términos con coeficientes por debajo del umbral de truncamiento atol.

  3. Evoluciona el resultado a través de la parte Clifford usando el soporte de Clifford integrado de Qiskit.

  4. Extrae el valor esperado sumando los coeficientes de los términos de Pauli diagonales (que contienen solo II y ZZ).

Umbral de truncamiento​

El parámetro atol en propagate_through_circuit controla con qué agresividad se podan las ramas pequeñas de Pauli. Un umbral muy estricto (por ejemplo, 1e-12) conserva casi todas las ramas y da resultados exactos, pero el tiempo de simulación crece marcadamente con la profundidad del circuito; la simulación de 120 qubits en el artículo tardó aproximadamente 8.5 horas con la configuración predeterminada. Aumentar el umbral (por ejemplo, a 1e-6 o 1e-3) descarta términos cuyos coeficientes caen por debajo de ese valor, reduciendo drásticamente el número de términos rastreados y acelerando el cálculo. La contrapartida es un error de aproximación pequeño y controlable que puedes validar comparando los resultados con diferentes umbrales.

import time
from pauli_prop import evolve_through_cliffords, propagate_through_circuit

# ── PPM Configuration ──
# Truncation threshold: controls the speed/accuracy trade-off.
PPM_THRESHOLD = 1e-3

# Maximum Pauli terms to track per observable (hard cap on memory/time)
PPM_MAX_TERMS = 66_000

print(f"PPM settings: atol={PPM_THRESHOLD}, max_terms={PPM_MAX_TERMS}")

# We propagate each single-qubit Z observable through each circuit.
# For PPM, we work with the un-transpiled circuits (ideal noiseless simulation).

observables_pp = [
    SparsePauliOp("I" * i + "Z" + "I" * (num_qubits_hw - i - 1))
    for i in range(num_qubits_hw)
]

def ppm_expectation_values(
    circuit, observables, max_terms=PPM_MAX_TERMS, atol=PPM_THRESHOLD
):
    """Compute expectation values of single-qubit Z observables
    via Pauli propagation.

    Args:
        circuit: The quantum circuit to simulate.
        observables: List of single-qubit Z observables.
        max_terms: Maximum number of Pauli terms to retain (hard cap).
        atol: Absolute tolerance — Pauli terms with coefficients below this
              value are discarded during propagation. Larger values give
              faster simulation at the cost of approximation accuracy.
    """
    circuit = circuit.decompose(["swap"])  # decompose SWAPs into 3 CX gates
    cliff, non_cliff = evolve_through_cliffords(circuit)

    evs = []
    for obs in observables:
        evolved_obs = propagate_through_circuit(
            obs, non_cliff, max_terms=max_terms, atol=atol, frame="h"
        )[0]
        evolved_obs.paulis = evolved_obs.paulis.evolve(cliff, frame="h")
        diagonal_mask = ~evolved_obs.paulis.x.any(axis=1)
        ev = float(evolved_obs.coeffs[diagonal_mask].sum().real)
        evs.append(ev)
    return np.array(evs)

# Run PPM for each Trotter step and record wall-clock time
pp_expvals_mid = []
pp_expvals = []
pp_times = []

for idx, d in enumerate(trotter_steps_hw):
    t_start = time.perf_counter()

    # Meson circuit
    evs_mid = ppm_expectation_values(circuits_mid_hw[idx], observables_pp)

    # Vacuum circuit
    evs_vac = ppm_expectation_values(circuits_hw[idx], observables_pp)

    elapsed = time.perf_counter() - t_start
    pp_times.append(elapsed)

    pp_expvals_mid.append(evs_mid[::-1])
    pp_expvals.append(evs_vac[::-1])

    print(f"Trotter step {d:2d}: {elapsed:.1f} s")

print(f"\nTotal PPM simulation time: {sum(pp_times):.1f} s")
print(f"Truncation threshold used: {PPM_THRESHOLD}")
PPM settings: atol=0.001, max_terms=66000
Trotter step  1: 5.0 s
Trotter step  2: 7.5 s
Trotter step  3: 11.2 s
Trotter step  4: 14.7 s
Trotter step  5: 18.3 s
Trotter step  6: 22.1 s
Trotter step  7: 25.6 s
Trotter step  8: 29.4 s
Trotter step  9: 33.2 s
Trotter step 10: 36.6 s

Total PPM simulation time: 203.6 s
Truncation threshold used: 0.001
# --- PPM simulation time vs. Trotter steps ---
fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(
    list(trotter_steps_hw),
    pp_times,
    "o-",
    color="tab:blue",
    linewidth=2,
    markersize=6,
)
ax.set_xlabel("Trotter step", fontsize=13)
ax.set_ylabel("Wall-clock time (s)", fontsize=13)
ax.set_title(
    "Pauli Propagation simulation time vs. Trotter steps\n(30-site lattice, 60 qubits)",
    fontsize=13,
)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

Output of the previous code cell

# --- PPM heatmap and comparison with hardware ---
N_mid_pp = get_number(pp_expvals_mid, num_lattice_point_hw)
N_pp = get_number(pp_expvals, num_lattice_point_hw)
N_diff_pp = calculate_difference(N_mid_pp, N_pp, num_lattice_point_hw)

N_diff_pp_arr = np.array(N_diff_pp)

fig, axes = plt.subplots(1, 2, figsize=(18, 6))
vmax = np.abs(N_hw).max()
vmin = -vmax
norm = TwoSlopeNorm(vmin=vmin, vcenter=0, vmax=vmax)

# PPM result
im0 = axes[0].imshow(
    N_diff_pp_arr,
    aspect="auto",
    origin="lower",
    cmap="RdBu_r",
    norm=norm,
    extent=[0, num_lattice_point_hw, 0.5, len(trotter_steps_hw) + 0.5],
)
axes[0].set_xlabel("Lattice site $r$", fontsize=12)
axes[0].set_ylabel("Trotter step $t$", fontsize=12)
axes[0].set_title(
    "Pauli Propagation\n(classical noiseless simulation)", fontsize=12
)
plt.colorbar(im0, ax=axes[0], label="$n_f(r,t)$")

# Hardware result
im1 = axes[1].imshow(
    N_diff_hw_arr,
    aspect="auto",
    origin="lower",
    cmap="RdBu_r",
    norm=norm,
    extent=[0, num_lattice_point_hw, 0.5, len(trotter_steps_hw) + 0.5],
)
axes[1].set_xlabel("Lattice site $r$", fontsize=12)
axes[1].set_ylabel("Trotter step $t$", fontsize=12)
axes[1].set_title(
    "Quantum Simulation\n(IBM Hardware, readout error mitigation only)",
    fontsize=12,
)
plt.colorbar(im1, ax=axes[1], label="$n_f(r,t)$")

plt.suptitle(
    "Staggered Fermionic Occupation Number Dynamics — 30-site lattice",
    fontsize=14,
    y=1.02,
)
plt.tight_layout()
plt.show()

Output of the previous code cell

Próximos pasos​

Si te ha parecido interesante este trabajo, considera explorar el siguiente material:

Recomendaciones

Referencias​

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

Fuente: IBM Quantum docs — actualizado el 9 sept 2026
Versión en inglés en doQumentation — actualizado el 9 sept 2026
¿Problema con el sitio o la traducción?¿Problema de contenido? Editar en IBM Quantum docs