Simular la dispersión de neutrones en materiales cuánticos con circuitos cuánticos

Estimación de uso: 13 minutos en un procesador Heron r2 (NOTA: Esto es solo una estimación. Tu tiempo de ejecución puede variar.)

¿Qué tutorial debo usar?

Usa este tutorial para aprender la implementación paso a paso, con el cómputo clásico ejecutándose localmente en tu portátil y la ejecución en hardware en una QPU de IBM Quantum®. El ejemplo a gran escala requiere memoria sustancial, y la compilación cuántica aproximada (AQC) puede tardar varias horas con recursos típicos de portátil. Para descargar la compresión en los recursos de cómputo y memoria de Qiskit Serverless, usa el tutorial complementario de Serverless.

Resultados de aprendizaje

Cómo se relacionan los espectros de dispersión inelástica de neutrones (INS) con los factores de estructura dinámica (DSF) de los modelos de espín cuántico.
Cómo preparar un estado fundamental, aplicar una perturbación local y realizar una evolución temporal de Trotter en un circuito cuántico.
Cómo usar la compilación cuántica aproximada (AQC) con qiskit-addon-aqc-tensor para comprimir circuitos de Trotter profundos para su ejecución en hardware.
Cómo extraer la función de Green retardada (RGF) a partir de valores esperados de qubits y transformarla mediante Fourier en un DSF.

Requisitos previos

Antecedentes

En este tutorial, reproducimos los resultados de Lee et al., arXiv:2603.15608.

La dispersión de neutrones ayuda a los investigadores a caracterizar materiales relevantes para la sostenibilidad, incluyendo electrodos de baterías. Este tutorial explora cómo los circuitos cuánticos pueden simular excitaciones magnéticas y conectar modelos teóricos con mediciones de dispersión. Su ejemplo de imán cuántico desarrolla métodos para estudiar el comportamiento de los materiales, ampliando las herramientas de investigación para el descubrimiento de materiales.

Dispersión inelástica de neutrones y el factor de estructura dinámica

La dispersión inelástica de neutrones (INS) es una de las sondas experimentales más potentes de las excitaciones magnéticas en materiales cuánticos. Cuando un haz de neutrones térmicos o fríos incide sobre un cristal, los neutrones individuales intercambian tanto momento $\mathbf{q}$ como energía $\omega$ con el subsistema magnético. La intensidad de dispersión medida es proporcional al factor de estructura 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 las correlaciones espacio-temporales completas de los grados de libertad de espín.

KCuF $_3$ : un imán canónico de líquido de Luttinger

El fluoruro de cobre y potasio (KCuF $_3$ ) es un antiferromagneto cuasi unidimensional en el que las cadenas de iones de espín $\frac{1}{2}$ Cu $^{2+}$ interactúan mediante un intercambio de Heisenberg de vecinos más cercanos $J$ , mientras que el acoplamiento entre cadenas es solo $\sim 2.7\%$ de $J$ . A $T = 6\;\mathrm{K}$ , donde hay datos de INS disponibles, el espectro está dominado por excitaciones de espinones fraccionalizados características de un líquido de Tomonaga-Luttinger. Debido a que la dinámica intracadena se describe bien mediante el hamiltoniano XXZ unidimensional de espín $\frac{1}{2}$ en el punto 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$ sirve como un banco de pruebas ideal para la simulación cuántica: el hamiltoniano es lo bastante simple como para implementarlo en un procesador cuántico, pero el estado fundamental está fuertemente entrelazado y el espectro de excitación presenta un amplio continuo de dos espinones.

Nota: este tutorial establece $J = 1$ como unidad de energía y adopta la normalización $H = J\sum_i[\ldots]$ , que corresponde al hamiltoniano local $H_\mathrm{loc} = J(S^X S^X + S^Y S^Y + S^Z S^Z)$ utilizado en la implementación de circuito del artículo (Fig. S3 del material suplementario). El hamiltoniano completo del artículo (Ec. 3) lleva un factor adicional de 2, así que la $J$ del artículo es el doble de la $J$ usada aquí.

Qué simulamos y medimos

La cantidad física que calculamos es la función de Green retardada (RGF), definida como la función de correlación espín-espín dependiente del tiempo

$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,$

