Simula lo scattering di neutroni in materiali quantistici con circuiti quantistici

Stima di utilizzo: 13 minuti su un processore Heron r2 (NOTA: questa è solo una stima. La tua durata effettiva potrebbe variare.)

Quale tutorial dovrei usare?

Usa questo tutorial per apprendere l'implementazione passo dopo passo, con il calcolo classico eseguito localmente sul tuo laptop e l'esecuzione hardware su una QPU IBM Quantum®. L'esempio su larga scala richiede memoria considerevole, e la compilazione quantistica approssimata (AQC) può richiedere diverse ore su risorse tipiche di un laptop. Per delegare la compressione alle risorse di calcolo e memoria di Qiskit Serverless, usa il tutorial Serverless di accompagnamento.

Risultati di apprendimento

Come gli spettri di scattering di neutroni anelastico (INS) si collegano ai fattori di struttura dinamici (DSF) dei modelli di spin quantistici.
Come preparare uno stato fondamentale, applicare una perturbazione locale ed eseguire un'evoluzione temporale di Trotter su un circuito quantistico.
Come usare la compilazione quantistica approssimata (AQC) con qiskit-addon-aqc-tensor per comprimere circuiti di Trotter profondi per l'esecuzione hardware.
Come estrarre la funzione di Green ritardata (RGF) dai valori di aspettazione dei qubit e trasformarla con Fourier in un DSF.

Prerequisiti

Contesto

In questo tutorial, riproduciamo i risultati di Lee et al., arXiv:2603.15608.

Lo scattering di neutroni aiuta i ricercatori a caratterizzare materiali rilevanti per la sostenibilità, inclusi gli elettrodi delle batterie. Questo tutorial esplora come i circuiti quantistici possano simulare eccitazioni magnetiche e collegare modelli teorici a misurazioni di scattering. Il suo esempio di magnete quantistico sviluppa metodi per studiare il comportamento dei materiali, contribuendo al più ampio insieme di strumenti di ricerca per la scoperta di nuovi materiali.

Scattering di neutroni anelastico e fattore di struttura dinamico

Lo scattering di neutroni anelastico (INS) è una delle sonde sperimentali più potenti per le eccitazioni magnetiche nei materiali quantistici. Quando un fascio di neutroni termici o freddi colpisce un cristallo, i singoli neutroni scambiano sia momento $\mathbf{q}$ che energia $\omega$ con il sottosistema magnetico. L'intensità di scattering misurata è proporzionale al fattore di struttura dinamico (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,$

che codifica le correlazioni spazio-temporali complete dei gradi di libertà di spin.

KCuF $_3$ : un magnete liquido di Luttinger canonico

