Simule o espalhamento de nêutrons em materiais quânticos com circuitos quânticos

Estimativa de uso: 13 minutos em um processador Heron r2 (NOTA: Esta é apenas uma estimativa. O seu tempo de execução pode variar.)

Qual tutorial devo usar?

Use este tutorial para aprender a implementação passo a passo, com a computação clássica sendo executada localmente em seu laptop e a execução em hardware em uma QPU IBM Quantum®. O exemplo em larga escala requer memória substancial, e a compilação quântica aproximada (AQC) pode levar várias horas em recursos típicos de laptop. Para transferir a compressão para os recursos de computação e memória do Qiskit Serverless, use o tutorial complementar do Serverless.

Resultados de aprendizagem

Como os espectros de espalhamento inelástico de nêutrons (INS) se conectam aos fatores de estrutura dinâmica (DSFs) de modelos de spin quânticos.
Como preparar um estado fundamental, aplicar uma perturbação local e realizar evolução temporal de Trotter em um circuito quântico.
Como usar compilação quântica aproximada (AQC) com qiskit-addon-aqc-tensor para comprimir circuitos de Trotter profundos para execução em hardware.
Como extrair a função de Green retardada (RGF) de valores esperados de qubits e transformá-la por Fourier em um DSF.

Pré-requisitos

Contexto

Neste tutorial, reproduzimos os resultados de Lee et al., arXiv:2603.15608.

O espalhamento de nêutrons ajuda os pesquisadores a caracterizar materiais relevantes para a sustentabilidade, incluindo eletrodos de baterias. Este tutorial explora como circuitos quânticos podem simular excitações magnéticas e conectar modelos teóricos a medições de espalhamento. Seu exemplo de ímã quântico desenvolve métodos para estudar o comportamento de materiais, contribuindo para o conjunto mais amplo de ferramentas de pesquisa voltado à descoberta de materiais.

Espalhamento inelástico de nêutrons e o fator de estrutura dinâmica

O espalhamento inelástico de nêutrons (INS) é uma das sondas experimentais mais poderosas de excitações magnéticas em materiais quânticos. Quando um feixe de nêutrons térmicos ou frios incide sobre um cristal, nêutrons individuais trocam momento $\mathbf{q}$ e energia $\omega$ com o subsistema magnético. A intensidade de espalhamento medida é proporcional ao fator de estrutura dinâmica (DSF),

$S^{\alpha\beta}(q,\omega) = \sum_{j} e^{-iq\,j}\int_{-\infty}^{\infty} dt\; e^{i\omega t}\, \langle S_0^{\alpha}(0)\, S_j^{\beta}(t)\rangle,$

que codifica as correlações espaço-temporais completas dos graus de liberdade de spin.

KCuF $_3$ : um ímã canônico de líquido de Luttinger

O fluoreto de cobre e potássio (KCuF $_3$ ) é um antiferromagneto quase unidimensional em que cadeias de íons Cu $^{2+}$ de spin- $\frac{1}{2}$ interagem por meio de uma troca de Heisenberg entre vizinhos mais próximos $J$ , enquanto o acoplamento entre cadeias é de apenas $\sim 2,7\%$ de $J$ . Em $T = 6\;\mathrm{K}$ , onde dados de INS estão disponíveis, o espectro é dominado por excitações fracionalizadas de spinons características de um líquido de Tomonaga-Luttinger. Como a dinâmica intracadeia é bem descrita pelo Hamiltoniano XXZ unidimensional de spin- $\frac{1}{2}$ no ponto isotrópico ( $\epsilon = 1$ ),

$H = J\sum_{i}\left[S_i^Z S_{i+1}^Z + \epsilon\left(S_i^X S_{i+1}^X + S_i^Y S_{i+1}^Y\right)\right],$

KCuF $_3$ serve como um benchmark ideal para simulação quântica: o Hamiltoniano é simples o suficiente para ser implementado em um processador quântico, mas o estado fundamental é fortemente emaranhado e o espectro de excitação apresenta um amplo contínuo de dois spinons.

Nota: este tutorial define $J = 1$ como unidade de energia e adota a normalização $H = J\sum_i[\ldots]$ , que corresponde ao Hamiltoniano local $H_\mathrm{loc} = J(S^X S^X + S^Y S^Y + S^Z S^Z)$ usado na implementação de circuito do artigo (Fig. S3 do suplementar). O Hamiltoniano completo do artigo (Eq. 3) carrega um fator adicional geral de 2, então o $J$ do artigo é o dobro do $J$ usado aqui.

O que simulamos e medimos

A quantidade física que calculamos é a função de Green retardada (RGF), definida como a função de correlação spin-spin dependente do tempo

$G^R_{\alpha,\beta}(j, j_c, t) = -\frac{i}{2}\,\langle\psi_{\mathrm{GS}}|\,S_j^\alpha(t)\,S_{j_c}^\beta(0) - S_{j_c}^\beta(0)\,S_j^\alpha(t)\,|\psi_{\mathrm{GS}}\rangle,$

onde $j_c$ é um sítio de referência (o centro da cadeia) e $S_j^\alpha(t) = e^{iHt}S_j^\alpha e^{-iHt}$ é o operador de spin na representação de Heisenberg. Neste tutorial, focamos na componente $zz$ ( $\alpha = \beta = z$ ). Em um computador quântico, a RGF é acessada preparando o estado fundamental, aplicando uma perturbação local em $j_c$ , evoluindo no tempo o estado perturbado e medindo o valor esperado de qubit único $\langle\sigma_j^z\rangle$ em cada sítio $j$ para cada passo de tempo. O insight principal é que cada $\langle\sigma_j^z\rangle$ fornece a diferença em relação à magnetização do estado fundamental; como o antiferromagneto de Heisenberg isotrópico tem magnetização líquida nula por sítio ( $\langle\sigma_j^z\rangle_\mathrm{GS} = 0$ ), o valor bruto medido produz diretamente $G^R(j, j_c, t)$ sem qualquer subtração explícita.

Ao coletar $G^R(j, j_c, t)$ em todos os sítios e passos de tempo, montamos um conjunto de dados bidimensional que é então transformado por Fourier tanto no espaço quanto no tempo para produzir o fator de estrutura dinâmica $S(q,\omega)$ . O DSF é a quantidade diretamente medida em um experimento de INS: ela nos diz quais excitações magnéticas existem em cada momento $q$ e energia $\omega$ . Para a cadeia de Heisenberg isotrópica, o espectro de excitação exato é um contínuo de dois spinons, uma ampla banda de intensidade de espalhamento cuja forma serve como um rigoroso benchmark de ponta a ponta para a simulação quântica: ela valida a preparação do estado fundamental, a perturbação, a evolução temporal de Trotter e o protocolo de medição, tudo de uma vez.

Fluxo de trabalho de simulação quântica

O fluxo de trabalho espelha a física de um evento de INS. Nós (1) preparamos o estado fundamental de muitos corpos $|\psi_{\mathrm{GS}}\rangle$ em $n$ qubits, (2) aplicamos uma perturbação local de inversão de spin $U_{j_c} = \frac{1}{\sqrt{2}}(I - i\sigma^z_{j_c})$ no centro da cadeia para imitar a transferência de spin do nêutron, (3) evoluímos sob $H$ para passos de tempo discretos usando trotterização de segunda ordem, e (4) medimos $\langle\sigma_i^z\rangle$ em cada qubit em cada passo para obter a RGF. Uma transformada de Fourier discreta bidimensional então produz o DSF $S(q,\omega)$ .

