Pular para o conteúdo principal

Observação de dinâmica de hádrons não abeliana robusta e coerente em processadores quânticos ruidosos

Estimativa de uso: 6 minutos em um processador Heron (ibm_boston ou equivalente) (NOTA: Isto é apenas uma estimativa. Seu tempo de execução pode variar.)

Resultados de aprendizado​

  • Como teorias de calibre de rede não abelianas (especificamente SU(2)) podem ser reformuladas usando o framework Loop-String-Hadron (LSH) para simulação quântica eficiente

  • Como construir circuitos de evolução temporal trotterizados para um hamiltoniano aproximado da teoria de calibre SU(2) e mapeá-los em qubits

  • Como executar esses circuitos em hardware do IBM Quantum® usando a primitiva Qiskit Estimator com mitigação de erro de leitura

Pré-requisitos​

Contexto​

Motivação​

A Cromodinâmica Quântica (QCD), a teoria de calibre SU(3) da força forte, liga quarks em hádrons e governa o confinamento e a quebra de cordas. Métodos clássicos de QCD em rede são excelentes para propriedades estáticas, mas não conseguem simular dinâmicas em tempo real devido ao problema do sinal. Computadores quânticos oferecem uma forma de contornar essa barreira codificando os graus de liberdade do campo de calibre diretamente em qubits.

Este tutorial demonstra tal simulação: usar hardware do IBM Quantum para simular a propagação de hádrons em tempo real em uma teoria de calibre de rede SU(2) (1+1)-dimensional — a teoria de calibre não abeliana mais simples e um passo em direção à QCD completa.

O hamiltoniano de Kogut-Susskind​

A teoria é formulada em uma rede espacial 1D com férmions escalonados (matéria) nos sítios e campos de calibre SU(2) nos elos. Após reescalonamento para a forma adimensional, o hamiltoniano é:

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

onde HEH_E é a energia do campo cromoelétrico, HMH_M é o termo de massa escalonada, HIH_I é o termo de interação matéria-calibre (hopping), μ=2mgx\mu = 2\frac{m}{g}\sqrt{x} codifica a massa do férmion, e x=1g2a2x = \frac{1}{g^2 a^2} é a força de interação. O limite contínuo da teoria situa-se em N→∞N \to \infty e x→∞x \to \infty.

O framework Loop-String-Hadron (LSH)​

Um desafio importante é que o espaço de Hilbert do campo de calibre em cada elo é de dimensão infinita. O framework Loop-String-Hadron (LSH) aborda isso reformulando a teoria em termos de variáveis invariantes de calibre — laços de fluxo, cordas conectando cargas separadas, e hádrons (pares de férmions singletos de calibre em um sítio). Na base LSH, a lei de Gauss é satisfeita automaticamente por construção, então todo estado de base é físico. Cada sítio da rede é caracterizado por três números quânticos (nl,ni,no)(n_l, n_i, n_o) representando o número de laço, a corda de entrada, e a corda de saída, onde ni,no∈{0,1}n_i, n_o \in \{0,1\} são fermiônicos e nl≥0n_l \geq 0 é bosônico. O número fermiônico local é definido a partir destes como nf(r)=ni(r)+no(r)n_f(r) = n_i(r) + n_o(r) para sítios pares e nf(r)=2−[ni(r)+no(r)]n_f(r) = 2 - [n_i(r) + n_o(r)] para sítios ímpares.

Do hamiltoniano completo ao circuito quântico: três aproximações chave​

O circuito quântico não simula o hamiltoniano SU(2) completo de forma exata. Em vez disso, implementa uma série controlada de aproximações que são válidas no regime de acoplamento fraco (x≫1x \gg 1). Entender o que é e o que não é aproximado é essencial:

Aproximação 1 — Limite de acoplamento fraco para HIH_I: O hamiltoniano de interação completo HI(LSH)H_I^{\text{(LSH)}} (Eq. 16 em [1]) contém prefatores que dependem do número quântico bosônico nln_l por meio de termos como 1/nl+11/\sqrt{n_l+1}. No regime de acoplamento fraco (x≫1x \gg 1), a dinâmica é dominada pelo termo elétrico HEH_E, que favorece estados com nln_l grande. Para nl≫1n_l \gg 1, a razão nl/(nl+1)→1n_l/(n_l+1) \to 1 e todos esses prefatores se simplificam para a unidade. O hamiltoniano de interação então se reduz a um hopping puramente local entre vizinhos próximos:

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 é independente de nln_l e atua apenas nos qubits fermiônicos (ni,no)(n_i, n_o).

