Simuliere Neutronenstreuung in Quantenmaterialien mit Quanten-Circuits

Nutzungsschätzung: 13 Minuten auf einem Heron-r2-Prozessor (HINWEIS: Dies ist nur eine Schätzung. Deine tatsächliche Laufzeit kann abweichen.)

Welches Tutorial sollte ich verwenden?

Verwende dieses Tutorial, um die Implementierung Schritt für Schritt zu lernen, wobei die klassische Berechnung lokal auf deinem Laptop läuft und die Hardware-Ausführung auf einer IBM Quantum®-QPU erfolgt. Das großskalige Beispiel benötigt viel Arbeitsspeicher, und approximatives Quantenkompilieren (AQC) kann auf einem typischen Laptop mehrere Stunden dauern. Um die Kompression auf die Rechen- und Speicherressourcen von Qiskit Serverless auszulagern, verwende das begleitende Serverless-Tutorial.

Lernziele

Wie inelastische Neutronenstreuungsspektren (INS) mit dynamischen Strukturfaktoren (DSFs) von Quantenspinmodellen zusammenhängen.
Wie du einen Grundzustand vorbereitest, eine lokale Störung anwendest und eine Trotter-Zeitentwicklung auf einem Quanten-Circuit durchführst.
Wie du approximatives Quantenkompilieren (AQC) mit qiskit-addon-aqc-tensor verwendest, um tiefe Trotter-Circuits für die Hardwareausführung zu komprimieren.
Wie du die retardierte Greensche Funktion (RGF) aus Qubit-Erwartungswerten extrahierst und mittels Fourier-Transformation in einen DSF umwandelst.

Voraussetzungen

Hintergrund

In diesem Tutorial reproduzieren wir die Ergebnisse von Lee et al., arXiv:2603.15608.

Neutronenstreuung hilft Forschenden dabei, Materialien zu charakterisieren, die für Nachhaltigkeit relevant sind, einschließlich Batterieelektroden. Dieses Tutorial untersucht, wie Quantum Circuits magnetische Anregungen simulieren und theoretische Modelle mit Streumessungen verbinden können. Sein Quantenmagnet-Beispiel entwickelt Methoden zur Untersuchung des Materialverhaltens und trägt so zum breiteren Forschungswerkzeugkasten für die Materialentdeckung bei.

Inelastische Neutronenstreuung und der dynamische Strukturfaktor