Il fluoruro di potassio e rame (KCuF $_3$ ) è un antiferromagnete quasi unidimensionale in cui catene di ioni Cu $^{2+}$ di spin- $\frac{1}{2}$ interagiscono tramite uno scambio di Heisenberg tra primi vicini $J$ , mentre l'accoppiamento tra catene è solo $\sim 2.7\%$ di $J$ . A $T = 6\;\mathrm{K}$ , dove sono disponibili dati INS, lo spettro è dominato da eccitazioni di spinoni frazionati caratteristiche di un liquido di Tomonaga-Luttinger. Poiché la dinamica intracatena è ben descritta dall'Hamiltoniana XXZ unidimensionale di spin- $\frac{1}{2}$ al punto isotropo ( $\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$ funge da benchmark ideale per la simulazione quantistica: l'Hamiltoniana è abbastanza semplice da implementare su un processore quantistico, ma lo stato fondamentale è fortemente entangled e lo spettro di eccitazione presenta un ampio continuo a due spinoni.

Nota: questo tutorial imposta $J = 1$ come unità di energia e adotta la normalizzazione $H = J\sum_i[\ldots]$ , che corrisponde all'Hamiltoniana locale $H_\mathrm{loc} = J(S^X S^X + S^Y S^Y + S^Z S^Z)$ usata nell'implementazione del circuito dell'articolo (Fig. S3 del materiale supplementare). L'Hamiltoniana completa dell'articolo (Eq. 3) porta un ulteriore fattore complessivo di 2, quindi il $J$ dell'articolo è il doppio del $J$ usato qui.

Cosa simuliamo e misuriamo

La grandezza fisica che calcoliamo è la funzione di Green ritardata (RGF), definita come la funzione di correlazione spin-spin dipendente dal 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,$

dove $j_c$ è un sito di riferimento (il centro della catena) e $S_j^\alpha(t) = e^{iHt}S_j^\alpha e^{-iHt}$ è l'operatore di spin nella rappresentazione di Heisenberg. In questo tutorial ci concentriamo sulla componente $zz$ ( $\alpha = \beta = z$ ). Su un computer quantistico l'RGF si ottiene preparando lo stato fondamentale, applicando una perturbazione locale a $j_c$ , evolvendo temporalmente lo stato perturbato, e misurando il valore di aspettazione a singolo qubit $\langle\sigma_j^z\rangle$ a ogni sito $j$ per ogni passo temporale. L'intuizione chiave è che ogni $\langle\sigma_j^z\rangle$ dà la differenza rispetto alla magnetizzazione dello stato fondamentale; poiché l'antiferromagnete di Heisenberg isotropo ha magnetizzazione netta per sito nulla ( $\langle\sigma_j^z\rangle_\mathrm{GS} = 0$ ), il valore grezzo misurato fornisce direttamente $G^R(j, j_c, t)$ senza alcuna sottrazione esplicita.

Raccogliendo $G^R(j, j_c, t)$ su tutti i siti e passi temporali, assembliamo un dataset bidimensionale che viene poi trasformato con Fourier sia nello spazio che nel tempo per produrre il fattore di struttura dinamico $S(q,\omega)$ . Il DSF è la grandezza misurata direttamente in un esperimento INS: ci dice quali eccitazioni magnetiche esistono a ogni momento $q$ ed energia $\omega$ . Per la catena di Heisenberg isotropa lo spettro di eccitazione esatto è un continuo a due spinoni, un'ampia banda di intensità di scattering la cui forma funge da rigoroso benchmark end-to-end per la simulazione quantistica: convalida contemporaneamente la preparazione dello stato fondamentale, la perturbazione, l'evoluzione temporale di Trotter e il protocollo di misura.

Workflow di simulazione quantistica

Il workflow rispecchia la fisica di un evento INS. Noi (1) prepariamo lo stato fondamentale a molti corpi $|\psi_{\mathrm{GS}}\rangle$ su $n$ qubit, (2) applichiamo una perturbazione locale di spin-flip $U_{j_c} = \frac{1}{\sqrt{2}}(I - i\sigma^z_{j_c})$ al centro della catena per imitare il trasferimento di spin del neutrone, (3) evolviamo sotto $H$ per passi temporali discreti usando la trotterizzazione del secondo ordine, e (4) misuriamo $\langle\sigma_i^z\rangle$ su ogni qubit a ogni passo per ottenere l'RGF. Una trasformata di Fourier discreta bidimensionale produce quindi il DSF $S(q,\omega)$ .

L'osservabile che misuriamo a ogni passo temporale è $\sigma_i^z$ su ogni qubit $i$ . In Qiskit questo è rappresentato come un elenco di operatori SparsePauliOp: uno $Z$ a singolo qubit incorporato nella stringa identità a $n$ qubit per ogni sito. Questi osservabili vengono costruiti una sola volta per istanza del problema durante la fase di mappatura del problema (Passo 1) e riutilizzati per ogni circuito a quella scala.

Compilazione quantistica approssimata (AQC)

I circuiti di Trotter profondi possono essere compressi con la compilazione quantistica approssimata (AQC), che sostituisce i primi livelli di Trotter con un ansatz parametrizzato più breve, i cui parametri vengono ottimizzati classicamente per massimizzare la fedeltà a livello MPS con il circuito profondo originale. I restanti passi Trotter vengono aggiunti esattamente, producendo un circuito "misto" AQC + Trotter con un numero sostanzialmente inferiore di gate a due qubit.

Simulazione MPS

Per un sistema unidimensionale, i metodi a stato-di-prodotto-matriciale (MPS) possono simulare efficientemente sia la preparazione dello stato fondamentale (con il gruppo di rinormalizzazione della matrice densità, o DMRG) sia l'evoluzione temporale a livello di circuito. Controllando la dimensione di legame $\chi$ , bilanciamo l'accuratezza rispetto al costo computazionale. In questo tutorial usiamo la simulazione MPS all'interno di qiskit-addon-aqc-tensor per calcolare ansätze AQC ad alta fedeltà che comprimono i circuiti di Trotter profondi per l'esecuzione hardware.

Requisiti

Prima di iniziare questo tutorial, assicurati di avere installato quanto segue:

Qiskit SDK con supporto per la visualizzazione
Qiskit Runtime (pip install qiskit-ibm-runtime)
qiskit-addon-aqc-tensor con gli extra quimb e JAX (pip install 'qiskit-addon-aqc-tensor[quimb-jax]')

Configurazione

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

Esempio simulatore su piccola scala

Prima dimostriamo l'intero workflow su 10 qubit, ottimizzando l'ansatz HVA dello stato fondamentale usando la simulazione MPS e utilizzando il simulatore statevector di Qiskit per l'evoluzione temporale. Il DMRG fornisce un'energia di riferimento dello stato fondamentale e un MPS di riferimento. Questo esempio su piccola scala ci permette di convalidare ogni passo prima di scalare.

Passo 1: Mappare gli input classici in un problema quantistico

Iniziamo definendo il modello fisico e costruendo i circuiti quantistici.

Hamiltoniana. KCuF $_3$ è modellato dall'Hamiltoniana XXZ 1D al punto isotropo ( $J = 1$ , $\epsilon = 1$ , impostando $J = 1$ come unità di energia).

Stato fondamentale. Usiamo un circuito ansatz variazionale hamiltoniano (HVA), build_ground_state_ansatz, come circuito di preparazione dello stato fondamentale. I parametri HVA vengono ottimizzati classicamente massimizzando la fedeltà dello stato $|\langle\psi_{\mathrm{HVA}}(\theta)|\psi_{\mathrm{DMRG}}\rangle|^2$ con l'MPS dello stato fondamentale DMRG, dove lo stato HVA viene valutato simulando il circuito come stato di prodotto matriciale con CircuitMPS di quimb. Usiamo scipy.optimize.minimize per l'ottimizzazione. Viene calcolata anche l'energia $\langle H \rangle$ dell'ansatz ottimizzato per riferimento.

Gate di Trotter. Ogni termine di interazione tra primi vicini $e^{-i\Delta t\, H_{\mathrm{pair}}}$ con $H_{\mathrm{pair}} = (J/4)(XX + YY + ZZ)$ è costruito con PauliEvolutionGate(H_pair, time=time_step). Qiskit lo sintetizza nella decomposizione ottimale a tre CNOT durante la transpilazione.

Perturbazione. Un gate $R_z(\pi/2)$ sul qubit centrale implementa $U_{j_c} = \frac{1}{\sqrt{2}}(I - i\sigma^z_{j_c})$ , imitando lo spin-flip locale prodotto da un neutrone diffuso.

Osservabili. Costruiamo un osservabile $\sigma^z$ per ogni sito di qubit. Questi oggetti SparsePauliOp vengono passati alla primitiva Estimator nel Passo 3 per estrarre $\langle\sigma_i^z\rangle$ a ogni passo temporale.

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

Passo 2: Ottimizzare il problema per l'esecuzione su hardware quantistico

Su hardware reale, i circuiti di Trotter sopra sarebbero troppo profondi. La compilazione quantistica approssimata (AQC) risolve questo problema sostituendo i primi $k$ livelli di Trotter (incluso il circuito dello stato fondamentale) con un ansatz parametrizzato più breve ottimizzato per massimizzare la fedeltà a livello MPS con il circuito profondo originale. I restanti passi Trotter vengono aggiunti esattamente, producendo un circuito "AQC + Trotter" meno profondo.

Per bilanciare l'espressività rispetto alla profondità del circuito, vengono usati due ansätze: un ansatz a un livello (generato da un singolo passo di Trotter) comprime i primi passi temporali, e un più profondo ansatz a due livelli (generato da due passi di Trotter) comprime i pochi passi successivi dove è necessaria una fedeltà maggiore.

Il workflow AQC ha quattro sotto-passi:

Costruire i circuiti target — i primi $k_1 + k_2$ circuiti dal Passo 1 fungono direttamente da target AQC.
Calcolare l'MPS target — simula ogni circuito target come stato di prodotto matriciale usando quimb.tensor.CircuitMPS.
Generare e ottimizzare gli ansätze — generate_ansatz_from_circuit crea un ansatz a un livello e uno a due livelli; i parametri vengono ottimizzati usando L-BFGS-B con gradienti accelerati da JAX per minimizzare $1 - |\langle\psi_{\mathrm{ansatz}}|\psi_{\mathrm{target}}\rangle|^2$ . I primi $k_1$ passi usano l'ansatz a un livello e i successivi $k_2$ passi usano l'ansatz a due livelli; all'interno di ogni fase, ogni passo parte a caldo dai parametri ottimizzati del passo precedente, e i parametri vengono reimpostati ai valori predefiniti della fase al confine tra le fasi.
Assemblare i circuiti misti — per i passi temporali oltre il checkpoint AQC, aggiungi livelli di Trotter esatti al circuito AQC a due livelli ottimizzato 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

Passo 3: Eseguire usando le primitive Qiskit

Simuliamo ogni circuito compilato AQC usando la primitiva StatevectorEstimator.

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

Passo 4: Post-processare e restituire il risultato nel formato classico desiderato

Ora estraiamo il valore di aspettazione $\langle\sigma_i^z\rangle$ per ogni qubit $i$ a ogni passo temporale $t$ . Questi valori formano la matrice della funzione di Green ritardata $G^R(j, j_c, t)$ . Poi trasformiamo con Fourier l'RGF nel fattore di struttura dinamico $S(q, \omega)$ e tracciamo sia l'RGF che il DSF. get_dsf applica internamente la simmetria a specchio e taglia i valori negativi.

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

Esecuzione hardware su larga scala

Ora scaliamo a 50 qubit. A questa scala, le ottimizzazioni raggiungono fedeltà inferiori rispetto all'esempio su piccola scala: la fedeltà dell'ansatz dello stato fondamentale scende a circa 0,65, e le fedeltà AQC scendono a circa 0,7 nei checkpoint successivi. È un comportamento atteso, e puoi migliorare queste fedeltà aumentando il numero di livelli dell'ansatz dello stato fondamentale (gs_n_layers) o il numero di iterazioni di ottimizzazione (maxiter) a costo classico aggiuntivo. Nota anche che le fedeltà AQC a due livelli risultano inferiori a quelle a un livello. Questa non è una regressione: i passi temporali successivi generano più entanglement e sono semplicemente più difficili da comprimere, motivo per cui l'ansatz a due livelli, più espressivo, viene usato per essi.

L'ottimizzazione AQC può anche richiedere diverse ore di tempo di calcolo classico (circa sei ore nell'esecuzione mostrata qui, la maggior parte spesa sui checkpoint a due livelli nel Passo 2c). Per ridurre il tempo di orologio, considera di eseguire questo notebook su hardware classico più potente, come un sistema di high-performance computing (HPC). In alternativa, puoi ridurre la scala a un'istanza del problema più piccola (ad esempio, meno qubit o meno passi temporali); i tuoi risultati differiranno quindi da quelli mostrati qui.

