Skip to main content
IBM Quantum Platform

Observação de dinâmica de hádrons não-abelianos robusta e coerente em processadores quânticos sujeitos a ruído

Estimativa de tempo de execução: 6 minutos em um processador Heron (ibm_boston ou equivalente) (NOTA: Trata-se apenas de uma estimativa. (O tempo de execução pode variar.)

Resultados do aprendizado

Ao concluir este tutorial, você terá aprendido o seguinte:

  • Como as teorias de gauge em reticulados não abelianos (especificamente SU(2)) podem ser reformuladas utilizando a estrutura Loop-String-Hadron (LSH) para uma simulação quântica eficiente
  • Como construir circuitos de evolução temporal de Trotter para um hamiltoniano aproximado de uma teoria de gauge SU(2) e mapeá-los para qubits
  • Como executar esses circuitos em um hardware d IBM Quantum®, utilizando a primitiva Qiskit Estimator com mitigação de erros de leitura

Pré-requisitos

Recomenda-se que você se familiarize com estes tópicos:

Segundo plano

Motivação

A Cromodinâmica Quântica (QCD), a teoria de gauge SU(3) da força forte, agrupa os quarks em hádrons e rege o confinamento e a quebra de cordas. Os métodos clássicos de QCD em rede se destacam no estudo de propriedades estáticas, mas não conseguem simular a dinâmica em tempo real devido ao problema do sinal. Os computadores quânticos oferecem uma maneira de contornar essa barreira, codificando os graus de liberdade dos campos de calibre diretamente nos qubits.

Este tutorial demonstra uma simulação desse tipo: o uso do hardwar IBM Quantum para simular a propagação de hádrons em tempo real em uma teoria de gauge em rede SU(2) de dimensão (1+1) — a teoria de gauge não-abeliana mais simples e um passo importante rumo à QCD completa.

O hamiltoniano de Kogut-Susskind

A teoria é formulada em uma rede espacial de tipo “ 1D ”, com férmions (matéria) dispostos de forma escalonada nos vértices e campos de calibre SU(2) nos arcos. Após a reescalonagem 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 (salto), μ=2mgx\mu = 2\frac{m}{g}\sqrt{x} representa a massa do férmion e x=1g2a2x = \frac{1}{g^2 a^2} é a intensidade da interação. O limite contínuo da teoria está em N→∞N \to \infty e x→∞x \to \infty.

A estrutura Loop-String-Hadron (LSH)

Um dos principais desafios é que o espaço de Hilbert do campo de calibre em cada elo é de dimensão infinita. A estrutura Loop-String-Hadron (LSH) aborda essa questão reformulando a teoria em termos de variáveis invariantes de calibre — laços de fluxo, cordas que conectam cargas separadas e hádrons (pares de férmions singletos de calibre em um nó). Na base LSH, a lei de Gauss é satisfeita automaticamente por definição; portanto, todo estado da base é físico. Cada vértice da rede é caracterizado por três números quânticos (nl,ni,no)(n_l, n_i, n_o), que representam o número de laço, a corda de entrada e a corda de saída, sendo que ni,no∈{0,1}n_i, n_o \in \{0,1\} são fermiónicos e nl≥0n_l \geq 0 é bosônico. A partir disso, o número de férmions local é definido 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 fundamentais

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

Aproximação 1 — Limite de acoplamento fraco para o modelo “ HIH_I ”: O hamiltoniano de interação completa HI(LSH)H_I^{\text{(LSH)}} (Eq. 16 em [1] ) contém pré-fatores 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 grande nln_l. Para nl≫1n_l \gg 1, a razão nl/(nl+1)→1n_l/(n_l+1) \to 1 e todos esses pré-fatores se simplificam para a unidade. O hamiltoniano de interação reduz-se, então, a um salto entre vizinhos mais próximos de natureza 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 é 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, a densidade de fluxo de energia ( nln_l ) é grande e aproximadamente uniforme. Substitua os valores de nln_l, que dependem do site, por um único valor médio global nˉl\bar{n}_l, tornando HEH_E uma fase diagonal proporcional à configuração dos férmions em cada site:

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'\} representa a 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 pode ser ignorada.

Aproximação 3 — Trotterização: O operador de evolução temporal para um passo de duração δτ\delta_\tau é decomposto da seguinte forma:

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). Essa decomposição de Trotter de primeira ordem introduz um erro que se anula à medida que δτ→0\delta_\tau \to 0. Fixamos δτ=0.0015\delta_\tau = 0.0015 em toda a análise.

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 pelos parâmetros efetivos. Isso resulta em um circuito compacto com 2N2N qubits para NN pontos da rede, em que cada etapa de Trotter possui profundidade constante de porta de dois qubits (13 por etapa).

O que este tutorial simula