O observável que medimos em cada passo de tempo é $\sigma_i^z$ em cada qubit $i$ . No Qiskit, isso é representado como uma lista de operadores SparsePauliOp: um único $Z$ de qubit único embutido na string identidade de $n$ qubits para cada sítio. Esses observáveis são construídos uma vez por instância do problema durante a etapa de mapeamento do problema (Etapa 1) e reutilizados para cada circuito nessa escala.

Compilação quântica aproximada (AQC)

Circuitos de Trotter profundos podem ser comprimidos com compilação quântica aproximada (AQC), que substitui as primeiras camadas de Trotter por um ansatz parametrizado mais curto cujos parâmetros são otimizados classicamente para maximizar a fidelidade em nível de MPS com o circuito profundo original. Os passos de Trotter restantes são anexados exatamente, produzindo um circuito "misto" AQC + Trotter com substancialmente menos portas de dois qubits.

Simulação MPS

Para um sistema unidimensional, métodos de estado de produto de matrizes (MPS) podem simular eficientemente tanto a preparação do estado fundamental (com o grupo de renormalização de matriz de densidade, ou DMRG) quanto a evolução temporal em nível de circuito. Ao controlar a dimensão de vínculo $\chi$ , equilibramos precisão e custo computacional. Neste tutorial, usamos simulação MPS dentro de qiskit-addon-aqc-tensor para calcular ansätze AQC de alta fidelidade que comprimem os circuitos de Trotter profundos para execução em hardware.

Requisitos

Antes de iniciar este tutorial, certifique-se de ter instalado o seguinte:

Qiskit SDK com suporte a visualização
Qiskit Runtime (pip install qiskit-ibm-runtime)
qiskit-addon-aqc-tensor com os extras quimb e JAX (pip install 'qiskit-addon-aqc-tensor[quimb-jax]')

Setup

# Added by doQumentation — installs the packages this notebook needs if they are missing
import importlib.util

_needed = {"matplotlib": "matplotlib", "numpy": "numpy", "qiskit": "qiskit", "qiskit_addon_aqc_tensor": "qiskit-addon-aqc-tensor", "qiskit_ibm_runtime": "qiskit-ibm-runtime", "qiskit_quimb": "qiskit-quimb", "quimb": "quimb", "scipy": "scipy"}
_missing = [pip for module, pip in _needed.items()
            if importlib.util.find_spec(module) is None]
# One at a time, so a package that fails to install does not block the others
for _pip in _missing:
    %pip install -q {_pip}
if not _missing:
    print("\u2713 All required packages are installed")

import timeit
import warnings
from collections.abc import Iterator, Sequence
from functools import partial

import matplotlib.pyplot as plt
import numpy as np
import quimb.tensor as qtn
import scipy.optimize
from numpy.typing import NDArray
from qiskit import QuantumCircuit
from qiskit.circuit import CircuitInstruction, ParameterVector, Qubit
from qiskit.circuit.library import PauliEvolutionGate
from qiskit.primitives import StatevectorEstimator
from qiskit.quantum_info import SparsePauliOp
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
from qiskit_addon_aqc_tensor import generate_ansatz_from_circuit
from qiskit_addon_aqc_tensor.objective import MaximizeStateFidelity
from qiskit_addon_aqc_tensor.simulation import tensornetwork_from_circuit
from qiskit_addon_aqc_tensor.simulation.quimb import QuimbSimulator
from qiskit_ibm_runtime import EstimatorV2 as Estimator
from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_quimb import quimb_circuit
from scipy.sparse import SparseEfficiencyWarning

# scipy.sparse.linalg.expm (reached via PauliEvolutionGate.to_matrix) prefers
# CSC and warns when handed the CSR matrix from SparsePauliOp.to_matrix; the
# result is unaffected, so silence the cosmetic warning.
warnings.filterwarnings("ignore", category=SparseEfficiencyWarning)

def xxz_hamiltonian_mpo(
    n_qubits: int, interaction: float = 1.0, anisotropy: float = 1.0
) -> qtn.MatrixProductOperator:
    """1D XXZ Hamiltonian as a quimb MPO.

    Builds the Hamiltonian using ``qtn.SpinHam1D``.

    Args:
        n_qubits: Number of sites.
        interaction: Overall interaction strength (J in the paper).
        anisotropy: XY/Z anisotropy (ε in the paper). ``ε = 1`` is the
            isotropic Heisenberg point; ``ε = 0`` is the Ising limit.

    Returns:
        The Hamiltonian as a matrix product operator.
    """
    builder = qtn.SpinHam1D(S=1 / 2)
    builder += interaction * anisotropy * 0.5, "+", "-"
    builder += interaction * anisotropy * 0.5, "-", "+"
    builder += interaction, "Z", "Z"
    return builder.build_mpo(L=n_qubits)