Il codice seguente segue la stessa struttura a quattro passi dell'esempio su piccola scala. I parametri HVA dello stato fondamentale vengono nuovamente ottimizzati usando la simulazione MPS. Sulla QPU, abilitiamo dynamical decoupling (DD), Pauli twirling e twirled readout error extinction (TREX) per la soppressione e mitigazione degli errori. La tabella seguente riassume come l'esperimento su larga scala differisce da quello su piccola scala:

	Piccola scala	Grande scala
Qubit	10	50
Passi temporali	10	20
Checkpoint AQC (1 livello + 2 livelli)	3 + 2 = 5	6 + 4 = 10
Livelli ansatz stato fondamentale	3	5
Dimensione massima di legame MPS	32	128
Estimator	`StatevectorEstimator`	QPU con DD, Pauli twirling e TREX

Messaggi di compilazione lenta XLA

Durante l'ottimizzazione AQC (Passo 2c seguente) potresti vedere messaggi stderr come:

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

Queste sono diagnostiche benigne di XLA, il compilatore che supporta l'autodiff JAX di qiskit-addon-aqc-tensor. A 50 qubit con dimensione di legame MPS 128, XLA impiega un paio di minuti per compilare la funzione del gradiente la prima volta che viene tracciata. La compilazione ha successo, e i risultati dell'ottimizzazione non sono influenzati.

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