O tutorial simula a propagação de hádrons : partindo do vácuo de acoplamento forte (um estado de produto), coloque um méson no centro da rede e faça a evolução no tempo. O protocolo de medição diferencial — que consiste em operar o circuito com e sem o méson central e, em seguida, subtrair os resultados — isola o sinal coerente do hádrão tanto do ruído do hardware quanto dos efeitos de contorno. O resultado é um padrão em forma de cone de luz de oscilações na 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 versão posterior, com suporte à visualização
  • Qiskit Runtime v0.22 ou posterior (pip install qiskit-ibm-runtime)
  • Pacote de propagação de Pauli (pip install pauli-prop)
  • NumPy (pip install numpy)
  • Matplotlib (pip install matplotlib)

Instalação

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

  1. pair_hamiltonian_circuit: Implementa o UIU_I o unitário de dois qubits para o hamiltoniano de interação aproximado entre sites vizinhos. A decomposição em 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 o operador unitário de dois qubits UEU_E para a energia aproximada do campo elétrico em cada local. A decomposição em 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, combinando termos de interação, elétricos e de massa com portas SWAP para gerenciar a conectividade dos 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

Exemplo de simulador em pequena escala

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

Etapa 1: Mapeamento de 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 do 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 etapa do algoritmo de Trotter, crie dois circuitos : um para inicializar um méson no centro (inverse_mid=True) e outro para preparar 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ádrão.

# 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

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

Defina as grandezas observáveis: medições de “ ZZ ” de um único qubit em cada qubit. A partir de ⟨Z⟩\langle Z \rangle, é possível extrair as probabilidades de ocupação e, em seguida, o número de férmions escalonado nf(r)n_f(r) em cada vértice da 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)}")

Output:

Number of observables: 12

Etapa 3: Executar usando Qiskit primitives

Use StatevectorEstimator para uma simulação exata e 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")

Output:

Computed expectation values for 10 Trotter steps

Etapa 4: Realizar o pós-processamento e apresentar o resultado no formato clássico desejado

Converta os valores de expectativa para o 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 gerar o mapa de calor da propagação do hádrão. Isso reproduz a estrutura da Figura 3 do artigo de referência: o local na rede rr no eixo x, o passo de Trotter (tempo) tt no eixo y e nf(r,t)n_f(r,t) como 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:

Output of the previous code cell

Exemplo de hardware em grande escala

Agora ampliamos para uma rede de 30 nós (60 qubits) no hardwar IBM Quantum. Nessa escala, o circuito de 10 passos de Trotter compreende mais de 3.400 portas de dois qubits e 14.000 portas de um qubit.

Etapas 1 a 4 (reunidas em um único bloco de código)

Principais aspectos do fluxo de trabalho de hardware:

  • 10 passos do algoritmo de Trotter para os circuitos do méson e do vácuo (intercalados para minimizar o desvio)
  • Transpilacão com optimization_level=1 — o layout do circuito já é isomórfico à topologia do dispositivo (uma cadeia linear), portanto, não são necessárias alterações de roteamento. O transpiler é utilizado exclusivamente para selecionar uma cadeia de qubits físicos com baixo ruído e decompor os portões no conjunto de portões nativos.
  • EstimatorV2 com mitigação de erros de leitura do TREX e giro de Pauli
  • Batch sessão para enviar todos os trabalhos de uma só 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

Avaliação comparativa clássica por meio da propagação de Pauli

O Método de Propagação de Pauli (PPM) fornece uma simulação clássica sem ruído do circuito quântico por meio da retropropagação de observáveis medidos ao longo do circuito na perspectiva de Heisenberg. Nas camadas de Clifford (portas CNOT, H, S e X), os operadores de Pauli são mapeados para outros operadores de Pauli sem aumentar o número de termos. Camadas não-Clifford (as portas do tipo “ RzR_z ” no circuito) podem causar ramificações — no pior dos casos, dobrando o número de termos —, mas muitas ramificações têm coeficientes pequenos e podem ser truncadas.

O fluxo de trabalho com pauli-prop é o seguinte:

  1. Divida o circuito em suas partes Clifford e não-Clifford usando evolve_through_cliffords.
  2. atolPropague cada observável pela parte não-Clifford usando propagate_through_circuit, mantendo até max_terms termos de Pauli e descartando termos com coeficientes abaixo do limiar de truncamento.
  3. Aplique a operação Clifford ao resultado utilizando o suporte integrado ao Clifford do Qiskit.
  4. Calcule o valor esperado somando os coeficientes dos termos diagonais de Pauli (que contêm apenas II e ZZ ).

Limite de truncamento

O atol parâmetro em propagate_through_circuit determina o grau de agressividade com que os pequenos ramos de Pauli são podados. Um limiar muito restrito (por exemplo, 1e-12) mantém quase todos os ramos e fornece resultados exatos, mas o tempo de simulação aumenta acentuadamente com a profundidade do circuito; a simulação de 120 qubits apresentada 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 monitorados e acelerando o cálculo. A contrapartida é um erro de aproximação pequeno e controlável, que você pode validar comparando os 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}")

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óximas etapas

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

Recomendações

Referências

[1] O artigo original: Ilčić, Majumdar, Mathew et al. “Observação de dinâmica hadrônica não abeliana robusta e coerente em processadores quânticos sujeitos a ruído” arXiv:2602.18080 (2026)

Esta página foi útil?
Relate um bug, erro de digitação ou solicite conteúdo no GitHub.