Inelastische Neutronenstreuung (INS) ist eine der leistungsfähigsten experimentellen Sonden für magnetische Anregungen in Quantenmaterialien. Wenn ein Strahl thermischer oder kalter Neutronen auf einen Kristall trifft, tauschen einzelne Neutronen sowohl Impuls $\mathbf{q}$ als auch Energie $\omega$ mit dem magnetischen Subsystem aus. Die gemessene Streuintensität ist proportional zum dynamischen Strukturfaktor (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,$

der die vollständigen Raum-Zeit-Korrelationen der Spin-Freiheitsgrade kodiert.

KCuF $_3$ : ein kanonischer Luttinger-Flüssigkeits-Magnet

Kaliumkupferfluorid (KCuF $_3$ ) ist ein quasi-eindimensionaler Antiferromagnet, in dem Ketten von Spin- $\frac{1}{2}$ -Cu $^{2+}$ -Ionen über einen Nächste-Nachbar-Heisenberg-Austausch $J$ wechselwirken, während die Kopplung zwischen den Ketten nur $\sim 2.7\%$ von $J$ beträgt. Bei $T = 6\;\mathrm{K}$ , wo INS-Daten vorliegen, wird das Spektrum von fraktionalisierten Spinon-Anregungen dominiert, die charakteristisch für eine Tomonaga-Luttinger-Flüssigkeit sind. Da die Intrachain-Dynamik gut durch den eindimensionalen Spin- $\frac{1}{2}$ -XXZ-Hamiltonian am isotropen Punkt ( $\epsilon = 1$ ) beschrieben wird,

$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$ dient als idealer Benchmark für Quantensimulation: Der Hamiltonian ist einfach genug, um auf einem Quantenprozessor implementiert zu werden, dennoch ist der Grundzustand stark verschränkt und das Anregungsspektrum weist ein breites Zwei-Spinon-Kontinuum auf.

Hinweis: Dieses Tutorial setzt $J = 1$ als Energieeinheit und übernimmt die Normierung $H = J\sum_i[\ldots]$ , was dem lokalen Hamiltonian $H_\mathrm{loc} = J(S^X S^X + S^Y S^Y + S^Z S^Z)$ entspricht, der in der Circuit-Implementierung des Papers verwendet wird (Abb. S3 der Zusatzinformationen). Der vollständige Hamiltonian des Papers (Gl. 3) trägt einen zusätzlichen Gesamtfaktor 2, sodass das $J$ des Papers doppelt so groß ist wie das hier verwendete $J$ .

Was wir simulieren und messen

Die physikalische Größe, die wir berechnen, ist die retardierte Greensche Funktion (RGF), definiert als die zeitabhängige Spin-Spin-Korrelationsfunktion

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

Dabei ist $j_c$ eine Referenzstelle (das Kettenzentrum) und $S_j^\alpha(t) = e^{iHt}S_j^\alpha e^{-iHt}$ der Spinoperator im Heisenberg-Bild. In diesem Tutorial konzentrieren wir uns auf die $zz$ -Komponente ( $\alpha = \beta = z$ ). Auf einem Quantencomputer wird auf die RGF zugegriffen, indem der Grundzustand vorbereitet, eine lokale Störung bei $j_c$ angewendet, der gestörte Zustand zeitlich entwickelt und der Einzel-Qubit-Erwartungswert $\langle\sigma_j^z\rangle$ an jeder Stelle $j$ für jeden Zeitschritt gemessen wird. Die entscheidende Erkenntnis ist, dass jedes $\langle\sigma_j^z\rangle$ die Abweichung von der Grundzustandsmagnetisierung liefert; da der isotrope Heisenberg-Antiferromagnet eine verschwindende Netto-Magnetisierung pro Stelle hat ( $\langle\sigma_j^z\rangle_\mathrm{GS} = 0$ ), ergibt der rohe Messwert direkt $G^R(j, j_c, t)$ ohne explizite Subtraktion.

Indem wir $G^R(j, j_c, t)$ über alle Stellen und Zeitschritte hinweg sammeln, stellen wir einen zweidimensionalen Datensatz zusammen, der dann sowohl räumlich als auch zeitlich per Fourier-Transformation in den dynamischen Strukturfaktor $S(q,\omega)$ umgewandelt wird. Der DSF ist die Größe, die in einem INS-Experiment direkt gemessen wird: Er zeigt uns, welche magnetischen Anregungen bei jedem Impuls $q$ und jeder Energie $\omega$ existieren. Für die isotrope Heisenberg-Kette ist das exakte Anregungsspektrum ein Zwei-Spinon-Kontinuum, ein breites Band von Streuintensität, dessen Form als strenger Ende-zu-Ende-Benchmark für die Quantensimulation dient: Es validiert gleichzeitig die Grundzustandsvorbereitung, die Störung, die Trotter-Zeitentwicklung und das Messprotokoll.

Quantensimulations-Workflow

Der Workflow spiegelt die Physik eines INS-Ereignisses wider. Wir (1) bereiten den Vielteilchen-Grundzustand $|\psi_{\mathrm{GS}}\rangle$ auf $n$ Qubits vor, (2) wenden eine lokale Spin-Flip-Störung $U_{j_c} = \frac{1}{\sqrt{2}}(I - i\sigma^z_{j_c})$ am Kettenzentrum an, um den Spin-Transfer des Neutrons nachzubilden, (3) entwickeln unter $H$ für diskrete Zeitschritte mittels Trotterisierung zweiter Ordnung, und (4) messen $\langle\sigma_i^z\rangle$ an jedem Qubit bei jedem Schritt, um die RGF zu erhalten. Eine zweidimensionale diskrete Fourier-Transformation liefert dann den DSF $S(q,\omega)$ .

Die Observable, die wir bei jedem Zeitschritt messen, ist $\sigma_i^z$ auf jedem Qubit $i$ . In Qiskit wird dies als Liste von SparsePauliOp-Operatoren dargestellt: ein Einzel-Qubit- $Z$ , eingebettet in die $n$ -Qubit-Identitätskette für jede Stelle. Diese Observablen werden einmal pro Probleminstanz während des Problem-Mapping-Schritts (Schritt 1) konstruiert und für jeden Circuit auf dieser Skala wiederverwendet.

Approximatives Quantenkompilieren (AQC)

Tiefe Trotter-Circuits können mit approximativem Quantenkompilieren (AQC) komprimiert werden, wobei die ersten Trotter-Schichten durch einen kürzeren parametrisierten Ansatz ersetzt werden, dessen Parameter klassisch optimiert werden, um die Fidelity auf MPS-Ebene mit dem ursprünglichen tiefen Circuit zu maximieren. Die verbleibenden Trotter-Schritte werden exakt angehängt, wodurch ein „gemischter“ AQC + Trotter-Circuit mit deutlich weniger Zwei-Qubit-Gates entsteht.

MPS-Simulation

Für ein eindimensionales System können Matrix-Produkt-Zustand-Methoden (MPS) sowohl die Grundzustandsvorbereitung (mit der Dichtematrix-Renormierungsgruppe, DMRG) als auch die Zeitentwicklung auf Circuit-Ebene effizient simulieren. Durch die Steuerung der Bindungsdimension $\chi$ wägen wir Genauigkeit gegen Rechenaufwand ab. In diesem Tutorial verwenden wir MPS-Simulation innerhalb von qiskit-addon-aqc-tensor, um hochpräzise AQC-Ansätze zu berechnen, die die tiefen Trotter-Circuits für die Hardware-Ausführung komprimieren.

Anforderungen

Stelle vor Beginn dieses Tutorials sicher, dass Folgendes installiert ist:

Qiskit SDK mit Visualisierungs-Unterstützung
Qiskit Runtime (pip install qiskit-ibm-runtime)
qiskit-addon-aqc-tensor mit quimb- und JAX-Extras (pip install 'qiskit-addon-aqc-tensor[quimb-jax]')

Setup

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    Returns:
        The Scipy OptimizeResult.
    """

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Kleinskaliges Simulator-Beispiel

Wir demonstrieren den vollständigen Workflow zunächst auf 10 Qubits, wobei wir den HVA-Grundzustandsansatz mittels MPS-Simulation optimieren und für die Zeitentwicklung den Qiskit-Statevector-Simulator verwenden. DMRG liefert eine Referenz-Grundzustandsenergie und einen Referenz-MPS. Dieses kleinskalige Beispiel erlaubt es uns, jeden Schritt zu validieren, bevor wir hochskalieren.

Schritt 1: Klassische Eingaben auf ein Quantenproblem abbilden

Wir beginnen damit, das physikalische Modell zu definieren und die Quanten-Circuits zu konstruieren.

Hamiltonian. KCuF $_3$ wird durch den 1D-XXZ-Hamiltonian am isotropen Punkt modelliert ( $J = 1$ , $\epsilon = 1$ , wobei $J = 1$ als Energieeinheit gesetzt wird).

Grundzustand. Wir verwenden einen Hamiltonian-Variational-Ansatz(HVA)-Circuit, build_ground_state_ansatz, als Grundzustandsvorbereitungs-Circuit. Die HVA-Parameter werden klassisch optimiert, indem die Zustandsfidelity $|\langle\psi_{\mathrm{HVA}}(\theta)|\psi_{\mathrm{DMRG}}\rangle|^2$ mit dem DMRG-Grundzustands-MPS maximiert wird, wobei der HVA-Zustand ausgewertet wird, indem der Circuit mit quimbs CircuitMPS als Matrix-Produkt-Zustand simuliert wird. Wir verwenden scipy.optimize.minimize für die Optimierung. Die Energie $\langle H \rangle$ des optimierten Ansatzes wird ebenfalls zu Referenzzwecken berechnet.

Trotter-Gates. Jeder Wechselwirkungsterm zwischen nächsten Nachbarn $e^{-i\Delta t\, H_{\mathrm{pair}}}$ mit $H_{\mathrm{pair}} = (J/4)(XX + YY + ZZ)$ wird mit PauliEvolutionGate(H_pair, time=time_step) konstruiert. Qiskit synthetisiert dies während der Transpilation in die optimale Drei-CNOT-Zerlegung.

Störung. Ein $R_z(\pi/2)$ -Gate auf dem zentralen Qubit implementiert $U_{j_c} = \frac{1}{\sqrt{2}}(I - i\sigma^z_{j_c})$ und ahmt den lokalen Spin-Flip nach, der durch ein gestreutes Neutron erzeugt wird.

Observablen. Wir konstruieren für jede Qubit-Stelle eine $\sigma^z$ -Observable. Diese SparsePauliOp-Objekte werden in Schritt 3 an die Estimator-Primitive übergeben, um $\langle\sigma_i^z\rangle$ zu jedem Zeitschritt zu extrahieren.

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

Schritt 2: Problem für die Ausführung auf Quantenhardware optimieren

Für reale Hardware wären die obigen Trotter-Circuits zu tief. Approximatives Quantenkompilieren (AQC) löst dieses Problem, indem die ersten $k$ Trotter-Schichten (einschließlich des Grundzustands-Circuits) durch einen kürzeren parametrisierten Ansatz ersetzt werden, der so optimiert ist, dass die Fidelity auf MPS-Ebene mit dem ursprünglichen tiefen Circuit maximiert wird. Die verbleibenden Trotter-Schritte werden exakt angehängt, wodurch ein flacherer „AQC + Trotter“-Circuit entsteht.

Um Ausdrucksstärke gegen Circuit-Tiefe abzuwägen, werden zwei Ansätze verwendet: Ein einschichtiger Ansatz (erzeugt aus einem einzigen Trotter-Schritt) komprimiert die frühesten Zeitschritte, und ein tieferer zweischichtiger Ansatz (erzeugt aus zwei Trotter-Schritten) komprimiert die nächsten Schritte, bei denen eine höhere Fidelity benötigt wird.

Der AQC-Workflow besteht aus vier Teilschritten:

Ziel-Circuits erstellen — die ersten $k_1 + k_2$ Circuits aus Schritt 1 dienen direkt als AQC-Ziele.
Ziel-MPS berechnen — simuliere jeden Ziel-Circuit als Matrixproduktzustand mit quimb.tensor.CircuitMPS.
Ansätze erzeugen und optimieren — generate_ansatz_from_circuit erstellt einen einschichtigen und einen zweischichtigen parametrisierten Ansatz; die Parameter werden mit L-BFGS-B unter Verwendung von JAX-beschleunigten Gradienten optimiert, um $1 - |\langle\psi_{\mathrm{ansatz}}|\psi_{\mathrm{target}}\rangle|^2$ zu minimieren. Die ersten $k_1$ Schritte verwenden den einschichtigen Ansatz und die nächsten $k_2$ Schritte den zweischichtigen Ansatz; innerhalb jeder Stufe startet jeder Schritt mit den optimierten Parametern des vorherigen Schritts (Warmstart), und die Parameter werden an der Stufengrenze auf die Standardwerte der Stufe zurückgesetzt.
Gemischte Circuits zusammenstellen — für Zeitschritte jenseits des AQC-Checkpoints werden mit trotter_evolution exakte Trotter-Schichten an den optimierten zweischichtigen AQC-Circuit angehängt.

# 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

Ausgabe der vorherigen Code-Zelle

Schritt 3: Ausführung mit Qiskit-Primitiven

Wir simulieren jeden AQC-kompilierten Circuit mit der StatevectorEstimator-Primitive.

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

Schritt 4: Nachbearbeitung und Rückgabe des Ergebnisses im gewünschten klassischen Format

Wir extrahieren nun den Erwartungswert $\langle\sigma_i^z\rangle$ für jedes Qubit $i$ zu jedem Zeitschritt $t$ . Diese Werte bilden die Matrix der retardierten Greenschen Funktion $G^R(j, j_c, t)$ . Anschließend transformieren wir die RGF per Fourier-Transformation in den dynamischen Strukturfaktor $S(q, \omega)$ und stellen sowohl die RGF als auch den DSF grafisch dar. get_dsf wendet intern Spiegelsymmetrie an und schneidet negative Werte bei null ab.

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)",
)

Ausgabe der vorherigen Code-Zelle

Großskalige Hardware-Ausführung

Wir skalieren nun auf 50 Qubits. Bei diesem Umfang erreichen die Optimierungen niedrigere Fidelities als im kleinskaligen Beispiel: Die Fidelity des Grundzustands-Ansatzes sinkt auf etwa 0,65, und die AQC-Fidelities sinken bei den späteren Checkpoints auf etwa 0,7. Das ist zu erwarten, und du kannst diese Fidelities verbessern, indem du die Anzahl der Grundzustands-Ansatz-Schichten (gs_n_layers) oder die Anzahl der Optimierungsiterationen (maxiter) erhöhst, was zusätzliche klassische Rechenkosten verursacht. Beachte außerdem, dass die zweischichtigen AQC-Fidelities niedriger ausfallen als die einschichtigen. Das ist keine Verschlechterung: Spätere Zeitschritte erzeugen mehr Verschränkung und sind daher einfach schwerer zu komprimieren, weshalb für sie der ausdrucksstärkere zweischichtige Ansatz verwendet wird.

Die AQC-Optimierung kann außerdem mehrere Stunden klassische Rechenzeit in Anspruch nehmen (rund sechs Stunden im hier gezeigten Lauf, wobei der Großteil davon auf die zweischichtigen Checkpoints in Schritt 2c entfällt). Um die tatsächliche Laufzeit zu verkürzen, kannst du dieses Notebook auf leistungsfähigerer klassischer Hardware ausführen, etwa auf einem Hochleistungsrechensystem (HPC). Alternativ kannst du die Problemgröße verkleinern (zum Beispiel weniger Qubits oder weniger Zeitschritte); deine Ergebnisse werden sich dann von den hier gezeigten unterscheiden.

Der folgende Code folgt derselben vierschrittigen Struktur wie das kleinskalige Beispiel. Die HVA-Grundzustandsparameter werden erneut mithilfe von MPS-Simulation optimiert. Auf der QPU aktivieren wir Dynamical Decoupling (DD), Pauli-Twirling und Twirled Readout Error Extinction (TREX) zur Fehlerunterdrückung und -minderung. Die folgende Tabelle fasst zusammen, wie sich das großskalige Experiment vom kleinskaligen unterscheidet:

	Kleiner Maßstab	Großer Maßstab
Qubits	10	50
Zeitschritte	10	20
AQC-Checkpoints (1-schichtig + 2-schichtig)	3 + 2 = 5	6 + 4 = 10
Grundzustands-Ansatz-Schichten	3	5
Maximale MPS-Bindungsdimension	32	128
Estimator	`StatevectorEstimator`	QPU mit DD, Pauli-Twirling und TREX

XLA-Meldungen zu langsamer Kompilierung

Während der AQC-Optimierung (Schritt 2c unten) siehst du möglicherweise stderr-Meldungen wie:

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

Dabei handelt es sich um harmlose Diagnosemeldungen von XLA, dem Compiler, der die JAX-Autodiff-Funktionalität von qiskit-addon-aqc-tensor unterstützt. Bei 50 Qubits mit einer MPS-Bindungsdimension von 128 verbringt XLA beim ersten Tracing der Gradientenfunktion einige Minuten mit deren Kompilierung. Die Kompilierung ist erfolgreich, und die Optimierungsergebnisse sind davon nicht betroffen.

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

Ausgabe der vorherigen Code-Zelle

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

Ausgabe der vorherigen Code-Zelle

Die Hardware-Ergebnisse geben die zentralen Merkmale des Zwei-Spinon-Kontinuums wieder: Die Streuintensität konzentriert sich bei niedriger Energie in der Nähe des antiferromagnetischen Wellenvektors $q = \pi$ und ist nach unten durch die sinusförmige Spinon-Dispersion begrenzt, mit einem breiten Kontinuum spektralen Gewichts darüber statt eines einzelnen scharfen Modus. Dies ist dieselbe Struktur, die durch inelastische Neutronenstreuung an KCuF $_3$ gemessen wird, und sie bestätigt den vollständigen Workflow aus Grundzustandspräparation, Störung, AQC-komprimierter Trotter-Evolution und fehlergeminderter Messung im Maßstab von 50 Qubits.

Nächste Schritte

Wenn du diese Arbeit interessant fandest, könnte dich auch folgendes Material interessieren:

Empfehlungen

Neutronenstreuung mit einem AQC + Trotter-Dynamics-Serverless-Workflow simulieren — ein Begleit-Tutorial, das dasselbe Experiment über eine bereitgestellte Qiskit-Serverless-Funktionsvorlage ausführt
Lee et al., "Benchmarking quantum simulation with neutron-scattering experiments" (arXiv:2603.15608) — das Referenzpapier, auf dem dieses Tutorial basiert
Fehlerminderungs- und -unterdrückungstechniken — DD, Pauli-Twirling und TREX, die in den Hardware-Experimenten verwendet werden
Approximatives Quantenkompilieren für Zeitentwicklungs-Circuits — Tutorial zu AQC-Tensor
AQC-Tensor-Dokumentation

Lernziele​

Voraussetzungen​

Hintergrund​

Inelastische Neutronenstreuung und der dynamische Strukturfaktor​

KCuF3_33​: ein kanonischer Luttinger-Flüssigkeits-Magnet​

Was wir simulieren und messen​

Quantensimulations-Workflow​

Approximatives Quantenkompilieren (AQC)​

MPS-Simulation​

Anforderungen​

Setup​

Kleinskaliges Simulator-Beispiel​

Schritt 1: Klassische Eingaben auf ein Quantenproblem abbilden​

Schritt 2: Problem für die Ausführung auf Quantenhardware optimieren​

Schritt 3: Ausführung mit Qiskit-Primitiven​

Schritt 4: Nachbearbeitung und Rückgabe des Ergebnisses im gewünschten klassischen Format​

Großskalige Hardware-Ausführung​

Nächste Schritte​

Lernziele

Voraussetzungen

Hintergrund

Inelastische Neutronenstreuung und der dynamische Strukturfaktor

KCuF $_3$ : ein kanonischer Luttinger-Flüssigkeits-Magnet

Was wir simulieren und messen

Quantensimulations-Workflow

Approximatives Quantenkompilieren (AQC)

MPS-Simulation

Anforderungen

Setup

Kleinskaliges Simulator-Beispiel

Schritt 1: Klassische Eingaben auf ein Quantenproblem abbilden

Schritt 2: Problem für die Ausführung auf Quantenhardware optimieren

Schritt 3: Ausführung mit Qiskit-Primitiven

Schritt 4: Nachbearbeitung und Rückgabe des Ergebnisses im gewünschten klassischen Format

Großskalige Hardware-Ausführung

Nächste Schritte