I risultati hardware riproducono le caratteristiche chiave del continuo a due spinoni: l'intensità di scattering è concentrata vicino al vettore d'onda antiferromagnetico $q = \pi$ a bassa energia ed è delimitata dal basso dalla dispersione sinusoidale degli spinoni, con un ampio continuo di peso spettrale sopra di essa piuttosto che un singolo modo netto. Questa è la stessa struttura misurata dallo scattering di neutroni anelastico su KCuF $_3$ , e convalida l'intero workflow di preparazione dello stato fondamentale, perturbazione, evoluzione di Trotter compressa AQC e misura con mitigazione degli errori alla scala di 50 qubit.

Passi successivi

Se hai trovato interessante questo lavoro, potresti essere interessato al seguente materiale:

Raccomandazioni

Simula lo scattering di neutroni con un workflow Serverless di dinamica AQC + Trotter — un tutorial complementare che esegue questo stesso esperimento tramite un modello di funzione Qiskit Serverless distribuito
Lee et al., "Benchmarking quantum simulation with neutron-scattering experiments" (arXiv:2603.15608) — l'articolo di riferimento su cui si basa questo tutorial
Tecniche di mitigazione e soppressione degli errori — DD, Pauli twirling e TREX usati negli esperimenti hardware
Compilazione quantistica approssimata per circuiti di evoluzione temporale — tutorial su AQC-Tensor
Documentazione AQC-Tensor