def build_ground_state_ansatz(n_qubits: int, n_layers: int) -> QuantumCircuit:
    """Hamiltonian variational ansatz (HVA) circuit for ground-state preparation.

    Starts from a product of singlet pairs and applies alternating
    odd/even layers of parameterized XXZ pair-evolution gates.

    The returned circuit is parameterized: it carries a ``ParameterVector``
    named ``"theta"`` of length ``2 * n_layers`` whose values must be
    assigned (e.g. via ``circuit.assign_parameters``) before simulation.
    Element ``2 * r`` is the odd-layer angle and element ``2 * r + 1`` is the
    even-layer angle of layer ``r``.

    Args:
        n_qubits: Number of qubits (must be even).
        n_layers: Number of HVA layers.

    Returns:
        The parameterized HVA preparation circuit.
    """
    theta = ParameterVector("theta", 2 * n_layers)
    circuit = QuantumCircuit(n_qubits)
    # Initial singlet product state
    for i in range(n_qubits // 2):
        circuit.x(2 * i)
        circuit.x(2 * i + 1)
        circuit.h(2 * i + 1)
        circuit.cx(2 * i + 1, 2 * i)
    # Variational layers: each pair gate is exp(-i theta (XX+YY+ZZ)/2)
    pair_ham = SparsePauliOp(
        ["XX", "YY", "ZZ"], coeffs=[0.5, 0.5, 0.5]
    )  # H_pair (HVA form)
    for r in range(n_layers):
        for i in range(1, (n_qubits + 1) // 2):  # odd layer
            circuit.append(
                PauliEvolutionGate(pair_ham, time=theta[2 * r]),
                [2 * i - 1, 2 * i],
            )
        for i in range(n_qubits // 2):  # even layer
            circuit.append(
                PauliEvolutionGate(pair_ham, time=theta[2 * r + 1]),
                [2 * i, 2 * i + 1],
            )
    return circuit

def optimize_ground_state_ansatz(
    ansatz: QuantumCircuit,
    x0: NDArray[np.floating],
    target_mps: qtn.MatrixProductState,
    *,
    max_bond: int | None = None,
    cutoff: float = 1e-10,
    method: str = "COBYQA",
    options: dict | None = None,
) -> scipy.optimize.OptimizeResult:
    """Optimize HVA parameters by maximizing fidelity with a target MPS.

    The HVA circuit is simulated as a matrix product state with the given
    bond-dimension truncation, and the parameters are optimized to maximize
    the state fidelity ``|<psi_HVA | target_mps>|**2`` with the DMRG ground
    state. Both states are normalized, so the minimized objective is the
    infidelity ``1 - |<psi_HVA | target_mps>|**2``.

    Args:
        ansatz: Parameterized HVA circuit from ``build_ground_state_ansatz``.
            The length of ``x0`` must equal ``ansatz.num_parameters``.
        x0: Initial parameters.
        target_mps: Target MPS (DMRG ground state) to maximize fidelity with.
        max_bond: Maximum MPS bond dimension during gate application.
        cutoff: Singular-value cutoff during gate application.
        method: ``scipy.optimize.minimize`` method.
        options: Options dict forwarded to ``scipy.optimize.minimize``.

    Returns:
        The Scipy OptimizeResult.
    """

    def infidelity(params: NDArray[np.floating]) -> float:
        circuit = ansatz.assign_parameters(params)
        circuit_mps = quimb_circuit(
            circuit.decompose(["PauliEvolution"]),
            quimb_circuit_class=qtn.CircuitMPS,
            max_bond=max_bond,
            cutoff=cutoff,
        )
        return 1 - abs(circuit_mps.psi.H @ target_mps) ** 2

    return scipy.optimize.minimize(
        infidelity, np.asarray(x0), method=method, options=options
    )

def trotter_evolution(
    qubits: Sequence[Qubit],
    interaction: float,
    anisotropy: float,
    time_step: float,
    n_steps: int,
) -> Iterator[CircuitInstruction]:
    """Second-order Trotter steps of the XXZ pair Hamiltonian.

    While the paper used a hand-optimized circuit for the Trotter steps, we use
    PauliEvolutionGate here for simplicity and generality. The final two-qubit gate
    count and gate depth are equivalent when transpiled with ``optimization_level=3``.

    Args:
        qubits: Qubits to act on (length ``n_qubits``).
        interaction: Overall interaction strength (J in the paper).
        anisotropy: XY/Z anisotropy (ε in the paper). ``ε = 1`` is the
            isotropic Heisenberg point; ``ε = 0`` is the Ising limit.
        time_step: Per-step Trotter time.
        n_steps: Number of Trotter steps.

    Yields:
        ``CircuitInstruction``s implementing the Trotter steps.
    """
    if n_steps == 0:
        return
    n_qubits = len(qubits)
    pair_ham = SparsePauliOp(
        ["XX", "YY", "ZZ"],
        coeffs=[
            0.25 * interaction * anisotropy,
            0.25 * interaction * anisotropy,
            0.25 * interaction,
        ],
    )
    half_evo = PauliEvolutionGate(pair_ham, time=time_step / 2)
    full_evo = PauliEvolutionGate(pair_ham, time=time_step)
    for i in range(n_qubits // 2):  # half even layer
        yield CircuitInstruction(half_evo, (qubits[2 * i], qubits[2 * i + 1]))
    for i in range(n_qubits // 2 - 1):  # full odd layer
        yield CircuitInstruction(
            full_evo, (qubits[2 * i + 1], qubits[2 * i + 2])
        )
    for _ in range(n_steps - 1):  # interior steps
        for i in range(n_qubits // 2):
            yield CircuitInstruction(
                full_evo, (qubits[2 * i], qubits[2 * i + 1])
            )
        for i in range(n_qubits // 2 - 1):
            yield CircuitInstruction(
                full_evo, (qubits[2 * i + 1], qubits[2 * i + 2])
            )
    for i in range(n_qubits // 2):  # half even layer
        yield CircuitInstruction(half_evo, (qubits[2 * i], qubits[2 * i + 1]))

def get_dsf(
    n_qubits: int,
    rgf_mat: NDArray[np.floating],
    time_step: float,
    n_steps: int,
    n_points_momentum: int,
    n_points_frequency: int,
) -> NDArray[np.floating]:
    """Compute the dynamical structure factor from the retarded Green's function.

    Uses the center-site approximation and a discrete Fourier transform.
    The result is symmetrized about the momentum axis and clipped to
    non-negative values, ready for plotting.

    Args:
        n_qubits: Number of qubits (sites).
        rgf_mat: RGF matrix of shape ``(n_steps, n_qubits)``.
        time_step: Trotter time-step size.
        n_steps: Number of time steps.
        n_points_momentum: Number of momentum points.
        n_points_frequency: Number of frequency points.

    Returns:
        DSF array of shape ``(n_points_frequency, n_points_momentum)``,
        symmetrized about the momentum axis and clipped to non-negative values.
    """
    max_frequency = np.pi / time_step
    momentum_range = np.linspace(0, 2 * np.pi, n_points_momentum)
    frequency_range = np.linspace(0, max_frequency, n_points_frequency)
    result = np.zeros((frequency_range.shape[0], momentum_range.shape[0]))
    center = n_qubits // 2 - 1
    for iw, w in enumerate(frequency_range):
        exponent = np.exp(1j * w * time_step * np.arange(1, n_steps + 1))
        # S = sigma/2, so two spin operators contribute a factor of (1/2)(1/2) = 1/4.
        rgf_omega = (
            np.dot(rgf_mat.T, exponent) * time_step / 4
        )  # S(omega): time Fourier slice of the Green's function
        for iq, q in enumerate(momentum_range):
            momentum_phases = np.exp(
                -1j * q * np.arange(-center, center + 2, 1)
            )
            result[iw, iq] = np.imag(np.dot(rgf_omega, momentum_phases))
    result = -(result + result[:, ::-1]) / 2
    result = np.clip(result, a_min=0, a_max=None)
    return result

def plot_dsf(
    dsf: NDArray[np.floating],
    time_step: float,
    n_points_momentum: int,
    n_points_frequency: int,
    title: str | None = None,
) -> None:
    """Heat-map of the dynamical structure factor.

    Args:
        dsf: DSF array of shape ``(n_points_frequency, n_points_momentum)``.
        time_step: Trotter time-step size.
        n_points_momentum: Number of momentum points.
        n_points_frequency: Number of frequency points.
        title: Optional plot title.
    """
    max_frequency = np.pi / time_step
    momentum_range = np.linspace(0, 2 * np.pi, n_points_momentum)
    frequency_range = np.linspace(0, max_frequency, n_points_frequency)
    x, y = np.meshgrid(momentum_range, frequency_range)
    fig, ax = plt.subplots(figsize=(8, 5))
    c = ax.pcolormesh(x, y, dsf / np.max(dsf), cmap="viridis", shading="auto")
    fig.colorbar(c, ax=ax, label="Normalized intensity")
    ax.set_ylim(0, 3.6)
    ax.set_xlim(0, 2 * np.pi)
    ax.set_xlabel(r"$q$", fontsize=16)
    ax.set_ylabel(r"$\tilde{\omega} = \omega / J$", fontsize=16)
    ax.set_xticks([0, np.pi / 2, np.pi, 3 * np.pi / 2, 2 * np.pi])
    ax.set_xticklabels(["0", r"$\pi/2$", r"$\pi$", r"$3\pi/2$", r"$2\pi$"])
    if title:
        ax.set_title(title, fontsize=14)
    plt.tight_layout()
    plt.show()

def plot_rgf(
    n_qubits: int,
    rgf_mat: NDArray[np.floating],
    time_step: float,
    n_steps: int,
    title: str | None = None,
) -> None:
    """Heat-map of the retarded Green's function in real space and time.

    Args:
        n_qubits: Number of qubits (sites).
        rgf_mat: RGF matrix of shape ``(n_steps, n_qubits)``.
        time_step: Trotter time-step size.
        n_steps: Number of time steps.
        title: Optional plot title.
    """
    fig, ax = plt.subplots(figsize=(8, 6))
    qubit_axis = np.arange(n_qubits)
    t_axis = np.arange(1, n_steps + 1) * time_step
    x, y = np.meshgrid(qubit_axis, t_axis)
    c = ax.pcolormesh(
        x,
        y,
        np.real(rgf_mat),
        cmap="RdBu",
        vmax=0.5,
        vmin=-0.5,
        shading="auto",
    )
    fig.colorbar(c, ax=ax, label=r"Re $G^R(j, j_c, t)$")
    ax.set_xlabel("Qubit", fontsize=16)
    ax.xaxis.set_major_locator(
        plt.matplotlib.ticker.MaxNLocator(integer=True)
    )
    ax.set_ylabel(r"Time", fontsize=16)
    if title:
        ax.set_title(title, fontsize=14)
    plt.tight_layout()
    plt.show()

def uniform_2q_depth(circuit: QuantumCircuit) -> int:
    """Two-qubit gate depth in a standardized basis."""
    pass_manager = generate_preset_pass_manager(
        optimization_level=0, basis_gates=["cz", "id", "rz", "sx", "x"]
    )
    return pass_manager.run(circuit).depth(
        lambda inst: inst.operation.num_qubits == 2
    )

Exemplo de simulador em pequena escala

Primeiro demonstramos o fluxo de trabalho completo em 10 qubits, otimizando o ansatz de estado fundamental HVA usando simulação MPS e usando o simulador statevector do Qiskit para a evolução temporal. O DMRG fornece uma energia de estado fundamental de referência e uma MPS de referência. Este exemplo em pequena escala nos permite validar cada etapa antes de escalar.

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

Começamos definindo o modelo físico e construindo os circuitos quânticos.

Hamiltoniano. KCuF $_3$ é modelado pelo Hamiltoniano XXZ 1D no ponto isotrópico ( $J = 1$ , $\epsilon = 1$ , definindo $J = 1$ como unidade de energia).

Estado fundamental. Usamos um circuito de ansatz variacional Hamiltoniano (HVA), build_ground_state_ansatz, como circuito de preparação do estado fundamental. Os parâmetros do HVA são otimizados classicamente maximizando a fidelidade do estado $|\langle\psi_{\mathrm{HVA}}(\theta)|\psi_{\mathrm{DMRG}}\rangle|^2$ com a MPS do estado fundamental do DMRG, onde o estado HVA é avaliado simulando o circuito como um estado de produto de matrizes com o CircuitMPS do quimb. Usamos scipy.optimize.minimize para a otimização. A energia $\langle H \rangle$ do ansatz otimizado também é calculada como referência.

Portas de Trotter. Cada termo de interação entre vizinhos mais próximos $e^{-i\Delta t\, H_{\mathrm{pair}}}$ com $H_{\mathrm{pair}} = (J/4)(XX + YY + ZZ)$ é construído com PauliEvolutionGate(H_pair, time=time_step). O Qiskit sintetiza isso na decomposição ótima de três CNOTs durante a transpilação.

Perturbação. Uma porta $R_z(\pi/2)$ no qubit central implementa $U_{j_c} = \frac{1}{\sqrt{2}}(I - i\sigma^z_{j_c})$ , imitando a inversão de spin local produzida por um nêutron espalhado.

Observáveis. Construímos um observável $\sigma^z$ para cada sítio de qubit. Esses objetos SparsePauliOp são passados para a primitiva Estimator na Etapa 3 para extrair $\langle\sigma_i^z\rangle$ em cada passo de tempo.

# -- Physical parameters --
n_qubits = 10
interaction = 1.0  # J
anisotropy = 1.0  # ε (isotropic point)
time_step = 0.6
n_steps = 10
mps_max_bond = 32
mps_cutoff = 1e-8
center = n_qubits // 2 - 1

# -- Hamiltonian MPO --
ham_mpo = xxz_hamiltonian_mpo(n_qubits, interaction, anisotropy)

# -- Reference ground-state energy via DMRG --
dmrg = qtn.DMRG2(ham_mpo)
dmrg.solve(tol=1e-8)
print(f"Ground-state energy (DMRG): {dmrg.energy:.6f}")

# -- Build ground state ansatz circuit --
gs_n_layers = 3
gs_ansatz = build_ground_state_ansatz(n_qubits, gs_n_layers)

# -- Optimize ground state ansatz parameters to maximize fidelity with the DMRG MPS --
# Initialize odd-layer angles near 0 (where the inter-pair gate is the
# identity) and even-layer angles near pi/2 (where the intra-pair gate
# is a SWAP, since 0.5*(XX + YY + ZZ) = SWAP - I/2).
rng = np.random.default_rng(12345)
x0 = np.tile([0, np.pi / 2], gs_n_layers) + rng.normal(
    scale=0.1, size=2 * gs_n_layers
)
print("Optimizing ground state ansatz...")
t0 = timeit.default_timer()
result = optimize_ground_state_ansatz(
    gs_ansatz,
    x0,
    dmrg.state,
    max_bond=mps_max_bond,
    cutoff=mps_cutoff,
    options=dict(maxiter=100),
)
t1 = timeit.default_timer()
print(f"Finished optimizing ground state ansatz in {t1 - t0} seconds.")
print(f"Ground state ansatz fidelity: {1 - result.fun:.6f}")

gs_circuit = gs_ansatz.assign_parameters(result.x)
gs_circuit_mps = quimb_circuit(
    gs_circuit.decompose(["PauliEvolution"]),
    quimb_circuit_class=qtn.CircuitMPS,
    max_bond=mps_max_bond,
    cutoff=mps_cutoff,
)
gs_ansatz_energy = qtn.expec_TN_1D(
    gs_circuit_mps.psi.H, ham_mpo, gs_circuit_mps.psi
)
print(f"Ground state ansatz energy: {gs_ansatz_energy:.6f}")

# -- Build circuits for each time step --
perturbed = gs_circuit.copy()
perturbed.rz(np.pi / 2, center)

circuits = []
for t in range(1, n_steps + 1):
    circuit = perturbed.copy()
    for instr in trotter_evolution(
        circuit.qubits, interaction, anisotropy, time_step, t
    ):
        circuit.append(instr)
    circuits.append(circuit)
print(
    f"Built {len(circuits)} circuits, deepest 2q depth (uniform basis) = "
    f"{uniform_2q_depth(circuits[-1])}"
)

# -- Observables: Z on each qubit site --
observables = [
    SparsePauliOp.from_sparse_list([("Z", [i], 1)], num_qubits=n_qubits)
    for i in range(n_qubits)
]
print(f"Defined {len(observables)} Z observables.")

Ground-state energy (DMRG): -4.258035
Optimizing ground state ansatz...
Finished optimizing ground state ansatz in 8.850687621976249 seconds.
Ground state ansatz fidelity: 0.984277
Ground state ansatz energy: -4.232565
Built 10 circuits, deepest 2q depth (uniform basis) = 163
Defined 10 Z observables.

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

Para hardware real, os circuitos de Trotter acima seriam muito profundos. A compilação quântica aproximada (AQC) resolve isso substituindo as primeiras $k$ camadas de Trotter (incluindo o circuito de estado fundamental) por um ansatz parametrizado mais curto, otimizado para maximizar a fidelidade em nível de MPS com o circuito profundo original. Os passos de Trotter restantes são anexados exatamente, produzindo um circuito "AQC + Trotter" mais raso.

Para equilibrar expressividade e profundidade do circuito, dois ansätze são usados: um ansatz de uma camada (gerado a partir de um único passo de Trotter) comprime os primeiros passos de tempo, e um ansatz de duas camadas mais profundo (gerado a partir de dois passos de Trotter) comprime os próximos passos, onde é necessária maior fidelidade.

O fluxo de trabalho AQC tem quatro subetapas:

Construir circuitos alvo — os primeiros $k_1 + k_2$ circuitos da Etapa 1 servem diretamente como alvos AQC.
Calcular a MPS alvo — simule cada circuito alvo como um estado de produto de matrizes usando quimb.tensor.CircuitMPS.
Gerar e otimizar ansätze — generate_ansatz_from_circuit cria um ansatz de uma camada e um de duas camadas; os parâmetros são otimizados usando L-BFGS-B com gradientes acelerados por JAX para minimizar $1 - |\langle\psi_{\mathrm{ansatz}}|\psi_{\mathrm{target}}\rangle|^2$ . Os primeiros $k_1$ passos usam o ansatz de uma camada e os próximos $k_2$ passos usam o ansatz de duas camadas; dentro de cada estágio, cada passo é inicializado com os parâmetros otimizados do passo anterior, e os parâmetros são redefinidos para os padrões do estágio no limite entre estágios.
Montar circuitos mistos — para passos de tempo além do ponto de verificação do AQC, anexe camadas de Trotter exatas ao circuito AQC de duas camadas otimizado usando trotter_evolution.

# Number of time steps to compress into an AQC ansatz with one layer
aqc_n_steps_1 = 3
# Number of time steps to compress into an AQC ansatz with two layers
aqc_n_steps_2 = 2
aqc_n_steps_total = aqc_n_steps_1 + aqc_n_steps_2

# ── Step 2a: Target circuits (first aqc_n_steps_total circuits from Step 1) ──
target_circuits = {
    k: circuits[k - 1] for k in range(1, aqc_n_steps_total + 1)
}
# ── Step 2b: Compute target MPS ──
# Decompose PauliEvolutionGate to RXX/RYY/RZZ before passing to the AQC MPS
# backend, which does not understand PauliEvolutionGate natively.
aqc_sim = QuimbSimulator(
    quimb_circuit_factory=partial(
        qtn.CircuitMPS,
        gate_opts=dict(max_bond=mps_max_bond, cutoff=mps_cutoff),
    ),
    autodiff_backend="jax",
)
print("\nStep 2b — target MPS:")
target_mps = {}
for k in range(1, aqc_n_steps_total + 1):
    target_mps[k] = tensornetwork_from_circuit(
        target_circuits[k].decompose(["PauliEvolution"]), aqc_sim
    )
    print(f"  k={k}: max bond = {target_mps[k].psi.max_bond()}")

# ── Step 2c: Generate ansätze (one-layer and two-layer) and optimize parameters ──
ansatz_1, initial_params_1 = generate_ansatz_from_circuit(
    target_circuits[1].decompose(["PauliEvolution"]),
    qubits_initially_zero=True,
)
initial_params_1 = np.array(initial_params_1)
ansatz_2, initial_params_2 = generate_ansatz_from_circuit(
    target_circuits[2].decompose(["PauliEvolution"]),
    qubits_initially_zero=True,
)
initial_params_2 = np.array(initial_params_2)
print(
    f"\nStep 2c — one-layer ansatz: {ansatz_1.num_parameters} parameters, "
    f"2q depth (uniform basis) = {uniform_2q_depth(ansatz_1)}"
)
print(
    f"Step 2c — two-layer ansatz: {ansatz_2.num_parameters} parameters, "
    f"2q depth (uniform basis) = {uniform_2q_depth(ansatz_2)}"
)

aqc_circuits = {}
aqc_params = {}
for k in range(1, aqc_n_steps_total + 1):
    if k <= aqc_n_steps_1:
        ansatz, base_params = ansatz_1, initial_params_1
    else:
        ansatz, base_params = ansatz_2, initial_params_2
    # Warm-start from the previous step only within the same stage
    same_stage = (k - 1 >= 1) and (
        (k - 1 <= aqc_n_steps_1) == (k <= aqc_n_steps_1)
    )
    x0 = aqc_params[k - 1] if same_stage else base_params
    obj = MaximizeStateFidelity(target_mps[k], ansatz, aqc_sim)
    t0 = timeit.default_timer()
    result = scipy.optimize.minimize(
        obj.loss_function,
        x0,
        method="L-BFGS-B",
        jac=True,
        options=dict(maxiter=100),
    )
    elapsed = timeit.default_timer() - t0
    aqc_params[k] = result.x
    aqc_circuits[k] = ansatz.assign_parameters(result.x)
    print(
        f"  k={k}: fidelity = {1 - result.fun:.4f}, "
        f"2q depth (uniform basis) = "
        f"{uniform_2q_depth(aqc_circuits[k])}, "
        f"{elapsed:.1f}s"
    )

# ── Step 2d: Assemble full circuit set (AQC + Trotter) ──
all_circuits = []
for k in range(1, aqc_n_steps_total + 1):
    all_circuits.append(aqc_circuits[k])
base = aqc_circuits[aqc_n_steps_total]
for k in range(1, n_steps - aqc_n_steps_total + 1):
    circuit = base.copy()
    for instr in trotter_evolution(
        circuit.qubits, interaction, anisotropy, time_step, k
    ):
        circuit.append(instr)
    all_circuits.append(circuit)

full_depths = [
    uniform_2q_depth(circuits[k - 1]) for k in range(1, n_steps + 1)
]
aqc_depths = [uniform_2q_depth(circuit) for circuit in all_circuits]
full_2q = [
    circuits[k - 1].decompose(["PauliEvolution"]).num_nonlocal_gates()
    for k in range(1, n_steps + 1)
]
aqc_2q = [
    circuit.decompose(["PauliEvolution"]).num_nonlocal_gates()
    for circuit in all_circuits
]
print(f"\nStep 2d — assembled {len(all_circuits)} circuits")
print(
    f"  At step {n_steps} (uniform basis): "
    f"Trotter 2q depth = {full_depths[-1]}, AQC+Trotter 2q depth = {aqc_depths[-1]}; "
    f"2q gates: Trotter = {full_2q[-1]}, AQC+Trotter = {aqc_2q[-1]}"
)

steps_axis = np.arange(1, n_steps + 1)
fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(steps_axis, full_depths, "-o", color="black", label="Full Trotter")
ax.plot(
    steps_axis, aqc_depths, "-o", color="cadetblue", label="AQC + Trotter"
)
ax.set_xlabel("Trotter steps", fontsize=13)
ax.xaxis.set_major_locator(plt.matplotlib.ticker.MaxNLocator(integer=True))
ax.set_ylabel("2q gate depth (uniform basis)", fontsize=13)
ax.set_title(f"AQC circuit-depth reduction ({n_qubits} qubits)", fontsize=13)
ax.legend(fontsize=11)
plt.tight_layout()
plt.show()

Step 2b — target MPS:
  k=1: max bond = 22
  k=2: max bond = 22
  k=3: max bond = 26
  k=4: max bond = 27
  k=5: max bond = 30

Step 2c — one-layer ansatz: 389 parameters, 2q depth (uniform basis) = 27
Step 2c — two-layer ansatz: 470 parameters, 2q depth (uniform basis) = 33
  k=1: fidelity = 1.0000, 2q depth (uniform basis) = 27, 11.8s
  k=2: fidelity = 0.9980, 2q depth (uniform basis) = 27, 13.2s
  k=3: fidelity = 0.9890, 2q depth (uniform basis) = 27, 13.7s
  k=4: fidelity = 0.9983, 2q depth (uniform basis) = 33, 16.4s
  k=5: fidelity = 0.9957, 2q depth (uniform basis) = 33, 16.4s

Step 2d — assembled 10 circuits
  At step 10 (uniform basis): Trotter 2q depth = 163, AQC+Trotter 2q depth = 99; 2q gates: Trotter = 371, AQC+Trotter = 200

Output of the previous code cell

Etapa 3: Executar usando primitivas Qiskit

Simulamos cada circuito compilado com AQC usando a primitiva StatevectorEstimator.

estimator = StatevectorEstimator()
pubs = [(circuit, observables) for circuit in all_circuits]
job = estimator.run(pubs)
result = job.result()

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

Agora extraímos o valor esperado $\langle\sigma_i^z\rangle$ para cada qubit $i$ em cada passo de tempo $t$ . Esses valores formam a matriz da função de Green retardada $G^R(j, j_c, t)$ . Em seguida, transformamos por Fourier a RGF no fator de estrutura dinâmica $S(q, \omega)$ e plotamos tanto a RGF quanto o DSF. get_dsf aplica simetria de espelhamento e recorta valores negativos internamente.

rgf_mat = np.stack([pub_result.data.evs for pub_result in result])

# -- Compute DSF --
n_points_momentum, n_points_frequency = 100, 100
spectrum = get_dsf(
    n_qubits,
    rgf_mat,
    time_step,
    n_steps,
    n_points_momentum,
    n_points_frequency,
)

# -- Plot retarded Green's function --
plot_rgf(
    n_qubits,
    rgf_mat,
    time_step,
    n_steps,
    title=f"Retarded Green's function — {n_qubits} qubits (simulation)",
)

# -- Plot DSF --
plot_dsf(
    spectrum,
    time_step,
    n_points_momentum,
    n_points_frequency,
    title=f"Dynamical structure factor — {n_qubits} qubits (simulation)",
)

Output of the previous code cell

Execução em hardware em larga escala

Agora escalamos para 50 qubits. Nessa escala, as otimizações atingem fidelidades menores do que no exemplo em pequena escala: a fidelidade do ansatz de estado fundamental cai para aproximadamente 0,65, e as fidelidades AQC diminuem para aproximadamente 0,7 nos pontos de verificação posteriores. Isso é esperado, e você pode melhorar essas fidelidades aumentando o número de camadas do ansatz de estado fundamental (gs_n_layers) ou o número de iterações de otimização (maxiter), com custo clássico adicional. Observe também que as fidelidades AQC de duas camadas saem menores do que as de uma camada. Isso não é uma regressão: passos de tempo posteriores geram mais emaranhamento e são simplesmente mais difíceis de comprimir, por isso o ansatz de duas camadas, mais expressivo, é usado para eles.

A otimização AQC também pode levar várias horas de tempo de computação clássica (aproximadamente seis horas na execução mostrada aqui, a maior parte gasta nos pontos de verificação de duas camadas na Etapa 2c). Para reduzir o tempo total, considere executar este notebook em hardware clássico mais potente, como um sistema de computação de alto desempenho (HPC). Alternativamente, você pode reduzir a escala para uma instância de problema menor (por exemplo, menos qubits ou menos passos de tempo); seus resultados serão então diferentes dos mostrados aqui.

O código abaixo segue a mesma estrutura de quatro etapas do exemplo em pequena escala. Os parâmetros do estado fundamental HVA são novamente otimizados usando simulação MPS. Na QPU, ativamos o desacoplamento dinâmico (DD), o twirling de Pauli e a extinção de erros de leitura por twirling (TREX) para supressão e mitigação de erros. A tabela a seguir resume como o experimento em larga escala difere do de pequena escala:

	Pequena escala	Grande escala
Qubits	10	50
Passos de tempo	10	20
Pontos de verificação AQC (1 camada + 2 camadas)	3 + 2 = 5	6 + 4 = 10
Camadas do ansatz de estado fundamental	3	5
Dimensão máxima de vínculo MPS	32	128
Estimator	`StatevectorEstimator`	QPU com DD, twirling de Pauli e TREX

Mensagens de compilação lenta do XLA

Durante a otimização AQC (Etapa 2c abaixo), você pode ver mensagens de stderr como:

[Compiling module jit_MakeArrayFn for CPU] Very slow compile? ...
The operation took 2m14s

Estes são diagnósticos benignos do XLA, o compilador que sustenta a autodiferenciação JAX do qiskit-addon-aqc-tensor. Com 50 qubits e dimensão de vínculo MPS 128, o XLA gasta alguns minutos compilando a função de gradiente na primeira vez que é rastreada. A compilação é bem-sucedida, e os resultados da otimização não são afetados.

# ── Parameters ──────────────────────────────────────────────────────────────
n_qubits = 50  # 10 → 50
interaction = 1.0  # J
anisotropy = 1.0  # ε (isotropic point)
time_step = 0.6
n_steps = 20  # 10 → 20
aqc_n_steps_1 = 6  # 3 → 6
aqc_n_steps_2 = 4  # 2 → 4
aqc_n_steps_total = aqc_n_steps_1 + aqc_n_steps_2
gs_n_layers = 5  # 3 → 5
mps_max_bond = 128  # 32 -> 128
mps_cutoff = 1e-8
center = n_qubits // 2 - 1

# ── Step 1: Map ──────────────────────────────────────────────────────────────
ham_mpo = xxz_hamiltonian_mpo(n_qubits, interaction, anisotropy)

dmrg = qtn.DMRG2(ham_mpo)
dmrg.solve(tol=1e-8)
print(f"Ground-state energy (DMRG): {dmrg.energy:.6f}")

gs_ansatz = build_ground_state_ansatz(n_qubits, gs_n_layers)

rng = np.random.default_rng(12345)
x0 = np.tile([0.0, np.pi / 2], gs_n_layers) + rng.normal(
    scale=0.1, size=2 * gs_n_layers
)
print("Optimizing ground state ansatz...")
t0 = timeit.default_timer()
result = optimize_ground_state_ansatz(
    gs_ansatz,
    x0,
    dmrg.state,
    max_bond=mps_max_bond,
    cutoff=mps_cutoff,
    options=dict(maxiter=100),
)
t1 = timeit.default_timer()
print(f"Finished optimizing ground state ansatz in {t1 - t0} seconds.")
print(f"Ground state ansatz fidelity: {1 - result.fun:.6f}")

gs_circuit = gs_ansatz.assign_parameters(result.x)
gs_circuit_mps = quimb_circuit(
    gs_circuit.decompose(["PauliEvolution"]),
    quimb_circuit_class=qtn.CircuitMPS,
    max_bond=mps_max_bond,
    cutoff=mps_cutoff,
)
gs_ansatz_energy = qtn.expec_TN_1D(
    gs_circuit_mps.psi.H, ham_mpo, gs_circuit_mps.psi
)
print(f"Ground state ansatz energy: {gs_ansatz_energy:.6f}")

perturbed = gs_circuit.copy()
perturbed.rz(np.pi / 2, center)

circuits = []
for t in range(1, n_steps + 1):
    circuit = perturbed.copy()
    for instr in trotter_evolution(
        circuit.qubits, interaction, anisotropy, time_step, t
    ):
        circuit.append(instr)
    circuits.append(circuit)
print(
    f"Built {len(circuits)} circuits, deepest 2q depth (uniform basis) = "
    f"{uniform_2q_depth(circuits[-1])}"
)

observables = [
    SparsePauliOp.from_sparse_list([("Z", [i], 1)], num_qubits=n_qubits)
    for i in range(n_qubits)
]
print(f"Defined {len(observables)} Z observables.")

# ── Step 2: AQC ──────────────────────────────────────────────────────────────
target_circuits = {
    k: circuits[k - 1] for k in range(1, aqc_n_steps_total + 1)
}

aqc_sim = QuimbSimulator(
    quimb_circuit_factory=partial(
        qtn.CircuitMPS,
        gate_opts=dict(max_bond=mps_max_bond, cutoff=mps_cutoff),
    ),
    autodiff_backend="jax",
)
print("\nStep 2b — target MPS:")
target_mps = {}
for k in range(1, aqc_n_steps_total + 1):
    target_mps[k] = tensornetwork_from_circuit(
        target_circuits[k].decompose(["PauliEvolution"]), aqc_sim
    )
    print(f"  k={k}: max bond = {target_mps[k].psi.max_bond()}")

ansatz_1, initial_params_1 = generate_ansatz_from_circuit(
    target_circuits[1].decompose(["PauliEvolution"]),
    qubits_initially_zero=True,
)
initial_params_1 = np.array(initial_params_1)
ansatz_2, initial_params_2 = generate_ansatz_from_circuit(
    target_circuits[2].decompose(["PauliEvolution"]),
    qubits_initially_zero=True,
)
initial_params_2 = np.array(initial_params_2)
print(
    f"\nStep 2c — one-layer ansatz: {ansatz_1.num_parameters} parameters, "
    f"2q depth (uniform basis) = {uniform_2q_depth(ansatz_1)}"
)
print(
    f"Step 2c — two-layer ansatz: {ansatz_2.num_parameters} parameters, "
    f"2q depth (uniform basis) = {uniform_2q_depth(ansatz_2)}"
)

aqc_circuits = {}
aqc_params = {}
for k in range(1, aqc_n_steps_total + 1):
    if k <= aqc_n_steps_1:
        ansatz, base_params = ansatz_1, initial_params_1
    else:
        ansatz, base_params = ansatz_2, initial_params_2
    # Warm-start from the previous step only within the same stage
    same_stage = (k - 1 >= 1) and (
        (k - 1 <= aqc_n_steps_1) == (k <= aqc_n_steps_1)
    )
    x0 = aqc_params[k - 1] if same_stage else base_params
    obj = MaximizeStateFidelity(target_mps[k], ansatz, aqc_sim)
    t0 = timeit.default_timer()
    result = scipy.optimize.minimize(
        obj.loss_function,
        x0,
        method="L-BFGS-B",
        jac=True,
        options=dict(maxiter=100),
    )
    elapsed = timeit.default_timer() - t0
    aqc_params[k] = result.x
    aqc_circuits[k] = ansatz.assign_parameters(result.x)
    print(
        f"  k={k}: fidelity = {1 - result.fun:.4f}, "
        f"2q depth (uniform basis) = "
        f"{uniform_2q_depth(aqc_circuits[k])}, "
        f"{elapsed:.1f}s"
    )

all_circuits = []
for k in range(1, aqc_n_steps_total + 1):
    all_circuits.append(aqc_circuits[k])
base = aqc_circuits[aqc_n_steps_total]
for k in range(1, n_steps - aqc_n_steps_total + 1):
    circuit = base.copy()
    for instr in trotter_evolution(
        circuit.qubits, interaction, anisotropy, time_step, k
    ):
        circuit.append(instr)
    all_circuits.append(circuit)

full_depths = [
    uniform_2q_depth(circuits[k - 1]) for k in range(1, n_steps + 1)
]
aqc_depths = [uniform_2q_depth(circuit) for circuit in all_circuits]
full_2q = [
    circuits[k - 1].decompose(["PauliEvolution"]).num_nonlocal_gates()
    for k in range(1, n_steps + 1)
]
aqc_2q = [
    circuit.decompose(["PauliEvolution"]).num_nonlocal_gates()
    for circuit in all_circuits
]
print(f"\nStep 2d — assembled {len(all_circuits)} circuits")
print(
    f"  At step {n_steps} (uniform basis): "
    f"Trotter 2q depth = {full_depths[-1]}, AQC+Trotter 2q depth = {aqc_depths[-1]}; "
    f"2q gates: Trotter = {full_2q[-1]}, AQC+Trotter = {aqc_2q[-1]}"
)

steps_axis = np.arange(1, n_steps + 1)
fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(steps_axis, full_depths, "-o", color="black", label="Full Trotter")
ax.plot(
    steps_axis, aqc_depths, "-o", color="cadetblue", label="AQC + Trotter"
)
ax.set_xlabel("Trotter steps", fontsize=13)
ax.xaxis.set_major_locator(plt.matplotlib.ticker.MaxNLocator(integer=True))
ax.set_ylabel("2q gate depth (uniform basis)", fontsize=13)
ax.set_title(f"AQC circuit-depth reduction ({n_qubits} qubits)", fontsize=13)
ax.legend(fontsize=11)
plt.tight_layout()
plt.show()

# ── Step 3: Execute on IBM Quantum hardware ───────────────────────────────────
# (replaces StatevectorEstimator)
service = QiskitRuntimeService()
backend = service.least_busy(
    min_num_qubits=n_qubits,
    operational=True,
    simulator=False,
    filters=lambda x: x.configuration().processor_type["family"] == "Heron",
)
print(f"Backend: {backend.name}")

pm = generate_preset_pass_manager(optimization_level=3, backend=backend)
isa_circuits = pm.run(all_circuits, num_processes=1)
isa_2q_depths = [
    isa_circuit.depth(lambda inst: inst.operation.num_qubits == 2)
    for isa_circuit in isa_circuits
]
print(
    f"Transpiled 2q depth (deepest, ISA on {backend.name}): "
    f"{max(isa_2q_depths)} "
    f"(2q gates: {max(isa_circuit.num_nonlocal_gates() for isa_circuit in isa_circuits)})"
)

estimator = Estimator(backend)
estimator.options.environment.job_tags = ["TUT_SNS"]
estimator.options.dynamical_decoupling.enable = True
estimator.options.dynamical_decoupling.sequence_type = "XY4"
estimator.options.twirling.enable_gates = True
estimator.options.twirling.num_randomizations = 1000
estimator.options.twirling.shots_per_randomization = 128
estimator.options.resilience.measure_mitigation = True
estimator.options.resilience.measure_noise_learning.num_randomizations = 32
estimator.options.resilience.measure_noise_learning.shots_per_randomization = 100

pubs = [
    (
        isa_circuit,
        [obs.apply_layout(isa_circuit.layout) for obs in observables],
    )
    for isa_circuit in isa_circuits
]
job = estimator.run(pubs)
print(f"Job ID: {job.job_id()}")

result = job.result()

# ── Step 4: Post-process ──────────────────────────────────────────────────────
rgf_mat = np.stack([pub_result.data.evs for pub_result in result])

n_points_momentum, n_points_frequency = 100, 100
spectrum = get_dsf(
    n_qubits,
    rgf_mat,
    time_step,
    n_steps,
    n_points_momentum,
    n_points_frequency,
)

plot_rgf(
    n_qubits,
    rgf_mat,
    time_step,
    n_steps,
    title=f"Retarded Green's function — {n_qubits} qubits (QPU)",
)
plot_dsf(
    spectrum,
    time_step,
    n_points_momentum,
    n_points_frequency,
    title=rf"KCuF$_3$ DSF — {n_qubits} qubits (QPU)"
    "\n(AQC + DD + Pauli twirling + TREX)",
)

Ground-state energy (DMRG): -21.972109
Optimizing ground state ansatz...
Finished optimizing ground state ansatz in 132.2076231740648 seconds.
Ground state ansatz fidelity: 0.645956
Ground state ansatz energy: -21.616744
Built 20 circuits, deepest 2q depth (uniform basis) = 307
Defined 50 Z observables.

Step 2b — target MPS:
  k=1: max bond = 44
  k=2: max bond = 46
  k=3: max bond = 53
  k=4: max bond = 62
  k=5: max bond = 75
  k=6: max bond = 96
  k=7: max bond = 118
  k=8: max bond = 128
  k=9: max bond = 128
  k=10: max bond = 128

Step 2c — one-layer ansatz: 2971 parameters, 2q depth (uniform basis) = 39
Step 2c — two-layer ansatz: 3412 parameters, 2q depth (uniform basis) = 45
  k=1: fidelity = 1.0000, 2q depth (uniform basis) = 39, 215.6s
  k=2: fidelity = 0.9676, 2q depth (uniform basis) = 39, 269.8s
  k=3: fidelity = 0.9014, 2q depth (uniform basis) = 39, 293.5s
  k=4: fidelity = 0.8755, 2q depth (uniform basis) = 39, 279.6s
  k=5: fidelity = 0.8450, 2q depth (uniform basis) = 39, 266.4s
  k=6: fidelity = 0.8146, 2q depth (uniform basis) = 39, 371.0s

E0626 02:59:07.229070 1508496 slow_operation_alarm.cc:73] 
********************************
[Compiling module jit_MakeArrayFn for CPU] Very slow compile? If you want to file a bug, run with envvar XLA_FLAGS=--xla_dump_to=/tmp/foo and attach the results.
********************************
E0626 02:59:21.713796 1508477 slow_operation_alarm.cc:140] The operation took 2m14.484932074s

********************************
[Compiling module jit_MakeArrayFn for CPU] Very slow compile? If you want to file a bug, run with envvar XLA_FLAGS=--xla_dump_to=/tmp/foo and attach the results.
********************************

k=7: fidelity = 0.7312, 2q depth (uniform basis) = 45, 10865.3s
  k=8: fidelity = 0.7720, 2q depth (uniform basis) = 45, 1450.3s
  k=9: fidelity = 0.7316, 2q depth (uniform basis) = 45, 3388.6s

E0626 07:21:05.148724 1508496 slow_operation_alarm.cc:73] 
********************************
[Compiling module jit_MakeArrayFn for CPU] Very slow compile? If you want to file a bug, run with envvar XLA_FLAGS=--xla_dump_to=/tmp/foo and attach the results.
********************************
E0626 07:21:24.286095 1508477 slow_operation_alarm.cc:140] The operation took 2m19.137566498s

********************************
[Compiling module jit_MakeArrayFn for CPU] Very slow compile? If you want to file a bug, run with envvar XLA_FLAGS=--xla_dump_to=/tmp/foo and attach the results.
********************************

k=10: fidelity = 0.6726, 2q depth (uniform basis) = 45, 3396.1s

Step 2d — assembled 20 circuits
  At step 20 (uniform basis): Trotter 2q depth = 307, AQC+Trotter 2q depth = 171; 2q gates: Trotter = 3775, AQC+Trotter = 1913

Output of the previous code cell

Backend: ibm_fez
Transpiled 2q depth (deepest, ISA on ibm_fez): 105 (2q gates: 2574)
Job ID: d8v39vhropqc738biotg

Output of the previous code cell

Os resultados de hardware reproduzem as características principais do contínuo de dois spinons: a intensidade de espalhamento se concentra próximo ao vetor de onda antiferromagnético $q = \pi$ em baixa energia e é limitada inferiormente pela dispersão senoidal do spinon, com um amplo contínuo de peso espectral acima dela em vez de um único modo agudo. Essa é a mesma estrutura medida pelo espalhamento inelástico de nêutrons em KCuF $_3$ , e valida o fluxo de trabalho completo de preparação do estado fundamental, perturbação, evolução de Trotter comprimida com AQC e medição com mitigação de erros na escala de 50 qubits.

Próximos passos

Se você achou este trabalho interessante, pode se interessar pelo seguinte material:

Recomendações

Simule o espalhamento de nêutrons com um fluxo de trabalho Serverless de dinâmica AQC + Trotter — um tutorial complementar que executa esse mesmo experimento por meio de um modelo de função Qiskit Serverless implantado
Lee et al., "Benchmarking quantum simulation with neutron-scattering experiments" (arXiv:2603.15608) — o artigo de referência no qual este tutorial se baseia
Técnicas de mitigação e supressão de erros — DD, twirling de Pauli e TREX usados nos experimentos em hardware
Compilação quântica aproximada para circuitos de evolução temporal — tutorial sobre AQC-Tensor
Documentação do AQC-Tensor

Resultados de aprendizagem​

Pré-requisitos​

Contexto​

Espalhamento inelástico de nêutrons e o fator de estrutura dinâmica​

KCuF3_33​: um ímã canônico de líquido de Luttinger​

O que simulamos e medimos​

Fluxo de trabalho de simulação quântica​

Compilação quântica aproximada (AQC)​

Simulação MPS​

Requisitos​

Setup​

Exemplo de simulador em pequena escala​

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

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

Etapa 3: Executar usando primitivas Qiskit​

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

Execução em hardware em larga escala​

Próximos passos​

Resultados de aprendizagem

Pré-requisitos

Contexto

Espalhamento inelástico de nêutrons e o fator de estrutura dinâmica

KCuF $_3$ : um ímã canônico de líquido de Luttinger

O que simulamos e medimos

Fluxo de trabalho de simulação quântica

Compilação quântica aproximada (AQC)

Simulação MPS

Requisitos

Setup

Exemplo de simulador em pequena escala

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

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

Etapa 3: Executar usando primitivas Qiskit

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

Execução em hardware em larga escala

Próximos passos