donde $j_c$ es un sitio de referencia (el centro de la cadena) y $S_j^\alpha(t) = e^{iHt}S_j^\alpha e^{-iHt}$ es el operador de espín en la representación de Heisenberg. En este tutorial nos centramos en la componente $zz$ ( $\alpha = \beta = z$ ). En un ordenador cuántico, se accede a la RGF preparando el estado fundamental, aplicando una perturbación local en $j_c$ , evolucionando en el tiempo el estado perturbado y midiendo el valor esperado de un solo qubit $\langle\sigma_j^z\rangle$ en cada sitio $j$ para cada paso temporal. La idea clave es que cada $\langle\sigma_j^z\rangle$ da la diferencia respecto a la magnetización del estado fundamental; debido a que el antiferromagneto de Heisenberg isotrópico tiene magnetización neta por sitio nula ( $\langle\sigma_j^z\rangle_\mathrm{GS} = 0$ ), el valor medido bruto produce directamente $G^R(j, j_c, t)$ sin necesidad de ninguna resta explícita.

Al recopilar $G^R(j, j_c, t)$ en todos los sitios y pasos temporales, ensamblamos un conjunto de datos bidimensional que luego se transforma mediante Fourier tanto en el espacio como en el tiempo para producir el factor de estructura dinámica $S(q,\omega)$ . El DSF es la cantidad medida directamente en un experimento de INS: nos indica qué excitaciones magnéticas existen en cada momento $q$ y energía $\omega$ . Para la cadena de Heisenberg isotrópica, el espectro de excitación exacto es un continuo de dos espinones, una amplia banda de intensidad de dispersión cuya forma sirve como un exigente banco de pruebas de extremo a extremo para la simulación cuántica: valida a la vez la preparación del estado fundamental, la perturbación, la evolución temporal de Trotter y el protocolo de medición.

Flujo de trabajo de simulación cuántica

El flujo de trabajo refleja la física de un evento de INS. (1) Preparamos el estado fundamental de muchos cuerpos $|\psi_{\mathrm{GS}}\rangle$ en $n$ qubits, (2) aplicamos una perturbación local de inversión de espín $U_{j_c} = \frac{1}{\sqrt{2}}(I - i\sigma^z_{j_c})$ en el centro de la cadena para imitar la transferencia de espín del neutrón, (3) evolucionamos bajo $H$ durante pasos temporales discretos usando trotterización de segundo orden, y (4) medimos $\langle\sigma_i^z\rangle$ en cada qubit en cada paso para obtener la RGF. Una transformada de Fourier discreta bidimensional produce entonces el DSF $S(q,\omega)$ .

El observable que medimos en cada paso temporal es $\sigma_i^z$ en cada qubit $i$ . En Qiskit, esto se representa como una lista de operadores SparsePauliOp: un único $Z$ de un qubit incrustado en la cadena de identidad de $n$ qubits para cada sitio. Estos observables se construyen una sola vez por instancia del problema durante el paso de mapeo del problema (Paso 1) y se reutilizan para cada circuito a esa escala.

Compilación cuántica aproximada (AQC)

Los circuitos de Trotter profundos se pueden comprimir con compilación cuántica aproximada (AQC), que reemplaza las primeras capas de Trotter con un ansatz parametrizado más corto cuyos parámetros se optimizan clásicamente para maximizar la fidelidad a nivel MPS con el circuito profundo original. Los pasos de Trotter restantes se añaden de forma exacta, produciendo un circuito "mixto" AQC + Trotter con sustancialmente menos compuertas de dos qubits.

Simulación MPS

Para un sistema unidimensional, los métodos de estado de producto matricial (MPS) pueden simular eficientemente tanto la preparación del estado fundamental (con el grupo de renormalización de la matriz de densidad, o DMRG) como la evolución temporal a nivel de circuito. Controlando la dimensión de enlace $\chi$ , equilibramos precisión y coste computacional. En este tutorial usamos la simulación MPS dentro de qiskit-addon-aqc-tensor para calcular ansätze AQC de alta fidelidad que comprimen los circuitos de Trotter profundos para su ejecución en hardware.

Requisitos

Antes de empezar este tutorial, asegúrate de tener instalado lo siguiente:

Qiskit SDK con soporte de visualización
Qiskit Runtime (pip install qiskit-ibm-runtime)
qiskit-addon-aqc-tensor con los extras quimb y JAX (pip install 'qiskit-addon-aqc-tensor[quimb-jax]')

Configuración

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

Ejemplo de simulador a pequeña escala

Primero demostramos el flujo de trabajo completo en 10 qubits, optimizando el ansatz del estado fundamental HVA mediante simulación MPS y usando el simulador de vector de estado de Qiskit para la evolución temporal. DMRG proporciona una energía de estado fundamental de referencia y una MPS de referencia. Este ejemplo a pequeña escala nos permite validar cada paso antes de escalar.

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

Comenzamos definiendo el modelo físico y construyendo los circuitos cuánticos.

Hamiltoniano. KCuF $_3$ se modela mediante el hamiltoniano XXZ 1D en el punto isotrópico ( $J = 1$ , $\epsilon = 1$ , estableciendo $J = 1$ como unidad de energía).

