Symulacja rozpraszania neutronów w materiałach kwantowych za pomocą circuits kwantowych

Szacowany czas użycia: 13 minut na procesorze Heron r2 (UWAGA: To tylko szacunek. Twój czas wykonania może się różnić.)

Którego samouczka powinienem użyć?

Użyj tego samouczka, aby poznać implementację krok po kroku, z obliczeniami klasycznymi działającymi lokalnie na laptopie i wykonaniem sprzętowym na QPU IBM Quantum®. Przykład na dużą skalę wymaga znacznej ilości pamięci, a przybliżona kompilacja kwantowa (AQC) może zająć kilka godzin przy typowych zasobach laptopa. Aby przenieść kompresję na zasoby obliczeniowe i pamięciowe Qiskit Serverless, skorzystaj z towarzyszącego samouczka Serverless.

Efekty uczenia się

Jak widma nieelastycznego rozpraszania neutronów (INS) łączą się z dynamicznymi czynnikami struktury (DSF) kwantowych modeli spinowych.
Jak przygotować stan podstawowy, zastosować lokalne zaburzenie i wykonać ewolucję czasową Trottera na circuit kwantowym.
Jak używać przybliżonej kompilacji kwantowej (AQC) z qiskit-addon-aqc-tensor do kompresji głębokich circuits Trottera pod wykonanie na sprzęcie.
Jak wyodrębnić retardowaną funkcję Greena (RGF) z wartości oczekiwanych qubitów i przekształcić ją Fourierem na DSF.

Wymagania wstępne

Tło

W tym samouczku odtwarzamy wyniki Lee i in., arXiv:2603.15608.

Rozpraszanie neutronów pomaga badaczom charakteryzować materiały istotne dla zrównoważonego rozwoju, w tym elektrody baterii. Ten samouczek analizuje, jak obwody kwantowe mogą symulować wzbudzenia magnetyczne i łączyć modele teoretyczne z pomiarami rozpraszania. Jego przykład magnesu kwantowego rozwija metody badania zachowania materiałów, przyczyniając się do szerszego zestawu narzędzi badawczych do odkrywania materiałów.

Nieelastyczne rozpraszanie neutronów i dynamiczny czynnik struktury