Aproximação 2 — Fluxo médio global para HEH_E: A energia elétrica depende de nln_l em cada elo. No vácuo de acoplamento fraco, nln_l é grande e aproximadamente uniforme. Substitua os valores de nln_l dependentes do sítio por uma única média global nˉl\bar{n}_l, tornando HEH_E uma fase diagonal proporcional à configuração fermiônica em cada sítio:

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)

onde {r′}\{r'\} soma sobre os sítios na configuração fermiônica (ni=0,no=1)(n_i=0, n_o=1), e hE0h_E^0 é uma fase global que você pode ignorar.

Aproximação 3 — Trotterização: O operador de evolução temporal para um passo de duração δτ\delta_\tau é decomposto 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}}}

onde c=δτxc = \delta_\tau x, m~=δτμ\tilde{m} = \delta_\tau \mu, e θ=−δτ(nˉl/2+3/4)\theta = -\delta_\tau(\bar{n}_l/2 + 3/4). Esta decomposição de Trotter de primeira ordem introduz um erro que se anula quando δτ→0\delta_\tau \to 0. Fixamos δτ=0.0015\delta_\tau = 0.0015 ao longo de todo o processo.

O resultado dessas três aproximações é que apenas os dois qubits fermiônicos por sítio (ni,no)(n_i, n_o) são dinâmicos — o grau de liberdade bosônico nln_l foi absorvido em parâmetros efetivos. Isso resulta em um circuito compacto com 2N2N qubits para NN sítios de rede, onde cada passo de Trotter tem profundidade constante de portas de dois qubits (13 por passo).

O que este tutorial simula​

O tutorial simula a propagação de hádrons: começando a partir do vácuo de acoplamento forte (um estado produto), coloque um méson no centro da rede e evolua no tempo. O protocolo de medição diferencial — executando o circuito com e sem o méson central, e então subtraindo — isola o sinal coerente do hádron tanto do ruído do hardware quanto dos efeitos de contorno. O resultado é um padrão de cone de luz de oscilações de densidade de férmions característico de um modo de respiração de méson confinado.

Requisitos​

Antes de iniciar este tutorial, instale o seguinte:

  • Qiskit SDK v2.0 ou posterior, com suporte a visualização

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

  • Pacote Pauli Propagation (pip install pauli-prop)

  • NumPy (pip install numpy)

  • Matplotlib (pip install matplotlib)

Configuração​

Comece importando as bibliotecas necessárias e definindo as funções auxiliares que constroem os circuitos quânticos para a evolução temporal LSH. Há três funções principais de construção de circuitos:

  1. pair_hamiltonian_circuit: Implementa a unitária de dois qubits UIU_I para o hamiltoniano de interação aproximado entre sítios vizinhos. A decomposição de portas é: 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 a unitária de dois qubits UEU_E para a energia aproximada do campo elétrico em cada sítio. A decomposição de portas é: 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 o circuito trotterizado completo, sobrepondo os termos de interação, elétrico e de massa com portas SWAP para gerenciar a conectividade dos 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

Exemplo de simulador em pequena escala​

Primeiro, demonstre o fluxo de trabalho em pequena escala usando uma rede de seis sítios (12 qubits), para que você possa verificar a construção do circuito e entender os observáveis físicos antes de executar no hardware.

Passo 1: Mapear entradas clássicas para um problema quântico​

Defina os parâmetros físicos correspondentes ao regime de acoplamento fraco estudado no artigo (x=100x = 100, m/g=1m/g = 1). Os parâmetros de circuito derivados são:

  • c=δτ⋅x=0.15c = \delta_\tau \cdot x = 0.15 (parâmetro de interação)

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

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

Para cada contagem de passos de Trotter, construa dois circuitos: um inicializando um méson no centro (inverse_mid=True) e outro preparando o vácuo de acoplamento forte (inverse_mid=False). O protocolo de medição diferencial subtrai a evolução do vácuo para isolar o sinal do hádron.