Risultati di apprendimento​

Prerequisiti​

Contesto​

Scattering di neutroni anelastico e fattore di struttura dinamico​

KCuF3_33​: un magnete liquido di Luttinger canonico​

Cosa simuliamo e misuriamo​

Workflow di simulazione quantistica​

Compilazione quantistica approssimata (AQC)​

Simulazione MPS​

Requisiti​

Configurazione​

Esempio simulatore su piccola scala​

Passo 1: Mappare gli input classici in un problema quantistico​

Passo 2: Ottimizzare il problema per l'esecuzione su hardware quantistico​

Passo 3: Eseguire usando le primitive Qiskit​

Passo 4: Post-processare e restituire il risultato nel formato classico desiderato​

Esecuzione hardware su larga scala​

Passi successivi​

Risultati di apprendimento

Prerequisiti

Contesto

Scattering di neutroni anelastico e fattore di struttura dinamico

KCuF $_3$ : un magnete liquido di Luttinger canonico

Cosa simuliamo e misuriamo

Workflow di simulazione quantistica

Compilazione quantistica approssimata (AQC)

Simulazione MPS

Requisiti

Configurazione

Esempio simulatore su piccola scala

Passo 1: Mappare gli input classici in un problema quantistico

Passo 2: Ottimizzare il problema per l'esecuzione su hardware quantistico

Passo 3: Eseguire usando le primitive Qiskit

Passo 4: Post-processare e restituire il risultato nel formato classico desiderato

Esecuzione hardware su larga scala

Passi successivi