Estado fundamental. Usamos un circuito de ansatz variacional hamiltoniano (HVA), build_ground_state_ansatz, como circuito de preparación del estado fundamental. Los parámetros de HVA se optimizan clásicamente maximizando la fidelidad del estado $|\langle\psi_{\mathrm{HVA}}(\theta)|\psi_{\mathrm{DMRG}}\rangle|^2$ con la MPS del estado fundamental de DMRG, donde el estado HVA se evalúa simulando el circuito como un estado de producto matricial con CircuitMPS de quimb. Usamos scipy.optimize.minimize para la optimización. También se calcula, como referencia, la energía $\langle H \rangle$ del ansatz optimizado.

Compuertas de Trotter. Cada término de interacción de vecinos más cercanos $e^{-i\Delta t\, H_{\mathrm{pair}}}$ con $H_{\mathrm{pair}} = (J/4)(XX + YY + ZZ)$ se construye con PauliEvolutionGate(H_pair, time=time_step). Qiskit sintetiza esto en la descomposición óptima de tres CNOT durante la transpilación.

Perturbación. Una compuerta $R_z(\pi/2)$ en el qubit central implementa $U_{j_c} = \frac{1}{\sqrt{2}}(I - i\sigma^z_{j_c})$ , imitando la inversión de espín local producida por un neutrón dispersado.

Observables. Construimos un observable $\sigma^z$ para cada sitio de qubit. Estos objetos SparsePauliOp se pasan a la primitiva Estimator en el Paso 3 para extraer $\langle\sigma_i^z\rangle$ en cada paso temporal.

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

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

Para hardware real, los circuitos de Trotter anteriores serían demasiado profundos. La compilación cuántica aproximada (AQC) aborda esto reemplazando las primeras $k$ capas de Trotter (incluido el circuito del estado fundamental) con un ansatz parametrizado más corto optimizado para maximizar la fidelidad a nivel MPS con el circuito profundo original. Los pasos de Trotter restantes se añaden de forma exacta, produciendo un circuito "AQC + Trotter" más superficial.

Para equilibrar la expresividad frente a la profundidad del circuito, se usan dos ansätze: un ansatz de una capa (generado a partir de un único paso de Trotter) comprime los primeros pasos temporales, y un ansatz de dos capas más profundo (generado a partir de dos pasos de Trotter) comprime los siguientes pasos donde se necesita mayor fidelidad.

El flujo de trabajo AQC tiene cuatro subpasos:

Construir circuitos objetivo — los primeros $k_1 + k_2$ circuitos del Paso 1 sirven directamente como objetivos de AQC.
Calcular la MPS objetivo — simula cada circuito objetivo como un estado de producto matricial usando quimb.tensor.CircuitMPS.
Generar y optimizar ansätze — generate_ansatz_from_circuit crea un ansatz de una capa y otro de dos capas parametrizados; los parámetros se optimizan usando L-BFGS-B con gradientes acelerados por JAX para minimizar $1 - |\langle\psi_{\mathrm{ansatz}}|\psi_{\mathrm{target}}\rangle|^2$ . Los primeros $k_1$ pasos usan el ansatz de una capa y los siguientes $k_2$ pasos usan el ansatz de dos capas; dentro de cada etapa, cada paso arranca en caliente a partir de los parámetros optimizados del paso anterior, y los parámetros se reinician a los valores predeterminados de la etapa en el límite de esta.
Ensamblar circuitos mixtos — para los pasos temporales más allá del punto de control de AQC, añade capas de Trotter exactas al circuito AQC de dos capas optimizado 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

Paso 3: Ejecutar usando primitivas de Qiskit

Simulamos cada circuito compilado con AQC usando la primitiva StatevectorEstimator.

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

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

Ahora extraemos el valor esperado $\langle\sigma_i^z\rangle$ para cada qubit $i$ en cada paso temporal $t$ . Estos valores forman la matriz de la función de Green retardada $G^R(j, j_c, t)$ . Luego, transformamos mediante Fourier la RGF en el factor de estructura dinámica $S(q, \omega)$ y graficamos tanto la RGF como el DSF. get_dsf aplica simetría de espejo y recorta los 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

Ejecución en hardware a gran escala

Ahora escalamos a 50 qubits. A esta escala, las optimizaciones alcanzan fidelidades más bajas que en el ejemplo a pequeña escala: la fidelidad del ansatz del estado fundamental cae a aproximadamente 0.65, y las fidelidades de AQC disminuyen a aproximadamente 0.7 en los últimos puntos de control. Esto es lo esperado, y puedes mejorar estas fidelidades aumentando el número de capas del ansatz del estado fundamental (gs_n_layers) o el número de iteraciones de optimización (maxiter), con un coste clásico adicional. Ten en cuenta también que las fidelidades de AQC de dos capas resultan más bajas que las de una capa. Esto no es una regresión: los pasos temporales posteriores generan más entrelazamiento y son simplemente más difíciles de comprimir, por lo que se usa el ansatz de dos capas, más expresivo, para ellos.

