Skip to main content
IBM Quantum Platform

Observación de una dinámica hadrónica no abeliana robusta y coherente en procesadores cuánticos con ruido

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.)

Resultados del aprendizaje

Al finalizar este tutorial, habrás aprendido lo siguiente:

  • 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
  • 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
  • Cómo ejecutar estos circuitos en un hardware de tip IBM Quantum® utilizando la primitiva «Qiskit Estimator» con mitigación de errores de lectura

Requisitos previos

Se recomienda que te familiarices con los siguientes temas:

En segundo plano

Motivació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.

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—.

El hamiltoniano de Kogut-Susskind

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

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 (salto), μ=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)

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 (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. A partir de estos, el número de fermiones local se define como « nf(r)=ni(r)+no(r)n_f(r) = n_i(r) + n_o(r) » para los sitios pares y « nf(r)=2−[ni(r)+no(r)]n_f(r) = 2 - [n_i(r) + n_o(r)] » para los sitios impares.

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

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≫1x \gg 1 ). Es fundamental comprender qué se aproxima y qué no:

Aproximación 1 — Límite de acoplamiento débil para « HIH_I »: El hamiltoniano de interacción completa HI(LSH)H_I^{\text{(LSH)}} (ecuación 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 viene dominada por el término eléctrico HEH_E, que favorece los estados con un nln_l elevado. Para nl≫1n_l \gg 1, la relació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 salto entre vecinos más cercanos de carácter puramente local:

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 únicamente sobre los qubits fermiónicos (ni,no)(n_i, n_o).

Aproximación 2 — Flujo medio 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. Sustituye los valores de nln_l, que dependen de cada sitio, por un único valor medio global nˉl\bar{n}_l, de modo 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'\} se suma por todos 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 de la siguiente manera:

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 integrado en los parámetros efectivos. Esto da como resultado un circuito compacto con 2N2N qubits para NN nodos de la red, en el que cada paso de Trotter tiene una profundidad constante de puertas de dos qubits (13 por paso).

Qué simula este tutorial

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.

Requisitos

Antes de empezar este tutorial, instala lo siguiente:

  • Qiskit SDK v2.0 o posterior, con soporte para visualización
  • Qiskit Runtime v0.22 o posterior (pip install qiskit-ibm-runtime)
  • Paquete de propagación de Pauli (pip install pauli-prop)
  • NumPy (pip install numpy)
  • Matplotlib (pip install matplotlib)

Configuració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:

  1. pair_hamiltonian_circuit: Implementa el operador unitario de dos qubits UIU_I para el hamiltoniano de interacción aproximado entre sitios vecinos. La descomposición de la puerta 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 el operador unitario de dos qubits UEU_E para la energía aproximada del campo eléctrico en cada sitio. La descomposición de la puerta 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: 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.

# 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 de simulador a pequeña escala

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.

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

Defina los parámetros físicos que se ajustan al régimen de acoplamiento débil estudiado en el artículo ( x=100x = 100, m/g=1m/g = 1 ). Los parámetros del circuito obtenidos 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 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.

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

print(f"Lattice sites: {num_lattice_point}, Qubits: {num_qubits}")
print(f"Parameters: c={c}, theta={theta}, m_tilde={m}")

Output:

Lattice sites: 6, Qubits: 12
Parameters: c=0.15, theta=0.01, m_tilde=0.03

# Build circuits: meson initial state and vacuum (SCV) initial state
circuits_mid = [
    construct_circuit(
        num_lattice_point,
        d,
        c,
        theta,
        m,
        barriers=False,
        measurement=False,
        add_init_state=True,
        inverse_mid=True,
    )
    for d in trotter_steps
]

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

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

Output:

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

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

Definir las magnitudes observables: mediciones de « ZZ » de un solo qubit en cada qubit. En ⟨Z⟩\langle Z \rangle se pueden extraer las probabilidades de ocupación y, a continuación, 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)}")

Output:

Number of observables: 12

Paso 3: Ejecutar utilizando Qiskit primitives

Utilízalo StatevectorEstimator para realizar simulaciones exactas y 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")

Output:

Computed expectation values for 10 Trotter steps

Paso 4: Realizar el posprocesamiento y obtener el resultado en el formato clásico deseado

Convertir los valores esperados al número de fermiones escalonado nf(r,t)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 rr en el eje x, el paso de Trotter (tiempo) tt en el eje y, y nf(r,t)n_f(r,t) como escala de colores.

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

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

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

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

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

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

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

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

Output:

Output of the previous code cell

Ejemplo de hardware a gran escala

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.

Pasos 1-4 (agrupados en un único 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 minimizar la deriva)
  • 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.
  • EstimatorV2 con mitigación de errores de lectura de TREX y giro de Pauli
  • Batch sesión para enviar todos los trabajos a la vez
# -------------------------Step 1: Define parameters & build circuits-------------------------

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

service = QiskitRuntimeService()

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

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

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

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

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

backend = service.backend("ibm_boston")

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Output:

Output of the previous code cell

Evaluación comparativa clásica mediante la propagación de Pauli

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 « RzR_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.

El proceso con pauli-prop es el siguiente:

  1. Divide el circuito en sus partes de Clifford y no Clifford utilizando evolve_through_cliffords.
  2. atolPropaga 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.
  3. Calcula el resultado mediante la operación de Clifford utilizando la compatibilidad integrada de Qiskit con las operaciones de Clifford.
  4. Se obtiene el valor esperado sumando los coeficientes de los términos diagonales de Pauli (que contienen únicamente II y ZZ ).

Umbral de truncamiento

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 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.

import time
from pauli_prop import evolve_through_cliffords, propagate_through_circuit

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

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

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

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

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


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

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

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


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

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

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

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

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

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

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

print(f"\nTotal PPM simulation time: {sum(pp_times):.1f} s")
print(f"Truncation threshold used: {PPM_THRESHOLD}")

Output:

PPM settings: atol=0.001, max_terms=66000
Trotter step  1: 5.0 s
Trotter step  2: 7.5 s
Trotter step  3: 11.2 s
Trotter step  4: 14.7 s
Trotter step  5: 18.3 s
Trotter step  6: 22.1 s
Trotter step  7: 25.6 s
Trotter step  8: 29.4 s
Trotter step  9: 33.2 s
Trotter step 10: 36.6 s

Total PPM simulation time: 203.6 s
Truncation threshold used: 0.001

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

Output:

Output of the previous code cell 
# --- PPM heatmap and comparison with hardware ---
N_mid_pp = get_number(pp_expvals_mid, num_lattice_point_hw)
N_pp = get_number(pp_expvals, num_lattice_point_hw)
N_diff_pp = calculate_difference(N_mid_pp, N_pp, num_lattice_point_hw)

N_diff_pp_arr = np.array(N_diff_pp)

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

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

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

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

Output:

Output of the previous code cell

Próximos pasos

Si este trabajo te ha parecido interesante, te recomendamos que eches un vistazo al siguiente material:

Recomendaciones

Referencias

[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 (2026)

¿Le ha resultado útil esta página?
Informe de un error, de una errata o solicite contenido en GitHub.