# 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

Passo 2: Otimizar o problema para execução em hardware quântico​

Defina os observáveis: medições ZZ de um único qubit em cada qubit. A partir de ⟨Z⟩\langle Z \rangle você pode extrair probabilidades de ocupação e depois o número de férmions escalonado nf(r)n_f(r) em cada sítio de rede 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

Passo 3: Executar usando primitivas Qiskit​

Use StatevectorEstimator para simulação exata sem ruído em pequena 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

Passo 4: Pós-processar e retornar o resultado no formato clássico desejado​

Converta os valores esperados no número de férmions escalonado nf(r,t)n_f(r, t) e aplique o protocolo de medição diferencial (méson −- vácuo) para produzir o mapa de calor de propagação do hádron. Isso reproduz a estrutura da Figura 3 do artigo de referência: sítio de rede rr no eixo x, passo de Trotter (tempo) tt no eixo y, e nf(r,t)n_f(r,t) como a escala de cores.

# 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

Exemplo de hardware em larga escala​

Agora escalamos para uma rede de 30 sítios (60 qubits) em hardware do IBM Quantum. Nessa escala, o circuito com 10 passos de Trotter compreende mais de 3400 portas de dois qubits e 14.000 portas de um qubit.

Passos 1-4 (comprimidos em um único bloco de código)​

Aspectos chave do fluxo de trabalho em hardware:

  • 10 passos de Trotter para os circuitos de méson e vácuo (intercalados para desvio mínimo)

  • Transpilação com optimization_level=1 — o layout do circuito já é isomorfo à topologia do dispositivo (uma cadeia linear), então nenhum SWAP de roteamento é necessário. O transpilador é usado apenas para selecionar uma cadeia de qubits físicos com baixo ruído e decompor as portas no conjunto de portas nativo.

  • EstimatorV2 com mitigação de erro de leitura TREX e Pauli twirling

  • Sessão Batch para submeter todos os jobs 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

Benchmarking clássico via Pauli Propagation​

O Pauli Propagation Method (PPM) fornece uma simulação clássica sem ruído do circuito quântico ao retropropagar os observáveis medidos através do circuito na imagem de Heisenberg. Sob camadas de Clifford (portas CNOT, H, S, X), os operadores de Pauli são mapeados em outros operadores de Pauli sem aumentar o número de termos. Camadas não-Clifford (as portas RzR_z no circuito) podem causar ramificação — no pior caso, dobrando o número de termos — mas muitos ramos têm coeficientes pequenos e podem ser truncados.

O fluxo de trabalho com pauli-prop é:

  1. Dividir o circuito em suas partes de Clifford e não-Clifford usando evolve_through_cliffords.

  2. Propagar cada observável através da parte não-Clifford usando propagate_through_circuit, mantendo até max_terms termos de Pauli e descartando termos com coeficientes abaixo do limite de truncamento atol.

  3. Evoluir o resultado através da parte de Clifford usando o suporte integrado do Qiskit para Clifford.

  4. Extrair o valor esperado somando os coeficientes dos termos de Pauli diagonais (contendo apenas II e ZZ).

Limite de truncamento​

O parâmetro atol em propagate_through_circuit controla o quão agressivamente pequenos ramos de Pauli são podados. Um limite muito rígido (por exemplo, 1e-12) mantém quase todos os ramos e fornece resultados exatos, mas o tempo de simulação cresce acentuadamente com a profundidade do circuito; a simulação de 120 qubits no artigo levou aproximadamente 8,5 horas com as configurações padrão. Aumentar o limite (por exemplo, para 1e-6 ou 1e-3) descarta termos cujos coeficientes ficam abaixo desse valor, reduzindo drasticamente o número de termos rastreados e acelerando o cálculo. O trade-off é um erro de aproximação pequeno e controlável, que você pode validar comparando resultados em diferentes limites.

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 passos​

Se você achou este trabalho interessante, considere explorar o seguinte material:

Recomendações

Referências​

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

Fonte: IBM Quantum docs — atualizado em 9 de set. de 2026
Versão em inglês no doQumentation — atualizada em 9 de set. de 2026
Problema com o site ou a tradução?Problema de conteúdo? Editar no IBM Quantum docs