La optimización de AQC también puede tardar varias horas de tiempo de cómputo clásico (aproximadamente seis horas en la ejecución mostrada aquí, la mayor parte dedicada a los puntos de control de dos capas en el Paso 2c). Para reducir el tiempo de reloj, considera ejecutar este notebook en hardware clásico más potente, como un sistema de computación de alto rendimiento (HPC). Alternativamente, puedes reducir a una instancia de problema más pequeña (por ejemplo, menos qubits o menos pasos temporales); tus resultados diferirán entonces de los mostrados aquí.

El código a continuación sigue la misma estructura de cuatro pasos que el ejemplo a pequeña escala. Los parámetros del estado fundamental HVA se optimizan de nuevo usando simulación MPS. En la QPU, activamos el desacoplamiento dinámico (DD), el twirling de Pauli y la extinción de errores de lectura mediante twirling (TREX) para la supresión y mitigación de errores. La siguiente tabla resume en qué se diferencia el experimento a gran escala del de pequeña escala:

	Escala pequeña	Escala grande
Qubits	10	50
Pasos temporales	10	20
Puntos de control de AQC (1 capa + 2 capas)	3 + 2 = 5	6 + 4 = 10
Capas del ansatz del estado fundamental	3	5
Dimensión de enlace máxima de MPS	32	128
Estimator	`StatevectorEstimator`	QPU con DD, twirling de Pauli y TREX

Mensajes de compilación lenta de XLA

Durante la optimización de AQC (Paso 2c más abajo), podrías ver mensajes de stderr como los siguientes:

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

Estos son diagnósticos benignos de XLA, el compilador que respalda la autodiferenciación de JAX de qiskit-addon-aqc-tensor. Con 50 qubits y una dimensión de enlace MPS de 128, XLA tarda un par de minutos en compilar la función de gradiente la primera vez que se traza. La compilación se completa con éxito y los resultados de la optimización no se ven afectados.

# ── 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

Los resultados de hardware reproducen las características clave del continuo de dos espinones: la intensidad de dispersión se concentra cerca del vector de onda antiferromagnético $q = \pi$ a baja energía y está acotada por debajo por la dispersión sinusoidal de espinones, con un amplio continuo de peso espectral por encima en lugar de un único modo agudo. Esta es la misma estructura medida mediante dispersión inelástica de neutrones en KCuF $_3$ , y valida todo el flujo de trabajo de preparación del estado fundamental, perturbación, evolución de Trotter comprimida con AQC y medición con mitigación de errores a la escala de 50 qubits.

Próximos pasos

Si te ha parecido interesante este trabajo, quizá te interese el siguiente material:

Recomendaciones

Simular la dispersión de neutrones con un flujo de trabajo Serverless de dinámica AQC + Trotter — un tutorial complementario que ejecuta este mismo experimento a través de una plantilla de función de Qiskit Serverless desplegada
Lee et al., "Benchmarking quantum simulation with neutron-scattering experiments" (arXiv:2603.15608) — el artículo de referencia en el que se basa este tutorial
Técnicas de mitigación y supresión de errores — DD, twirling de Pauli y TREX usados en los experimentos de hardware
Compilación cuántica aproximada para circuitos de evolución temporal — tutorial sobre AQC-Tensor
Documentación de AQC-Tensor

Resultados de aprendizaje​

Requisitos previos​

Antecedentes​

Dispersión inelástica de neutrones y el factor de estructura dinámica​

KCuF3_33​: un imán canónico de líquido de Luttinger​

Qué simulamos y medimos​

Flujo de trabajo de simulación cuántica​

Compilación cuántica aproximada (AQC)​

Simulación MPS​

Requisitos​

Configuración​

Ejemplo de simulador a pequeña escala​

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

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

Paso 3: Ejecutar usando primitivas de Qiskit​

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

Ejecución en hardware a gran escala​

Próximos pasos​

Resultados de aprendizaje

Requisitos previos

Antecedentes

Dispersión inelástica de neutrones y el factor de estructura dinámica

KCuF $_3$ : un imán canónico de líquido de Luttinger

Qué simulamos y medimos

Flujo de trabajo de simulación cuántica

Compilación cuántica aproximada (AQC)

Simulación MPS

Requisitos

Configuración

Ejemplo de simulador a pequeña escala

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

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

Paso 3: Ejecutar usando primitivas de Qiskit

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

Ejecución en hardware a gran escala

Próximos pasos