Nieelastyczne rozpraszanie neutronów (INS) jest jedną z najpotężniejszych eksperymentalnych sond wzbudzeń magnetycznych w materiałach kwantowych. Gdy wiązka termicznych lub zimnych neutronów uderza w kryształ, poszczególne neutrony wymieniają zarówno pęd $\mathbf{q}$ , jak i energię $\omega$ z podsystemem magnetycznym. Zmierzone natężenie rozpraszania jest proporcjonalne do dynamicznego czynnika struktury (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,$

który koduje pełne korelacje przestrzenno-czasowe stopni swobody spinu.

KCuF $_3$ : kanoniczny magnes cieczy Luttingera

Fluorek miedzi potasu (KCuF $_3$ ) jest quasi-jednowymiarowym antyferromagnetykiem, w którym łańcuchy jonów Cu $^{2+}$ o spinie $\frac{1}{2}$ oddziałują poprzez wymianę Heisenberga $J$ między najbliższymi sąsiadami, podczas gdy sprzężenie międzyłańcuchowe stanowi zaledwie $\sim 2,7\%$ wartości $J$ . W $T = 6\;\mathrm{K}$ , gdzie dostępne są dane INS, widmo jest zdominowane przez sfrakcjonowane wzbudzenia spinonowe charakterystyczne dla cieczy Tomonagi-Luttingera. Ponieważ dynamika wewnątrzłańcuchowa jest dobrze opisana jednowymiarowym hamiltonianem XXZ o spinie $\frac{1}{2}$ w punkcie izotropowym ( $\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$ służy jako idealny benchmark dla symulacji kwantowej: hamiltonian jest wystarczająco prosty, aby zaimplementować go na procesorze kwantowym, a mimo to stan podstawowy jest silnie splątany, a widmo wzbudzeń zawiera szerokie kontinuum dwuspinonowe.

Uwaga: ten samouczek przyjmuje $J = 1$ jako jednostkę energii i stosuje normalizację $H = J\sum_i[\ldots]$ , która odpowiada lokalnemu hamiltonianowi $H_\mathrm{loc} = J(S^X S^X + S^Y S^Y + S^Z S^Z)$ używanemu w implementacji circuit w artykule (Rys. S3 materiałów uzupełniających). Pełny hamiltonian z artykułu (Rów. 3) zawiera dodatkowy ogólny czynnik 2, więc $J$ z artykułu jest dwukrotnie większe od $J$ używanego tutaj.

Co symulujemy i mierzymy

Wielkością fizyczną, którą obliczamy, jest retardowana funkcja Greena (RGF), zdefiniowana jako zależna od czasu funkcja korelacji spin-spin

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

gdzie $j_c$ jest węzłem referencyjnym (środkiem łańcucha), a $S_j^\alpha(t) = e^{iHt}S_j^\alpha e^{-iHt}$ jest operatorem spinu w obrazie Heisenberga. W tym samouczku skupiamy się na składowej $zz$ ( $\alpha = \beta = z$ ). Na komputerze kwantowym RGF jest dostępna poprzez przygotowanie stanu podstawowego, zastosowanie lokalnego zaburzenia przy $j_c$ , ewolucję czasową zaburzonego stanu i pomiar jednoqubitowej wartości oczekiwanej $\langle\sigma_j^z\rangle$ w każdym węźle $j$ dla każdego kroku czasowego. Kluczowa obserwacja jest taka, że każde $\langle\sigma_j^z\rangle$ daje różnicę względem magnetyzacji stanu podstawowego; ponieważ izotropowy antyferromagnetyk Heisenberga ma zerową netto magnetyzację na węzeł ( $\langle\sigma_j^z\rangle_\mathrm{GS} = 0$ ), surowa zmierzona wartość bezpośrednio daje $G^R(j, j_c, t)$ bez żadnego jawnego odejmowania.

Zbierając $G^R(j, j_c, t)$ dla wszystkich site i kroków czasowych, tworzymy dwuwymiarowy zbiór danych, który jest następnie przekształcany Fourierem zarówno w przestrzeni, jak i w czasie, aby uzyskać dynamiczny czynnik struktury $S(q,\omega)$ . DSF jest wielkością bezpośrednio mierzoną w eksperymencie INS: mówi nam, jakie wzbudzenia magnetyczne istnieją przy każdym pędzie $q$ i energii $\omega$ . Dla izotropowego łańcucha Heisenberga dokładne widmo wzbudzeń to kontinuum dwuspinonowe, szerokie pasmo natężenia rozpraszania, którego kształt służy jako rygorystyczny benchmark end-to-end dla symulacji kwantowej: waliduje jednocześnie przygotowanie stanu podstawowego, zaburzenie, ewolucję czasową Trottera i protokół pomiaru.

Workflow symulacji kwantowej

Workflow odzwierciedla fizykę zdarzenia INS. (1) Przygotowujemy wielociałowy stan podstawowy $|\psi_{\mathrm{GS}}\rangle$ na $n$ qubitach, (2) stosujemy lokalne zaburzenie typu spin-flip $U_{j_c} = \frac{1}{\sqrt{2}}(I - i\sigma^z_{j_c})$ w centrum łańcucha, naśladując transfer spinu neutronu, (3) ewoluujemy pod $H$ przez dyskretne kroki czasowe za pomocą trotteryzacji drugiego rzędu, oraz (4) mierzymy $\langle\sigma_i^z\rangle$ na każdym qubicie w każdym kroku, aby uzyskać RGF. Dwuwymiarowe dyskretne przekształcenie Fouriera daje wówczas DSF $S(q,\omega)$ .

Obserwablą, którą mierzymy w każdym kroku czasowym, jest $\sigma_i^z$ na każdym qubicie $i$ . W Qiskit jest to reprezentowane jako lista operatorów SparsePauliOp: jeden jednoqubitowy $Z$ osadzony w $n$ -qubitowym ciągu tożsamości dla każdego site. Te obserwable są konstruowane raz na instancję problemu podczas kroku mapowania problemu (Krok 1) i ponownie wykorzystywane dla każdego circuit na tej skali.

Przybliżona kompilacja kwantowa (AQC)

Głębokie circuits Trottera mogą być skompresowane za pomocą przybliżonej kompilacji kwantowej (AQC), która zastępuje kilka pierwszych warstw Trottera krótszym sparametryzowanym ansatzem, którego parametry są klasycznie optymalizowane, aby zmaksymalizować wierność na poziomie MPS z oryginalnym głębokim circuit. Pozostałe kroki Trottera są dodawane dokładnie, tworząc "mieszany" circuit AQC + Trotter z znacznie mniejszą liczbą bramek dwuqubitowych.

Symulacja MPS

W przypadku układu jednowymiarowego metody macierzowego stanu produktowego (matrix-product-state, MPS) mogą efektywnie symulować zarówno przygotowanie stanu podstawowego (za pomocą grupy renormalizacji macierzy gęstości, czyli DMRG), jak i ewolucję czasową na poziomie obwodu. Kontrolując wymiar wiązania $\chi$ , wymieniamy dokładność na koszt obliczeniowy. W tym samouczku używamy symulacji MPS w ramach qiskit-addon-aqc-tensor do obliczenia wysokowiernych ansätze AQC, które kompresują głębokie obwody Trottera do wykonania na sprzęcie.

Wymagania

Zanim rozpoczniesz ten samouczek, upewnij się, że masz zainstalowane następujące elementy:

Qiskit SDK z obsługą wizualizacji
Qiskit Runtime (pip install qiskit-ibm-runtime)
qiskit-addon-aqc-tensor z dodatkami quimb i JAX (pip install 'qiskit-addon-aqc-tensor[quimb-jax]')

Konfiguracja

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

Przykład symulatora na małą skalę

Najpierw demonstrujemy pełny przepływ pracy na 10 kubitach, optymalizując ansatz stanu podstawowego HVA za pomocą symulacji MPS i wykorzystując symulator wektora stanu Qiskit do ewolucji czasowej. DMRG dostarcza referencyjną energię stanu podstawowego oraz referencyjny MPS. Ten przykład na małą skalę pozwala nam zweryfikować każdy krok przed skalowaniem.

Krok 1: Odwzorowanie klasycznych danych wejściowych na problem kwantowy

Zaczynamy od zdefiniowania modelu fizycznego i skonstruowania obwodów kwantowych.

Hamiltonian. KCuF $_3$ jest modelowany za pomocą jednowymiarowego hamiltonianu XXZ w punkcie izotropowym ( $J = 1$ , $\epsilon = 1$ , przyjmując $J = 1$ jako jednostkę energii).

Stan podstawowy. Używamy obwodu ansatzu wariacyjnego hamiltonianu (HVA), build_ground_state_ansatz, jako obwodu przygotowania stanu podstawowego. Parametry HVA są optymalizowane klasycznie poprzez maksymalizację wierności stanu $|\langle\psi_{\mathrm{HVA}}(\theta)|\psi_{\mathrm{DMRG}}\rangle|^2$ z referencyjnym stanem podstawowym MPS z DMRG, gdzie stan HVA jest oceniany poprzez symulację obwodu jako macierzowy stan produktowy za pomocą CircuitMPS z quimb. Do optymalizacji używamy scipy.optimize.minimize. Dla porównania obliczana jest również energia $\langle H \rangle$ zoptymalizowanego ansatzu.

Bramki Trottera. Każdy człon oddziaływania między najbliższymi sąsiadami $e^{-i\Delta t\, H_{\mathrm{pair}}}$ z $H_{\mathrm{pair}} = (J/4)(XX + YY + ZZ)$ jest konstruowany za pomocą PauliEvolutionGate(H_pair, time=time_step). Qiskit syntetyzuje to do optymalnej dekompozycji z trzema bramkami CNOT podczas transpilacji.

Perturbacja. Bramka $R_z(\pi/2)$ na centralnym kubicie implementuje $U_{j_c} = \frac{1}{\sqrt{2}}(I - i\sigma^z_{j_c})$ , naśladując lokalne odwrócenie spinu wywołane przez rozproszony neutron.

Obserwable. Konstruujemy obserwablę $\sigma^z$ dla każdego kubitu w danym miejscu. Te obiekty SparsePauliOp są przekazywane do prymitywu Estimator w kroku 3 w celu wyodrębnienia $\langle\sigma_i^z\rangle$ w każdym kroku czasowym.

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

Krok 2: Optymalizacja problemu pod kątem wykonania na sprzęcie kwantowym

Na prawdziwym sprzęcie powyższe obwody Trottera byłyby zbyt głębokie. Przybliżona kompilacja kwantowa (AQC) rozwiązuje ten problem, zastępując pierwsze $k$ warstw Trottera (włącznie z obwodem stanu podstawowego) krótszym, sparametryzowanym ansatzem zoptymalizowanym pod kątem maksymalizacji wierności na poziomie MPS względem oryginalnego głębokiego obwodu. Pozostałe kroki Trottera są dołączane dokładnie, co daje płytszy obwód „AQC + Trotter”.

Aby zrównoważyć ekspresyjność i głębokość obwodu, używane są dwa ansätze: jednowarstwowy ansatz (wygenerowany z pojedynczego kroku Trottera) kompresuje najwcześniejsze kroki czasowe, a głębszy dwuwarstwowy ansatz (wygenerowany z dwóch kroków Trottera) kompresuje kolejne kroki, w których potrzebna jest wyższa wierność.

Przepływ pracy AQC składa się z czterech podkroków:

Budowanie obwodów docelowych — pierwsze $k_1 + k_2$ obwodów z kroku 1 służą bezpośrednio jako cele AQC.
Obliczanie docelowego MPS — symulacja każdego obwodu docelowego jako macierzowego stanu produktowego za pomocą quimb.tensor.CircuitMPS.
Generowanie i optymalizacja ansätze — generate_ansatz_from_circuit tworzy jednowarstwowy i dwuwarstwowy sparametryzowany ansatz; parametry są optymalizowane przy użyciu L-BFGS-B z gradientami przyspieszonymi przez JAX w celu minimalizacji $1 - |\langle\psi_{\mathrm{ansatz}}|\psi_{\mathrm{target}}\rangle|^2$ . Pierwsze $k_1$ kroki wykorzystują ansatz jednowarstwowy, a kolejne $k_2$ kroki wykorzystują ansatz dwuwarstwowy; w ramach każdego etapu każdy krok rozpoczyna się z parametrów zoptymalizowanych w poprzednim kroku, a na granicy etapu parametry są resetowane do wartości domyślnych dla danego etapu.
Składanie obwodów mieszanych — dla kroków czasowych wykraczających poza punkt kontrolny AQC, dołącz dokładne warstwy Trottera do zoptymalizowanego dwuwarstwowego obwodu AQC za pomocą 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

Krok 3: Wykonanie za pomocą prymitywów Qiskit

Symulujemy każdy obwód skompilowany za pomocą AQC, wykorzystując prymityw StatevectorEstimator.

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

Krok 4: Postprzetwarzanie i zwrócenie wyniku w pożądanym formacie klasycznym

Teraz wyodrębniamy wartość oczekiwaną $\langle\sigma_i^z\rangle$ dla każdego kubitu $i$ w każdym kroku czasowym $t$ . Wartości te tworzą macierz retardowanej funkcji Greena $G^R(j, j_c, t)$ . Następnie transformujemy metodą Fouriera RGF na dynamiczny czynnik struktury $S(q, \omega)$ i wykreślamy zarówno RGF, jak i DSF. get_dsf stosuje symetrię lustrzaną i wewnętrznie obcina wartości ujemne.

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

Wykonanie na sprzęcie na dużą skalę

Teraz skalujemy do 50 kubitów. Na tę skalę optymalizacje osiągają niższe wierności niż w przykładzie na małą skalę: wierność ansatzu stanu podstawowego spada do około 0,65, a wierności AQC spadają do około 0,7 w późniejszych punktach kontrolnych. Jest to oczekiwane, a wierności te można poprawić, zwiększając liczbę warstw ansatzu stanu podstawowego (gs_n_layers) lub liczbę iteracji optymalizacji (maxiter) kosztem dodatkowych obliczeń klasycznych. Zauważ też, że dwuwarstwowe wierności AQC wypadają niżej niż jednowarstwowe. Nie jest to regresja: późniejsze kroki czasowe generują więcej splątania i są po prostu trudniejsze do skompresowania, dlatego dla nich stosuje się bardziej ekspresyjny ansatz dwuwarstwowy.

Optymalizacja AQC może również zająć kilka godzin klasycznego czasu obliczeniowego (w pokazanym tutaj uruchomieniu było to około sześciu godzin, z czego większość poświęcono na dwuwarstwowe punkty kontrolne w kroku 2c). Aby skrócić czas rzeczywisty, rozważ uruchomienie tego notebooka na wydajniejszym sprzęcie klasycznym, na przykład w systemie wysokowydajnych obliczeń (HPC). Alternatywnie możesz zmniejszyć skalę problemu (na przykład mniej kubitów lub mniej kroków czasowych); twoje wyniki będą się wtedy różnić od tych pokazanych tutaj.

Poniższy kod przebiega według tej samej czteroetapowej struktury co przykład na małą skalę. Parametry stanu podstawowego HVA są ponownie optymalizowane za pomocą symulacji MPS. Na QPU włączamy dynamiczne rozprzęganie (DD), Pauli twirling oraz twirlowane usuwanie błędów odczytu (TREX) w celu tłumienia i mitygacji błędów. Poniższa tabela podsumowuje, czym eksperyment na dużą skalę różni się od tego na małą skalę:

	Mała skala	Duża skala
Kubity	10	50
Kroki czasowe	10	20
Punkty kontrolne AQC (1-warstwowe + 2-warstwowe)	3 + 2 = 5	6 + 4 = 10
Warstwy ansatzu stanu podstawowego	3	5
Maksymalny wymiar wiązania MPS	32	128
Estimator	`StatevectorEstimator`	QPU z DD, Pauli twirling i TREX

Komunikaty o wolnej kompilacji XLA

Podczas optymalizacji AQC (krok 2c poniżej) możesz zobaczyć komunikaty stderr, takie jak:

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

Są to nieszkodliwe komunikaty diagnostyczne z XLA, kompilatora stojącego za automatycznym różniczkowaniem JAX w qiskit-addon-aqc-tensor. Przy 50 kubitach i wymiarze wiązania MPS równym 128, XLA spędza kilka minut na kompilowaniu funkcji gradientu przy pierwszym jej śledzeniu. Kompilacja kończy się powodzeniem, a wyniki optymalizacji pozostają niezmienione.

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

Wyniki uzyskane na sprzęcie odtwarzają kluczowe cechy kontinuum dwuspinonowego: natężenie rozpraszania koncentruje się w pobliżu antyferromagnetycznego wektora falowego $q = \pi$ przy niskiej energii i jest ograniczone od dołu przez sinusoidalną dyspersję spinonów, z szerokim kontinuum wagi widmowej powyżej niego, a nie pojedynczym ostrym modem. Jest to ta sama struktura zmierzona metodą nieelastycznego rozpraszania neutronów na KCuF $_3$ , co potwierdza poprawność całego przepływu pracy: przygotowania stanu podstawowego, perturbacji, ewolucji Trottera skompresowanej metodą AQC oraz pomiaru z mitygacją błędów w skali 50 kubitów.

Kolejne kroki

Jeśli ta praca Cię zainteresowała, mogą Cię zainteresować następujące materiały:

Rekomendacje

Symulacja rozpraszania neutronów za pomocą bezserwerowego przepływu pracy dynamiki AQC + Trotter — towarzyszący tutorial, który uruchamia ten sam eksperyment poprzez wdrożony szablon funkcji Qiskit Serverless
Lee i in., „Benchmarking quantum simulation with neutron-scattering experiments” (arXiv:2603.15608) — artykuł źródłowy, na którym oparty jest ten tutorial
Techniki mitygacji i tłumienia błędów — DD, Pauli twirling i TREX użyte w eksperymentach sprzętowych
Przybliżona kompilacja kwantowa dla obwodów ewolucji czasowej — tutorial na temat AQC-Tensor
Dokumentacja AQC-Tensor

Efekty uczenia się​

Wymagania wstępne​

Tło​

Nieelastyczne rozpraszanie neutronów i dynamiczny czynnik struktury​

KCuF3_33​: kanoniczny magnes cieczy Luttingera​

Co symulujemy i mierzymy​

Workflow symulacji kwantowej​

Przybliżona kompilacja kwantowa (AQC)​

Symulacja MPS​

Wymagania​

Konfiguracja​

Przykład symulatora na małą skalę​

Krok 1: Odwzorowanie klasycznych danych wejściowych na problem kwantowy​

Krok 2: Optymalizacja problemu pod kątem wykonania na sprzęcie kwantowym​

Krok 3: Wykonanie za pomocą prymitywów Qiskit​

Krok 4: Postprzetwarzanie i zwrócenie wyniku w pożądanym formacie klasycznym​

Wykonanie na sprzęcie na dużą skalę​

Kolejne kroki​

Efekty uczenia się

Wymagania wstępne

Tło

Nieelastyczne rozpraszanie neutronów i dynamiczny czynnik struktury

KCuF $_3$ : kanoniczny magnes cieczy Luttingera

Co symulujemy i mierzymy

Workflow symulacji kwantowej

Przybliżona kompilacja kwantowa (AQC)

Symulacja MPS

Wymagania

Konfiguracja

Przykład symulatora na małą skalę

Krok 1: Odwzorowanie klasycznych danych wejściowych na problem kwantowy

Krok 2: Optymalizacja problemu pod kątem wykonania na sprzęcie kwantowym

Krok 3: Wykonanie za pomocą prymitywów Qiskit

Krok 4: Postprzetwarzanie i zwrócenie wyniku w pożądanym formacie klasycznym

Wykonanie na sprzęcie na dużą skalę

Kolejne kroki