{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "8e82ead7",
      "metadata": {},
      "source": [
        "---\n",
        "title: \"Simulation der Neutronenstreuung in Quantenmaterialien mit Quantenschaltungen\"\n",
        "description: \"Berechne den dynamischen Strukturfaktor eines Quantenmagneten mithilfe von Trotter-Schaltungen und einer MPS-Simulation.\"\n",
        "---\n",
        "\n",
        "{/* cspell:ignore DSFs spinon DMRG viridis fontsize vmax vmin cadetblue antiferromagnet COBYQA */}\n",
        "\n",
        "<span id=\"simulate-neutron-scattering-in-quantum-materials-with-quantum-circuits\" />\n",
        "\n",
        "# Simulation der Neutronenstreuung in Quantenmaterialien mit Quantenschaltungen\n",
        "\n",
        "*Geschätzte Laufzeit: 13 Minuten auf einem Heron- r2 -Prozessor (HINWEIS: Dies ist lediglich eine Schätzung. (Die Laufzeit kann variieren.)*\n",
        "\n",
        "<Admonition type=\"note\" title=\"Welches Tutorial soll ich verwenden?\">\n",
        "  In diesem Tutorial erfahren Sie Schritt für Schritt, wie die Implementierung funktioniert – sowohl bei der klassischen Berechnung, die lokal auf Ihrem Laptop ausgeführt wird, als auch bei der Hardware-Ausführung auf einer QPU von IBM Quantum®. Das groß angelegte Beispiel benötigt viel Speicherplatz, und die approximative Quantenkompilierung (AQC) kann auf einem typischen Laptop mehrere Stunden dauern. Um die Komprimierung auf die Rechen- und Speicherressourcen von „ Qiskit Serverless “ auszulagern, nutzen Sie das dazugehörige [Serverless-Tutorial](/docs/tutorials/simulate-neutron-scattering-with-a-serverless-workflow).\n",
        "</Admonition>\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "11033666",
      "metadata": {},
      "source": [
        "<span id=\"learning-outcomes\" />\n",
        "\n",
        "## Lernziele\n",
        "\n",
        "* Wie sich Spektren der unelastischen Neutronenstreuung (INS) auf dynamische Strukturfaktoren (DSFs) von Quanten-Spin-Modellen beziehen.\n",
        "* Wie man einen Grundzustand vorbereitet, eine lokale Störung anwendet und die Trotter-Zeitentwicklung auf einem Quantenschaltkreis durchführt.\n",
        "* Wie man Approximative Quantum Compiling (AQC) einsetzt, `qiskit-addon-aqc-tensor` um Deep-Trotter-Schaltkreise für die Ausführung auf Hardware zu komprimieren.\n",
        "* Wie man die retardierte Green-Funktion (RGF) aus den Erwartungswerten der Qubits extrahiert und sie mittels Fourier-Transformation in eine DSF umwandelt.\n",
        "\n",
        "<span id=\"prerequisites\" />\n",
        "\n",
        "## Voraussetzungen\n",
        "\n",
        "* [Grundlagen der Quanteninformation](/learning/courses/basics-of-quantum-information)\n",
        "* [Entwurf eines Variationsalgorithmus](/learning/courses/variational-algorithm-design)\n",
        "* [Einführung in Qiskit primitives (Schätzer und Stichprobenverfahren)](/docs/guides/qiskit-runtime-primitives)\n",
        "\n",
        "<span id=\"background\" />\n",
        "\n",
        "## Hintergrund\n",
        "\n",
        "In diesem Tutorial reproduzieren wir die Ergebnisse von [Lee et al., arXiv:2603.15608](https://arxiv.org/abs/2603.15608).\n",
        "\n",
        "Mithilfe der Neutronenstreuung können Forscher Materialien charakterisieren, die für die Nachhaltigkeit von Bedeutung sind, darunter auch Batterieelektroden. In diesem Tutorial wird untersucht, wie Quantenschaltungen magnetische Anregungen simulieren können und wie theoretische Modelle mit Streumessungen in Verbindung gebracht werden können. Anhand des Beispiels eines Quantenmagneten werden Methoden zur Untersuchung des Materialverhaltens entwickelt, die das Forschungsinstrumentarium für die Materialforschung erweitern.\n",
        "\n",
        "<span id=\"inelastic-neutron-scattering-and-the-dynamical-structure-factor\" />\n",
        "\n",
        "### Unelastische Neutronenstreuung und der dynamische Strukturfaktor\n",
        "\n",
        "Die unelastische Neutronenstreuung (INS) ist eine der leistungsfähigsten experimentellen Methoden zur Untersuchung magnetischer Anregungen in Quantenmaterialien. Wenn ein Strahl aus thermischen oder kalten Neutronen auf einen Kristall trifft, tauschen einzelne Neutronen sowohl Impuls $\\mathbf{q}$ als auch Energie $\\omega$ mit dem magnetischen Teilsystem aus. Die gemessene Streuintensität ist proportional zum dynamischen Strukturfaktor (DSF),\n",
        "\n",
        "$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,$\n",
        "\n",
        "das die vollständigen Raum-Zeit-Korrelationen der Spin-Freiheitsgrade kodiert.\n",
        "\n",
        "<span id=\"kcuf$_3$-a-canonical-luttinger-liquid-magnet\" />\n",
        "\n",
        "### KCuF$_3$: ein kanonischer Luttinger-Flüssigkeits-Magnet\n",
        "\n",
        "Kaliumkupferfluorid ( KCuF$_3$ ) ist ein quasi-eindimensionaler Antiferromagnet, in dem Ketten aus Spin-$\\frac{1}{2}$ Cu$^{2+}$ -Ionen über eine Heisenberg-Austausch $J$ zwischen nächsten Nachbarn wechselwirken, während die Kopplung zwischen den Ketten nur $\\sim 2.7\\%$ von $J$ beträgt. An $T = 6\\;\\mathrm{K}$, für die INS-Daten vorliegen, wird das Spektrum von fraktionierten Spinon-Anregungen dominiert, die charakteristisch für eine Tomonaga-Luttinger-Flüssigkeit sind. Da sich die Intrachain-Dynamik am isotropen Punkt gut durch den eindimensionalen Spin- $\\frac{1}{2}$ -XXZ-Hamiltonian beschreiben lässt ( $\\epsilon = 1$ ),\n",
        "\n",
        "$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],$\n",
        "\n",
        "KCuF$_3$ dient als idealer Testfall für die Quantensimulation: Der Hamilton-Operator ist einfach genug, um auf einem Quantenprozessor implementiert zu werden, doch der Grundzustand ist stark verschränkt, und das Anregungsspektrum weist ein breites Zwei-Spinon-Kontinuum auf.\n",
        "\n",
        "Hinweis: In dieser Anleitung wird $J = 1$ als Energieeinheit festgelegt und die Normierung $H = J\\sum_i[\\ldots]$ verwendet, was dem lokalen Hamilton-Operator $H_\\mathrm{loc} = J(S^X S^X + S^Y S^Y + S^Z S^Z)$ entspricht, der in der Schaltungsimplementierung der Arbeit verwendet wird (Abb. S3 (des Anhangs). Der vollständige Hamilton-Operator der Arbeit (Gleichung 3) weist einen zusätzlichen Gesamtfaktor von 2 auf, sodass der Wert $J$ des Papiers doppelt so hoch ist wie der hier verwendete Wert $J$.\n",
        "\n",
        "<span id=\"what-we-simulate-and-measure\" />\n",
        "\n",
        "### Was wir simulieren und messen\n",
        "\n",
        "Die von uns berechnete physikalische Größe ist die **verzögerte Green-Funktion** (RGF), definiert als die zeitabhängige Spin-Spin-Korrelationsfunktion\n",
        "\n",
        "$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,$\n",
        "\n",
        "wobei $j_c$ eine Referenzstelle (das Kettenzentrum) und $S_j^\\alpha(t) = e^{iHt}S_j^\\alpha e^{-iHt}$ der Spinoperator im Heisenberg-Bild ist. In diesem Tutorial konzentrieren wir uns auf die $zz$-Komponente ( $\\alpha = \\beta = z$ ). Auf einem Quantencomputer wird auf das RGF zugegriffen, indem der Grundzustand vorbereitet wird, eine lokale Störung an $j_c$ angelegt wird, der gestörte Zustand zeitlich weiterentwickelt wird und der Erwartungswert des einzelnen Qubits $\\langle\\sigma_j^z\\rangle$ an jedem Standort $j$ für jeden Zeitschritt gemessen wird. Die zentrale Erkenntnis ist, dass jede $\\langle\\sigma_j^z\\rangle$ die Differenz zur Magnetisierung im Grundzustand angibt; da der isotrope Heisenberg-Antiferromagnet eine Nettomagnetisierung pro Gitterpunkt von Null aufweist ( $\\langle\\sigma_j^z\\rangle_\\mathrm{GS} = 0$ ), ergibt sich aus dem rohen Messwert direkt $G^R(j, j_c, t)$, ohne dass eine explizite Subtraktion erforderlich ist.\n",
        "\n",
        "Durch die Erfassung von $G^R(j, j_c, t)$ über alle Orte und Zeitschritte hinweg stellen wir einen zweidimensionalen Datensatz zusammen, der anschließend sowohl räumlich als auch zeitlich einer Fourier-Transformation unterzogen wird, um den **dynamischen Strukturfaktor**$S(q,\\omega)$ zu ermitteln. Der DSF ist die Größe, die in einem INS-Experiment direkt gemessen wird: Er gibt Auskunft darüber, welche magnetischen Anregungen bei jedem Impuls $q$ und jeder Energie $\\omega$ vorliegen. Für die isotrope Heisenberg-Kette ist das exakte Anregungsspektrum ein **Zwei-Spinon-Kontinuum**, ein breites Band der Streuintensität, dessen Form als strenger End-to-End-Maßstab für die Quantensimulation dient: Es validiert gleichzeitig die Grundzustandspräparation, die Störung, die Trotter-Zeitentwicklung und das Messprotokoll.\n",
        "\n",
        "<span id=\"quantum-simulation-workflow\" />\n",
        "\n",
        "### Ablauf der Quantensimulation\n",
        "\n",
        "Der Arbeitsablauf spiegelt die physikalischen Abläufe eines INS-Ereignisses wider. Wir (1) bereiten die Vielteilchen-Grundzustands $|\\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})$ in der Mitte der Kette an, um den Spintransfer des Neutrons nachzubilden, (3) führen die Zeitentwicklung unter $H$ für diskrete Zeitschritte unter Verwendung der Trotterisierung zweiter Ordnung durch 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 die DSF- $S(q,\\omega)$.\n",
        "\n",
        "Die Größe, die wir bei jedem Zeitschritt messen, ist $\\sigma_i^z$ für jedes Qubit $i$. In Qiskit wird dies als Liste von `SparsePauliOp` Operatoren dargestellt: ein Einzel-Qubit- $Z$, eingebettet in die $n$ -Qubit-Identitätszeichenfolge für jeden Standort. Diese Observablen werden einmal pro Probleminstanz während des Problem-Mapping-Schritts (Schritt 1) gebildet und für jede Schaltung auf dieser Ebene wiederverwendet.\n",
        "\n",
        "<span id=\"approximate-quantum-compiling-aqc\" />\n",
        "\n",
        "### Approximatives Quanten-Compiling (AQC)\n",
        "\n",
        "Tiefe Trotter-Schaltkreise lassen sich mittels **approximativer Quantenkompilierung (AQC)** komprimieren, bei der die ersten paar Trotter-Schichten durch einen kürzeren parametrisierten Ansatz ersetzt werden, dessen Parameter klassisch so optimiert werden, dass die MPS-Fidelity gegenüber dem ursprünglichen tiefen Schaltkreis maximiert wird. Die verbleibenden Trotter-Schritte werden exakt angehängt, wodurch eine „gemischte“ AQC-+-Trotter-Schaltung mit deutlich weniger Zwei-Qubit-Gattern entsteht.\n",
        "\n",
        "<span id=\"mps-simulation\" />\n",
        "\n",
        "### MPS-Simulation\n",
        "\n",
        "Bei einem eindimensionalen System lassen sich mit Matrix-Product-State-Methoden (MPS) sowohl die Vorbereitung des Grundzustands (mit der Dichte-Matrix-Renormierungsgruppe, DMRG) als auch die Zeitentwicklung auf Schaltkreisebene effizient simulieren. Durch die Steuerung der Bindungsdimension $\\chi$ gehen wir einen Kompromiss zwischen Genauigkeit und Rechenaufwand ein. In diesem Tutorial verwenden wir die MPS-Simulation innerhalb von `qiskit-addon-aqc-tensor` zur Berechnung hochgenauer AQC-Ansätze, die die tiefen Trotter-Schaltungen für die Hardware-Ausführung komprimieren.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "de1ad5a2",
      "metadata": {},
      "source": [
        "<span id=\"requirements\" />\n",
        "\n",
        "## Anforderungen\n",
        "\n",
        "Bevor Sie mit dieser Anleitung beginnen, stellen Sie bitte sicher, dass folgende Komponenten installiert sind:\n",
        "\n",
        "* Qiskit SDK mit [Visual](/docs/api/qiskit/visualization) isierungsunterstützung\n",
        "* Qiskit Runtime (`pip install qiskit-ibm-runtime`)\n",
        "* `qiskit-addon-aqc-tensor` mit `quimb` und JAX-Extras (`pip install 'qiskit-addon-aqc-tensor[quimb-jax]'`)\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "9938e4bd",
      "metadata": {},
      "source": [
        "<span id=\"setup\" />\n",
        "\n",
        "## Einrichtung\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "id": "71835ff6",
      "metadata": {},
      "outputs": [],
      "source": [
        "import timeit\n",
        "import warnings\n",
        "from collections.abc import Iterator, Sequence\n",
        "from functools import partial\n",
        "\n",
        "import matplotlib.pyplot as plt\n",
        "import numpy as np\n",
        "import quimb.tensor as qtn\n",
        "import scipy.optimize\n",
        "from numpy.typing import NDArray\n",
        "from qiskit import QuantumCircuit\n",
        "from qiskit.circuit import CircuitInstruction, ParameterVector, Qubit\n",
        "from qiskit.circuit.library import PauliEvolutionGate\n",
        "from qiskit.primitives import StatevectorEstimator\n",
        "from qiskit.quantum_info import SparsePauliOp\n",
        "from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager\n",
        "from qiskit_addon_aqc_tensor import generate_ansatz_from_circuit\n",
        "from qiskit_addon_aqc_tensor.objective import MaximizeStateFidelity\n",
        "from qiskit_addon_aqc_tensor.simulation import tensornetwork_from_circuit\n",
        "from qiskit_addon_aqc_tensor.simulation.quimb import QuimbSimulator\n",
        "from qiskit_ibm_runtime import EstimatorV2 as Estimator\n",
        "from qiskit_ibm_runtime import QiskitRuntimeService\n",
        "from qiskit_quimb import quimb_circuit\n",
        "from scipy.sparse import SparseEfficiencyWarning\n",
        "\n",
        "# scipy.sparse.linalg.expm (reached via PauliEvolutionGate.to_matrix) prefers\n",
        "# CSC and warns when handed the CSR matrix from SparsePauliOp.to_matrix; the\n",
        "# result is unaffected, so silence the cosmetic warning.\n",
        "warnings.filterwarnings(\"ignore\", category=SparseEfficiencyWarning)\n",
        "\n",
        "\n",
        "def xxz_hamiltonian_mpo(\n",
        "    n_qubits: int, interaction: float = 1.0, anisotropy: float = 1.0\n",
        ") -> qtn.MatrixProductOperator:\n",
        "    \"\"\"1D XXZ Hamiltonian as a quimb MPO.\n",
        "\n",
        "    Builds the Hamiltonian using ``qtn.SpinHam1D``.\n",
        "\n",
        "    Args:\n",
        "        n_qubits: Number of sites.\n",
        "        interaction: Overall interaction strength (J in the paper).\n",
        "        anisotropy: XY/Z anisotropy (ε in the paper). ``ε = 1`` is the\n",
        "            isotropic Heisenberg point; ``ε = 0`` is the Ising limit.\n",
        "\n",
        "    Returns:\n",
        "        The Hamiltonian as a matrix product operator.\n",
        "    \"\"\"\n",
        "    builder = qtn.SpinHam1D(S=1 / 2)\n",
        "    builder += interaction * anisotropy * 0.5, \"+\", \"-\"\n",
        "    builder += interaction * anisotropy * 0.5, \"-\", \"+\"\n",
        "    builder += interaction, \"Z\", \"Z\"\n",
        "    return builder.build_mpo(L=n_qubits)\n",
        "\n",
        "\n",
        "def build_ground_state_ansatz(n_qubits: int, n_layers: int) -> QuantumCircuit:\n",
        "    \"\"\"Hamiltonian variational ansatz (HVA) circuit for ground-state preparation.\n",
        "\n",
        "    Starts from a product of singlet pairs and applies alternating\n",
        "    odd/even layers of parameterized XXZ pair-evolution gates.\n",
        "\n",
        "    The returned circuit is parameterized: it carries a ``ParameterVector``\n",
        "    named ``\"theta\"`` of length ``2 * n_layers`` whose values must be\n",
        "    assigned (e.g. via ``circuit.assign_parameters``) before simulation.\n",
        "    Element ``2 * r`` is the odd-layer angle and element ``2 * r + 1`` is the\n",
        "    even-layer angle of layer ``r``.\n",
        "\n",
        "    Args:\n",
        "        n_qubits: Number of qubits (must be even).\n",
        "        n_layers: Number of HVA layers.\n",
        "\n",
        "    Returns:\n",
        "        The parameterized HVA preparation circuit.\n",
        "    \"\"\"\n",
        "    theta = ParameterVector(\"theta\", 2 * n_layers)\n",
        "    circuit = QuantumCircuit(n_qubits)\n",
        "    # Initial singlet product state\n",
        "    for i in range(n_qubits // 2):\n",
        "        circuit.x(2 * i)\n",
        "        circuit.x(2 * i + 1)\n",
        "        circuit.h(2 * i + 1)\n",
        "        circuit.cx(2 * i + 1, 2 * i)\n",
        "    # Variational layers: each pair gate is exp(-i theta (XX+YY+ZZ)/2)\n",
        "    pair_ham = SparsePauliOp(\n",
        "        [\"XX\", \"YY\", \"ZZ\"], coeffs=[0.5, 0.5, 0.5]\n",
        "    )  # H_pair (HVA form)\n",
        "    for r in range(n_layers):\n",
        "        for i in range(1, (n_qubits + 1) // 2):  # odd layer\n",
        "            circuit.append(\n",
        "                PauliEvolutionGate(pair_ham, time=theta[2 * r]),\n",
        "                [2 * i - 1, 2 * i],\n",
        "            )\n",
        "        for i in range(n_qubits // 2):  # even layer\n",
        "            circuit.append(\n",
        "                PauliEvolutionGate(pair_ham, time=theta[2 * r + 1]),\n",
        "                [2 * i, 2 * i + 1],\n",
        "            )\n",
        "    return circuit\n",
        "\n",
        "\n",
        "def optimize_ground_state_ansatz(\n",
        "    ansatz: QuantumCircuit,\n",
        "    x0: NDArray[np.floating],\n",
        "    target_mps: qtn.MatrixProductState,\n",
        "    *,\n",
        "    max_bond: int | None = None,\n",
        "    cutoff: float = 1e-10,\n",
        "    method: str = \"COBYQA\",\n",
        "    options: dict | None = None,\n",
        ") -> scipy.optimize.OptimizeResult:\n",
        "    \"\"\"Optimize HVA parameters by maximizing fidelity with a target MPS.\n",
        "\n",
        "    The HVA circuit is simulated as a matrix product state with the given\n",
        "    bond-dimension truncation, and the parameters are optimized to maximize\n",
        "    the state fidelity ``|<psi_HVA | target_mps>|**2`` with the DMRG ground\n",
        "    state. Both states are normalized, so the minimized objective is the\n",
        "    infidelity ``1 - |<psi_HVA | target_mps>|**2``.\n",
        "\n",
        "    Args:\n",
        "        ansatz: Parameterized HVA circuit from ``build_ground_state_ansatz``.\n",
        "            The length of ``x0`` must equal ``ansatz.num_parameters``.\n",
        "        x0: Initial parameters.\n",
        "        target_mps: Target MPS (DMRG ground state) to maximize fidelity with.\n",
        "        max_bond: Maximum MPS bond dimension during gate application.\n",
        "        cutoff: Singular-value cutoff during gate application.\n",
        "        method: ``scipy.optimize.minimize`` method.\n",
        "        options: Options dict forwarded to ``scipy.optimize.minimize``.\n",
        "\n",
        "    Returns:\n",
        "        The Scipy OptimizeResult.\n",
        "    \"\"\"\n",
        "\n",
        "    def infidelity(params: NDArray[np.floating]) -> float:\n",
        "        circuit = ansatz.assign_parameters(params)\n",
        "        circuit_mps = quimb_circuit(\n",
        "            circuit.decompose([\"PauliEvolution\"]),\n",
        "            quimb_circuit_class=qtn.CircuitMPS,\n",
        "            max_bond=max_bond,\n",
        "            cutoff=cutoff,\n",
        "        )\n",
        "        return 1 - abs(circuit_mps.psi.H @ target_mps) ** 2\n",
        "\n",
        "    return scipy.optimize.minimize(\n",
        "        infidelity, np.asarray(x0), method=method, options=options\n",
        "    )\n",
        "\n",
        "\n",
        "def trotter_evolution(\n",
        "    qubits: Sequence[Qubit],\n",
        "    interaction: float,\n",
        "    anisotropy: float,\n",
        "    time_step: float,\n",
        "    n_steps: int,\n",
        ") -> Iterator[CircuitInstruction]:\n",
        "    \"\"\"Second-order Trotter steps of the XXZ pair Hamiltonian.\n",
        "\n",
        "    While the paper used a hand-optimized circuit for the Trotter steps, we use\n",
        "    PauliEvolutionGate here for simplicity and generality. The final two-qubit gate\n",
        "    count and gate depth are equivalent when transpiled with ``optimization_level=3``.\n",
        "\n",
        "    Args:\n",
        "        qubits: Qubits to act on (length ``n_qubits``).\n",
        "        interaction: Overall interaction strength (J in the paper).\n",
        "        anisotropy: XY/Z anisotropy (ε in the paper). ``ε = 1`` is the\n",
        "            isotropic Heisenberg point; ``ε = 0`` is the Ising limit.\n",
        "        time_step: Per-step Trotter time.\n",
        "        n_steps: Number of Trotter steps.\n",
        "\n",
        "    Yields:\n",
        "        ``CircuitInstruction``s implementing the Trotter steps.\n",
        "    \"\"\"\n",
        "    if n_steps == 0:\n",
        "        return\n",
        "    n_qubits = len(qubits)\n",
        "    pair_ham = SparsePauliOp(\n",
        "        [\"XX\", \"YY\", \"ZZ\"],\n",
        "        coeffs=[\n",
        "            0.25 * interaction * anisotropy,\n",
        "            0.25 * interaction * anisotropy,\n",
        "            0.25 * interaction,\n",
        "        ],\n",
        "    )\n",
        "    half_evo = PauliEvolutionGate(pair_ham, time=time_step / 2)\n",
        "    full_evo = PauliEvolutionGate(pair_ham, time=time_step)\n",
        "    for i in range(n_qubits // 2):  # half even layer\n",
        "        yield CircuitInstruction(half_evo, (qubits[2 * i], qubits[2 * i + 1]))\n",
        "    for i in range(n_qubits // 2 - 1):  # full odd layer\n",
        "        yield CircuitInstruction(\n",
        "            full_evo, (qubits[2 * i + 1], qubits[2 * i + 2])\n",
        "        )\n",
        "    for _ in range(n_steps - 1):  # interior steps\n",
        "        for i in range(n_qubits // 2):\n",
        "            yield CircuitInstruction(\n",
        "                full_evo, (qubits[2 * i], qubits[2 * i + 1])\n",
        "            )\n",
        "        for i in range(n_qubits // 2 - 1):\n",
        "            yield CircuitInstruction(\n",
        "                full_evo, (qubits[2 * i + 1], qubits[2 * i + 2])\n",
        "            )\n",
        "    for i in range(n_qubits // 2):  # half even layer\n",
        "        yield CircuitInstruction(half_evo, (qubits[2 * i], qubits[2 * i + 1]))\n",
        "\n",
        "\n",
        "def get_dsf(\n",
        "    n_qubits: int,\n",
        "    rgf_mat: NDArray[np.floating],\n",
        "    time_step: float,\n",
        "    n_steps: int,\n",
        "    n_points_momentum: int,\n",
        "    n_points_frequency: int,\n",
        ") -> NDArray[np.floating]:\n",
        "    \"\"\"Compute the dynamical structure factor from the retarded Green's function.\n",
        "\n",
        "    Uses the center-site approximation and a discrete Fourier transform.\n",
        "    The result is symmetrized about the momentum axis and clipped to\n",
        "    non-negative values, ready for plotting.\n",
        "\n",
        "    Args:\n",
        "        n_qubits: Number of qubits (sites).\n",
        "        rgf_mat: RGF matrix of shape ``(n_steps, n_qubits)``.\n",
        "        time_step: Trotter time-step size.\n",
        "        n_steps: Number of time steps.\n",
        "        n_points_momentum: Number of momentum points.\n",
        "        n_points_frequency: Number of frequency points.\n",
        "\n",
        "    Returns:\n",
        "        DSF array of shape ``(n_points_frequency, n_points_momentum)``,\n",
        "        symmetrized about the momentum axis and clipped to non-negative values.\n",
        "    \"\"\"\n",
        "    max_frequency = np.pi / time_step\n",
        "    momentum_range = np.linspace(0, 2 * np.pi, n_points_momentum)\n",
        "    frequency_range = np.linspace(0, max_frequency, n_points_frequency)\n",
        "    result = np.zeros((frequency_range.shape[0], momentum_range.shape[0]))\n",
        "    center = n_qubits // 2 - 1\n",
        "    for iw, w in enumerate(frequency_range):\n",
        "        exponent = np.exp(1j * w * time_step * np.arange(1, n_steps + 1))\n",
        "        # S = sigma/2, so two spin operators contribute a factor of (1/2)(1/2) = 1/4.\n",
        "        rgf_omega = (\n",
        "            np.dot(rgf_mat.T, exponent) * time_step / 4\n",
        "        )  # S(omega): time Fourier slice of the Green's function\n",
        "        for iq, q in enumerate(momentum_range):\n",
        "            momentum_phases = np.exp(\n",
        "                -1j * q * np.arange(-center, center + 2, 1)\n",
        "            )\n",
        "            result[iw, iq] = np.imag(np.dot(rgf_omega, momentum_phases))\n",
        "    result = -(result + result[:, ::-1]) / 2\n",
        "    result = np.clip(result, a_min=0, a_max=None)\n",
        "    return result\n",
        "\n",
        "\n",
        "def plot_dsf(\n",
        "    dsf: NDArray[np.floating],\n",
        "    time_step: float,\n",
        "    n_points_momentum: int,\n",
        "    n_points_frequency: int,\n",
        "    title: str | None = None,\n",
        ") -> None:\n",
        "    \"\"\"Heat-map of the dynamical structure factor.\n",
        "\n",
        "    Args:\n",
        "        dsf: DSF array of shape ``(n_points_frequency, n_points_momentum)``.\n",
        "        time_step: Trotter time-step size.\n",
        "        n_points_momentum: Number of momentum points.\n",
        "        n_points_frequency: Number of frequency points.\n",
        "        title: Optional plot title.\n",
        "    \"\"\"\n",
        "    max_frequency = np.pi / time_step\n",
        "    momentum_range = np.linspace(0, 2 * np.pi, n_points_momentum)\n",
        "    frequency_range = np.linspace(0, max_frequency, n_points_frequency)\n",
        "    x, y = np.meshgrid(momentum_range, frequency_range)\n",
        "    fig, ax = plt.subplots(figsize=(8, 5))\n",
        "    c = ax.pcolormesh(x, y, dsf / np.max(dsf), cmap=\"viridis\", shading=\"auto\")\n",
        "    fig.colorbar(c, ax=ax, label=\"Normalized intensity\")\n",
        "    ax.set_ylim(0, 3.6)\n",
        "    ax.set_xlim(0, 2 * np.pi)\n",
        "    ax.set_xlabel(r\"$q$\", fontsize=16)\n",
        "    ax.set_ylabel(r\"$\\tilde{\\omega} = \\omega / J$\", fontsize=16)\n",
        "    ax.set_xticks([0, np.pi / 2, np.pi, 3 * np.pi / 2, 2 * np.pi])\n",
        "    ax.set_xticklabels([\"0\", r\"$\\pi/2$\", r\"$\\pi$\", r\"$3\\pi/2$\", r\"$2\\pi$\"])\n",
        "    if title:\n",
        "        ax.set_title(title, fontsize=14)\n",
        "    plt.tight_layout()\n",
        "    plt.show()\n",
        "\n",
        "\n",
        "def plot_rgf(\n",
        "    n_qubits: int,\n",
        "    rgf_mat: NDArray[np.floating],\n",
        "    time_step: float,\n",
        "    n_steps: int,\n",
        "    title: str | None = None,\n",
        ") -> None:\n",
        "    \"\"\"Heat-map of the retarded Green's function in real space and time.\n",
        "\n",
        "    Args:\n",
        "        n_qubits: Number of qubits (sites).\n",
        "        rgf_mat: RGF matrix of shape ``(n_steps, n_qubits)``.\n",
        "        time_step: Trotter time-step size.\n",
        "        n_steps: Number of time steps.\n",
        "        title: Optional plot title.\n",
        "    \"\"\"\n",
        "    fig, ax = plt.subplots(figsize=(8, 6))\n",
        "    qubit_axis = np.arange(n_qubits)\n",
        "    t_axis = np.arange(1, n_steps + 1) * time_step\n",
        "    x, y = np.meshgrid(qubit_axis, t_axis)\n",
        "    c = ax.pcolormesh(\n",
        "        x,\n",
        "        y,\n",
        "        np.real(rgf_mat),\n",
        "        cmap=\"RdBu\",\n",
        "        vmax=0.5,\n",
        "        vmin=-0.5,\n",
        "        shading=\"auto\",\n",
        "    )\n",
        "    fig.colorbar(c, ax=ax, label=r\"Re $G^R(j, j_c, t)$\")\n",
        "    ax.set_xlabel(\"Qubit\", fontsize=16)\n",
        "    ax.xaxis.set_major_locator(\n",
        "        plt.matplotlib.ticker.MaxNLocator(integer=True)\n",
        "    )\n",
        "    ax.set_ylabel(r\"Time\", fontsize=16)\n",
        "    if title:\n",
        "        ax.set_title(title, fontsize=14)\n",
        "    plt.tight_layout()\n",
        "    plt.show()\n",
        "\n",
        "\n",
        "def uniform_2q_depth(circuit: QuantumCircuit) -> int:\n",
        "    \"\"\"Two-qubit gate depth in a standardized basis.\"\"\"\n",
        "    pass_manager = generate_preset_pass_manager(\n",
        "        optimization_level=0, basis_gates=[\"cz\", \"id\", \"rz\", \"sx\", \"x\"]\n",
        "    )\n",
        "    return pass_manager.run(circuit).depth(\n",
        "        lambda inst: inst.operation.num_qubits == 2\n",
        "    )"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "d28d01a6",
      "metadata": {},
      "source": [
        "<span id=\"small-scale-simulator-example\" />\n",
        "\n",
        "## Beispiel für einen Simulator im kleinen Maßstab\n",
        "\n",
        "Zunächst demonstrieren wir den gesamten Arbeitsablauf an **10 Qubits**, wobei wir den HVA-Grundzustandsansatz mithilfe einer MPS-Simulation optimieren und den Qiskit-Zustandsvektor-Simulator für die Zeitentwicklung verwenden. DMRG liefert eine Referenzenergie für den Grundzustand und ein Referenz-MPS. Anhand dieses Beispiels im kleinen Maßstab können wir jeden Schritt überprüfen, bevor wir das Projekt in größerem Maßstab umsetzen.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "3818ed8f",
      "metadata": {},
      "source": [
        "<span id=\"step-1-map-classical-inputs-to-a-quantum-problem\" />\n",
        "\n",
        "### Schritt 1: Klassische Eingaben auf ein Quantenproblem abbilden\n",
        "\n",
        "Zunächst definieren wir das physikalische Modell und konstruieren die Quantenschaltungen.\n",
        "\n",
        "**Hamiltonian.** KCuF$_3$ wird durch den Hamilton-Operator 1D XXZ am isotropen Punkt modelliert ( $J = 1$, $\\epsilon = 1$, wobei $J = 1$ als Energieeinheit festgelegt wird).\n",
        "\n",
        "**Grundzustand.** Wir verwenden eine Schaltung nach dem `build_ground_state_ansatz`Hamilton-Variationsansatz (HVA) als Schaltung zur Vorbereitung des Grundzustands. `CircuitMPS`Die HVA-Parameter werden klassisch durch Maximierung der Zustandsgenauigkeit $|\\langle\\psi_{\\mathrm{HVA}}(\\theta)|\\psi_{\\mathrm{DMRG}}\\rangle|^2$ mit dem DMRG-Grundzustands-MPS optimiert, wobei der HVA-Zustand durch Simulation der Schaltung als Matrixproduktzustand mit quimb's ermittelt wird. Für die Optimierung verwenden `scipy.optimize.minimize` wir. Zur Vergleichszwecken wird auch die Energie $\\langle H \\rangle$ e des optimierten Ansatzes berechnet.\n",
        "\n",
        "**Trotter-Tore.** `PauliEvolutionGate(H_pair, time=time_step)`Jeder „Nearest-Neighbor“-Wechselwirkungsterm $e^{-i\\Delta t\\, H_{\\mathrm{pair}}}$ mit $H_{\\mathrm{pair}} = (J/4)(XX + YY + ZZ)$ wird wie folgt gebildet:. Qiskit synthetisiert dies während der Transpilation zu der optimalen Drei-CNOT-Zerlegung.\n",
        "\n",
        "**Störung.** Ein $R_z(\\pi/2)$-Gatter am zentralen Qubit implementiert $U_{j_c} = \\frac{1}{\\sqrt{2}}(I - i\\sigma^z_{j_c})$ und ahmt so den lokalen Spin-Flip nach, der durch ein gestreutes Neutron erzeugt wird.\n",
        "\n",
        "**Beobachtbare Größen.** Wir konstruieren für jede Qubit-Position eine Beobachtungsgröße vom Typ $\\sigma^z$. Diese `SparsePauliOp` Objekte werden in Schritt 3 an die Estimator-Primitive übergeben, um bei jedem Zeitschritt die Werte von $\\langle\\sigma_i^z\\rangle$ zu ermitteln.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "id": "58305fbc",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Ground-state energy (DMRG): -4.258035\n",
            "Optimizing ground state ansatz...\n",
            "Finished optimizing ground state ansatz in 8.850687621976249 seconds.\n",
            "Ground state ansatz fidelity: 0.984277\n",
            "Ground state ansatz energy: -4.232565\n",
            "Built 10 circuits, deepest 2q depth (uniform basis) = 163\n",
            "Defined 10 Z observables.\n"
          ]
        }
      ],
      "source": [
        "# -- Physical parameters --\n",
        "n_qubits = 10\n",
        "interaction = 1.0  # J\n",
        "anisotropy = 1.0  # ε (isotropic point)\n",
        "time_step = 0.6\n",
        "n_steps = 10\n",
        "mps_max_bond = 32\n",
        "mps_cutoff = 1e-8\n",
        "center = n_qubits // 2 - 1\n",
        "\n",
        "# -- Hamiltonian MPO --\n",
        "ham_mpo = xxz_hamiltonian_mpo(n_qubits, interaction, anisotropy)\n",
        "\n",
        "# -- Reference ground-state energy via DMRG --\n",
        "dmrg = qtn.DMRG2(ham_mpo)\n",
        "dmrg.solve(tol=1e-8)\n",
        "print(f\"Ground-state energy (DMRG): {dmrg.energy:.6f}\")\n",
        "\n",
        "# -- Build ground state ansatz circuit --\n",
        "gs_n_layers = 3\n",
        "gs_ansatz = build_ground_state_ansatz(n_qubits, gs_n_layers)\n",
        "\n",
        "# -- Optimize ground state ansatz parameters to maximize fidelity with the DMRG MPS --\n",
        "# Initialize odd-layer angles near 0 (where the inter-pair gate is the\n",
        "# identity) and even-layer angles near pi/2 (where the intra-pair gate\n",
        "# is a SWAP, since 0.5*(XX + YY + ZZ) = SWAP - I/2).\n",
        "rng = np.random.default_rng(12345)\n",
        "x0 = np.tile([0, np.pi / 2], gs_n_layers) + rng.normal(\n",
        "    scale=0.1, size=2 * gs_n_layers\n",
        ")\n",
        "print(\"Optimizing ground state ansatz...\")\n",
        "t0 = timeit.default_timer()\n",
        "result = optimize_ground_state_ansatz(\n",
        "    gs_ansatz,\n",
        "    x0,\n",
        "    dmrg.state,\n",
        "    max_bond=mps_max_bond,\n",
        "    cutoff=mps_cutoff,\n",
        "    options=dict(maxiter=100),\n",
        ")\n",
        "t1 = timeit.default_timer()\n",
        "print(f\"Finished optimizing ground state ansatz in {t1 - t0} seconds.\")\n",
        "print(f\"Ground state ansatz fidelity: {1 - result.fun:.6f}\")\n",
        "\n",
        "gs_circuit = gs_ansatz.assign_parameters(result.x)\n",
        "gs_circuit_mps = quimb_circuit(\n",
        "    gs_circuit.decompose([\"PauliEvolution\"]),\n",
        "    quimb_circuit_class=qtn.CircuitMPS,\n",
        "    max_bond=mps_max_bond,\n",
        "    cutoff=mps_cutoff,\n",
        ")\n",
        "gs_ansatz_energy = qtn.expec_TN_1D(\n",
        "    gs_circuit_mps.psi.H, ham_mpo, gs_circuit_mps.psi\n",
        ")\n",
        "print(f\"Ground state ansatz energy: {gs_ansatz_energy:.6f}\")\n",
        "\n",
        "# -- Build circuits for each time step --\n",
        "perturbed = gs_circuit.copy()\n",
        "perturbed.rz(np.pi / 2, center)\n",
        "\n",
        "circuits = []\n",
        "for t in range(1, n_steps + 1):\n",
        "    circuit = perturbed.copy()\n",
        "    for instr in trotter_evolution(\n",
        "        circuit.qubits, interaction, anisotropy, time_step, t\n",
        "    ):\n",
        "        circuit.append(instr)\n",
        "    circuits.append(circuit)\n",
        "print(\n",
        "    f\"Built {len(circuits)} circuits, deepest 2q depth (uniform basis) = \"\n",
        "    f\"{uniform_2q_depth(circuits[-1])}\"\n",
        ")\n",
        "\n",
        "# -- Observables: Z on each qubit site --\n",
        "observables = [\n",
        "    SparsePauliOp.from_sparse_list([(\"Z\", [i], 1)], num_qubits=n_qubits)\n",
        "    for i in range(n_qubits)\n",
        "]\n",
        "print(f\"Defined {len(observables)} Z observables.\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "39191552",
      "metadata": {},
      "source": [
        "<span id=\"step-2-optimize-problem-for-quantum-hardware-execution\" />\n",
        "\n",
        "### Schritt 2: Optimierung des Problems für die Ausführung auf Quantenhardware\n",
        "\n",
        "Für echte Hardware wären die oben genannten Trotter-Schaltungen zu komplex. **Das „Approximate Quantum Compiling“ (AQC)** löst dieses Problem, indem es die ersten $k$ Trotter-Schichten (einschließlich des Grundzustandsschaltkreises) durch einen kürzeren, parametrisierten Ansatz ersetzt, der so optimiert ist, dass die MPS-Fidelity gegenüber dem ursprünglichen tiefen Schaltkreis maximiert wird. Die verbleibenden Trotter-Schritte werden exakt angehängt, wodurch eine flachere „AQC + Trotter“-Schaltung entsteht.\n",
        "\n",
        "Um einen Ausgleich zwischen Ausdruckskraft und Schaltungstiefe zu schaffen, werden zwei Ansätze verwendet: Ein **einlagiger Ansatz** (der aus einem einzigen Trotter-Schritt generiert wird) komprimiert die frühesten Zeitschritte, während ein tieferer **zweilagiger Ansatz** (der aus zwei Trotter-Schritten generiert wird) die nächsten Schritte komprimiert, bei denen eine höhere Genauigkeit erforderlich ist.\n",
        "\n",
        "Der AQC-Arbeitsablauf umfasst vier Teilschritte:\n",
        "\n",
        "1. **Zielschaltungen aufbauen** – die ersten $k_1 + k_2$ -Schaltungen aus Schritt 1 dienen direkt als AQC-Ziele.\n",
        "2. **MPS des Ziels berechnen** – jede Zielschaltung als Matrixproduktzustand simulieren, indem man `quimb.tensor.CircuitMPS`… verwendet.\n",
        "3. **Ansätze generieren und optimieren** — `generate_ansatz_from_circuit` erstellt einen einschichtigen und einen zweischichtigen parametrisierten Ansatz; die Parameter werden mithilfe von L-BFGS-B mit JAX-beschleunigten Gradienten optimiert, um $1 - |\\langle\\psi_{\\mathrm{ansatz}}|\\psi_{\\mathrm{target}}\\rangle|^2$ zu minimieren. Die ersten Schritte $k_1$ verwenden den einlagigen Ansatz, die nächsten Schritte $k_2$ den zweilagigen Ansatz; innerhalb jeder Stufe erfolgt bei jedem Schritt ein Warmstart ausgehend von den optimierten Parametern des vorherigen Schritts, und an der Stufengrenze werden die Parameter auf die Standardwerte der jeweiligen Stufe zurückgesetzt.\n",
        "4. `trotter_evolution`**Mischschaltungen zusammenstellen** – für Zeitschritte jenseits des AQC-Checkpoints fügen Sie dem optimierten zweischichtigen AQC-Schaltkreis mithilfe von exakte Trotter-Schichten an.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "id": "e78f8d9a",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "\n",
            "Step 2b — target MPS:\n",
            "  k=1: max bond = 22\n",
            "  k=2: max bond = 22\n",
            "  k=3: max bond = 26\n",
            "  k=4: max bond = 27\n",
            "  k=5: max bond = 30\n",
            "\n",
            "Step 2c — one-layer ansatz: 389 parameters, 2q depth (uniform basis) = 27\n",
            "Step 2c — two-layer ansatz: 470 parameters, 2q depth (uniform basis) = 33\n",
            "  k=1: fidelity = 1.0000, 2q depth (uniform basis) = 27, 11.8s\n",
            "  k=2: fidelity = 0.9980, 2q depth (uniform basis) = 27, 13.2s\n",
            "  k=3: fidelity = 0.9890, 2q depth (uniform basis) = 27, 13.7s\n",
            "  k=4: fidelity = 0.9983, 2q depth (uniform basis) = 33, 16.4s\n",
            "  k=5: fidelity = 0.9957, 2q depth (uniform basis) = 33, 16.4s\n",
            "\n",
            "Step 2d — assembled 10 circuits\n",
            "  At step 10 (uniform basis): Trotter 2q depth = 163, AQC+Trotter 2q depth = 99; 2q gates: Trotter = 371, AQC+Trotter = 200\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/simulate-neutron-scattering/extracted-outputs/e78f8d9a-1.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "# Number of time steps to compress into an AQC ansatz with one layer\n",
        "aqc_n_steps_1 = 3\n",
        "# Number of time steps to compress into an AQC ansatz with two layers\n",
        "aqc_n_steps_2 = 2\n",
        "aqc_n_steps_total = aqc_n_steps_1 + aqc_n_steps_2\n",
        "\n",
        "# ── Step 2a: Target circuits (first aqc_n_steps_total circuits from Step 1) ──\n",
        "target_circuits = {\n",
        "    k: circuits[k - 1] for k in range(1, aqc_n_steps_total + 1)\n",
        "}\n",
        "# ── Step 2b: Compute target MPS ──\n",
        "# Decompose PauliEvolutionGate to RXX/RYY/RZZ before passing to the AQC MPS\n",
        "# backend, which does not understand PauliEvolutionGate natively.\n",
        "aqc_sim = QuimbSimulator(\n",
        "    quimb_circuit_factory=partial(\n",
        "        qtn.CircuitMPS,\n",
        "        gate_opts=dict(max_bond=mps_max_bond, cutoff=mps_cutoff),\n",
        "    ),\n",
        "    autodiff_backend=\"jax\",\n",
        ")\n",
        "print(\"\\nStep 2b — target MPS:\")\n",
        "target_mps = {}\n",
        "for k in range(1, aqc_n_steps_total + 1):\n",
        "    target_mps[k] = tensornetwork_from_circuit(\n",
        "        target_circuits[k].decompose([\"PauliEvolution\"]), aqc_sim\n",
        "    )\n",
        "    print(f\"  k={k}: max bond = {target_mps[k].psi.max_bond()}\")\n",
        "\n",
        "# ── Step 2c: Generate ansätze (one-layer and two-layer) and optimize parameters ──\n",
        "ansatz_1, initial_params_1 = generate_ansatz_from_circuit(\n",
        "    target_circuits[1].decompose([\"PauliEvolution\"]),\n",
        "    qubits_initially_zero=True,\n",
        ")\n",
        "initial_params_1 = np.array(initial_params_1)\n",
        "ansatz_2, initial_params_2 = generate_ansatz_from_circuit(\n",
        "    target_circuits[2].decompose([\"PauliEvolution\"]),\n",
        "    qubits_initially_zero=True,\n",
        ")\n",
        "initial_params_2 = np.array(initial_params_2)\n",
        "print(\n",
        "    f\"\\nStep 2c — one-layer ansatz: {ansatz_1.num_parameters} parameters, \"\n",
        "    f\"2q depth (uniform basis) = {uniform_2q_depth(ansatz_1)}\"\n",
        ")\n",
        "print(\n",
        "    f\"Step 2c — two-layer ansatz: {ansatz_2.num_parameters} parameters, \"\n",
        "    f\"2q depth (uniform basis) = {uniform_2q_depth(ansatz_2)}\"\n",
        ")\n",
        "\n",
        "aqc_circuits = {}\n",
        "aqc_params = {}\n",
        "for k in range(1, aqc_n_steps_total + 1):\n",
        "    if k <= aqc_n_steps_1:\n",
        "        ansatz, base_params = ansatz_1, initial_params_1\n",
        "    else:\n",
        "        ansatz, base_params = ansatz_2, initial_params_2\n",
        "    # Warm-start from the previous step only within the same stage\n",
        "    same_stage = (k - 1 >= 1) and (\n",
        "        (k - 1 <= aqc_n_steps_1) == (k <= aqc_n_steps_1)\n",
        "    )\n",
        "    x0 = aqc_params[k - 1] if same_stage else base_params\n",
        "    obj = MaximizeStateFidelity(target_mps[k], ansatz, aqc_sim)\n",
        "    t0 = timeit.default_timer()\n",
        "    result = scipy.optimize.minimize(\n",
        "        obj.loss_function,\n",
        "        x0,\n",
        "        method=\"L-BFGS-B\",\n",
        "        jac=True,\n",
        "        options=dict(maxiter=100),\n",
        "    )\n",
        "    elapsed = timeit.default_timer() - t0\n",
        "    aqc_params[k] = result.x\n",
        "    aqc_circuits[k] = ansatz.assign_parameters(result.x)\n",
        "    print(\n",
        "        f\"  k={k}: fidelity = {1 - result.fun:.4f}, \"\n",
        "        f\"2q depth (uniform basis) = \"\n",
        "        f\"{uniform_2q_depth(aqc_circuits[k])}, \"\n",
        "        f\"{elapsed:.1f}s\"\n",
        "    )\n",
        "\n",
        "# ── Step 2d: Assemble full circuit set (AQC + Trotter) ──\n",
        "all_circuits = []\n",
        "for k in range(1, aqc_n_steps_total + 1):\n",
        "    all_circuits.append(aqc_circuits[k])\n",
        "base = aqc_circuits[aqc_n_steps_total]\n",
        "for k in range(1, n_steps - aqc_n_steps_total + 1):\n",
        "    circuit = base.copy()\n",
        "    for instr in trotter_evolution(\n",
        "        circuit.qubits, interaction, anisotropy, time_step, k\n",
        "    ):\n",
        "        circuit.append(instr)\n",
        "    all_circuits.append(circuit)\n",
        "\n",
        "full_depths = [\n",
        "    uniform_2q_depth(circuits[k - 1]) for k in range(1, n_steps + 1)\n",
        "]\n",
        "aqc_depths = [uniform_2q_depth(circuit) for circuit in all_circuits]\n",
        "full_2q = [\n",
        "    circuits[k - 1].decompose([\"PauliEvolution\"]).num_nonlocal_gates()\n",
        "    for k in range(1, n_steps + 1)\n",
        "]\n",
        "aqc_2q = [\n",
        "    circuit.decompose([\"PauliEvolution\"]).num_nonlocal_gates()\n",
        "    for circuit in all_circuits\n",
        "]\n",
        "print(f\"\\nStep 2d — assembled {len(all_circuits)} circuits\")\n",
        "print(\n",
        "    f\"  At step {n_steps} (uniform basis): \"\n",
        "    f\"Trotter 2q depth = {full_depths[-1]}, AQC+Trotter 2q depth = {aqc_depths[-1]}; \"\n",
        "    f\"2q gates: Trotter = {full_2q[-1]}, AQC+Trotter = {aqc_2q[-1]}\"\n",
        ")\n",
        "\n",
        "steps_axis = np.arange(1, n_steps + 1)\n",
        "fig, ax = plt.subplots(figsize=(7, 4))\n",
        "ax.plot(steps_axis, full_depths, \"-o\", color=\"black\", label=\"Full Trotter\")\n",
        "ax.plot(\n",
        "    steps_axis, aqc_depths, \"-o\", color=\"cadetblue\", label=\"AQC + Trotter\"\n",
        ")\n",
        "ax.set_xlabel(\"Trotter steps\", fontsize=13)\n",
        "ax.xaxis.set_major_locator(plt.matplotlib.ticker.MaxNLocator(integer=True))\n",
        "ax.set_ylabel(\"2q gate depth (uniform basis)\", fontsize=13)\n",
        "ax.set_title(f\"AQC circuit-depth reduction ({n_qubits} qubits)\", fontsize=13)\n",
        "ax.legend(fontsize=11)\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "43ad6048",
      "metadata": {},
      "source": [
        "<span id=\"step-3-execute-using-qiskit-primitives\" />\n",
        "\n",
        "### Schritt 3: Mit Qiskit primitives ausführen\n",
        "\n",
        "Wir simulieren jede mit AQC kompilierte Schaltung mithilfe der `StatevectorEstimator` Primitive.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "id": "b8f9e41c",
      "metadata": {},
      "outputs": [],
      "source": [
        "estimator = StatevectorEstimator()\n",
        "pubs = [(circuit, observables) for circuit in all_circuits]\n",
        "job = estimator.run(pubs)\n",
        "result = job.result()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "cb12662c",
      "metadata": {},
      "source": [
        "<span id=\"step-4-post-process-and-return-result-in-desired-classical-format\" />\n",
        "\n",
        "### Schritt 4: Nachbearbeitung und Ausgabe des Ergebnisses im gewünschten klassischen Format\n",
        "\n",
        "Nun berechnen wir den Erwartungswert $\\langle\\sigma_i^z\\rangle$ für jedes Qubit $i$ bei jedem Zeitschritt $t$. Diese Werte bilden die Matrix der verzögerten Green-Funktion $G^R(j, j_c, t)$. Anschließend führen wir eine Fourier-Transformation der RGF in den dynamischen Strukturfaktor $S(q, \\omega)$ durch und stellen sowohl die RGF als auch den DSF grafisch dar. `get_dsf` wendet Spiegelsymmetrie an und beschneidet negative Werte intern.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "id": "4585d139",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/simulate-neutron-scattering/extracted-outputs/4585d139-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        },
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/simulate-neutron-scattering/extracted-outputs/4585d139-1.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "rgf_mat = np.stack([pub_result.data.evs for pub_result in result])\n",
        "\n",
        "# -- Compute DSF --\n",
        "n_points_momentum, n_points_frequency = 100, 100\n",
        "spectrum = get_dsf(\n",
        "    n_qubits,\n",
        "    rgf_mat,\n",
        "    time_step,\n",
        "    n_steps,\n",
        "    n_points_momentum,\n",
        "    n_points_frequency,\n",
        ")\n",
        "\n",
        "# -- Plot retarded Green's function --\n",
        "plot_rgf(\n",
        "    n_qubits,\n",
        "    rgf_mat,\n",
        "    time_step,\n",
        "    n_steps,\n",
        "    title=f\"Retarded Green's function — {n_qubits} qubits (simulation)\",\n",
        ")\n",
        "\n",
        "# -- Plot DSF --\n",
        "plot_dsf(\n",
        "    spectrum,\n",
        "    time_step,\n",
        "    n_points_momentum,\n",
        "    n_points_frequency,\n",
        "    title=f\"Dynamical structure factor — {n_qubits} qubits (simulation)\",\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "ls-intro",
      "metadata": {},
      "source": [
        "<span id=\"large-scale-hardware-execution\" />\n",
        "\n",
        "## Hardware-Ausführung in großem Maßstab\n",
        "\n",
        "Wir skalieren nun auf **50 Qubits**. In diesem Maßstab erreichen die Optimierungen geringere Fidelity-Werte als im Beispiel im kleinen Maßstab: Die Fidelity des Grundzustands-Ansatzes sinkt auf etwa 0.65, und die AQC-Fidelity-Werte sinken an den späteren Kontrollpunkten auf etwa 0.7. Dies ist zu erwarten, und Sie können diese Genauigkeiten verbessern, indem Sie die Anzahl der Ansatzschichten für den Grundzustand (`gs_n_layers`) oder die Anzahl der Optimierungsiterationen (`maxiter`) erhöhen – was jedoch mit zusätzlichen klassischen Rechenkosten verbunden ist. Beachten Sie außerdem, dass die AQC-Fidelity-Werte bei zweischichtigen Systemen niedriger ausfallen als bei einschichtigen. Dies ist kein Rückschritt: Spätere Zeitschritte erzeugen mehr Verschränkung und lassen sich einfach schwerer komprimieren, weshalb für sie der ausdrucksstärkere zweischichtige Ansatz verwendet wird.\n",
        "\n",
        "Die AQC-Optimierung kann zudem mehrere Stunden klassischer Rechenzeit in Anspruch nehmen (in dem hier gezeigten Durchlauf etwa sechs Stunden, wobei der größte Teil dieser Zeit auf die zweischichtigen Checkpoints in Schritt 2c entfällt). Um die Gesamtlaufzeit zu verkürzen, sollten Sie in Erwägung ziehen, dieses Notebook auf leistungsstärkerer herkömmlicher Hardware auszuführen, beispielsweise auf einem Hochleistungsrechnersystem (HPC). Alternativ können Sie das Problem auf eine kleinere Instanz reduzieren (beispielsweise mit weniger Qubits oder weniger Zeitschritten); Ihre Ergebnisse werden sich dann von den hier gezeigten unterscheiden.\n",
        "\n",
        "Der folgende Code folgt derselben vierstufigen Struktur wie das Beispiel im kleinen Maßstab. Die Grundzustandsparameter des HVA werden erneut mithilfe einer MPS-Simulation optimiert. Auf der QPU setzen wir dynamische Entkopplung (DD), Pauli-Twirling und Twirled Readout Error Extinction (TREX) zur Fehlerunterdrückung und -minderung ein. Die folgende Tabelle fasst zusammen, inwiefern sich das Großversuch vom Kleinversuch unterscheidet:\n",
        "\n",
        "|                                        | In kleinem Maßstab     | Großformat                        |\n",
        "| -------------------------------------- | ---------------------- | --------------------------------- |\n",
        "| Qubits                                 | 10                     | 50                                |\n",
        "| Zeitschritte                           | 10                     | 20 Jahre                          |\n",
        "| AQC-Kontrollpunkte (1-lagig + 2-lagig) | 3 + 2 = 5              | 6 + 4 = 10                        |\n",
        "| Ansatzschichten für den Grundzustand   | 3                      | 5                                 |\n",
        "| Maximale Abmessung der MPS-Verbindung  | 32                     | 128                               |\n",
        "| Schätzer                               | `StatevectorEstimator` | QPU mit DD, Pauli-Wirbel und TREX |\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "9e85fbd1",
      "metadata": {},
      "source": [
        "<Admonition type=\"note\" title=\"XLA-Meldungen zur langsamen Kompilierung\">\n",
        "  Während der AQC-Optimierung (Schritt 2c weiter unten) werden möglicherweise Meldungen wie die folgenden angezeigt `stderr` :\n",
        "\n",
        "  ```\n",
        "  [Compiling module jit_MakeArrayFn for CPU] Very slow compile? ...\n",
        "  The operation took 2m14s\n",
        "  ```\n",
        "\n",
        "  `qiskit-addon-aqc-tensor`Dies sind harmlose Diagnosen von XLA, dem Compiler, der hinter dem JAX-Autodiff von steht. Bei 50 Qubits und einer MPS-Bindungsdimension von 128 benötigt XLA beim ersten Durchlaufen der Gradientenfunktion einige Minuten für die Kompilierung. Die Kompilierung verläuft erfolgreich, und die Optimierungsergebnisse bleiben davon unberührt.\n",
        "</Admonition>\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "id": "fbe57e1a",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Ground-state energy (DMRG): -21.972109\n",
            "Optimizing ground state ansatz...\n",
            "Finished optimizing ground state ansatz in 132.2076231740648 seconds.\n",
            "Ground state ansatz fidelity: 0.645956\n",
            "Ground state ansatz energy: -21.616744\n",
            "Built 20 circuits, deepest 2q depth (uniform basis) = 307\n",
            "Defined 50 Z observables.\n",
            "\n",
            "Step 2b — target MPS:\n",
            "  k=1: max bond = 44\n",
            "  k=2: max bond = 46\n",
            "  k=3: max bond = 53\n",
            "  k=4: max bond = 62\n",
            "  k=5: max bond = 75\n",
            "  k=6: max bond = 96\n",
            "  k=7: max bond = 118\n",
            "  k=8: max bond = 128\n",
            "  k=9: max bond = 128\n",
            "  k=10: max bond = 128\n",
            "\n",
            "Step 2c — one-layer ansatz: 2971 parameters, 2q depth (uniform basis) = 39\n",
            "Step 2c — two-layer ansatz: 3412 parameters, 2q depth (uniform basis) = 45\n",
            "  k=1: fidelity = 1.0000, 2q depth (uniform basis) = 39, 215.6s\n",
            "  k=2: fidelity = 0.9676, 2q depth (uniform basis) = 39, 269.8s\n",
            "  k=3: fidelity = 0.9014, 2q depth (uniform basis) = 39, 293.5s\n",
            "  k=4: fidelity = 0.8755, 2q depth (uniform basis) = 39, 279.6s\n",
            "  k=5: fidelity = 0.8450, 2q depth (uniform basis) = 39, 266.4s\n",
            "  k=6: fidelity = 0.8146, 2q depth (uniform basis) = 39, 371.0s\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "E0626 02:59:07.229070 1508496 slow_operation_alarm.cc:73] \n",
            "********************************\n",
            "[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.\n",
            "********************************\n",
            "E0626 02:59:21.713796 1508477 slow_operation_alarm.cc:140] The operation took 2m14.484932074s\n",
            "\n",
            "********************************\n",
            "[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.\n",
            "********************************\n"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "  k=7: fidelity = 0.7312, 2q depth (uniform basis) = 45, 10865.3s\n",
            "  k=8: fidelity = 0.7720, 2q depth (uniform basis) = 45, 1450.3s\n",
            "  k=9: fidelity = 0.7316, 2q depth (uniform basis) = 45, 3388.6s\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "E0626 07:21:05.148724 1508496 slow_operation_alarm.cc:73] \n",
            "********************************\n",
            "[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.\n",
            "********************************\n",
            "E0626 07:21:24.286095 1508477 slow_operation_alarm.cc:140] The operation took 2m19.137566498s\n",
            "\n",
            "********************************\n",
            "[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.\n",
            "********************************\n"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "  k=10: fidelity = 0.6726, 2q depth (uniform basis) = 45, 3396.1s\n",
            "\n",
            "Step 2d — assembled 20 circuits\n",
            "  At step 20 (uniform basis): Trotter 2q depth = 307, AQC+Trotter 2q depth = 171; 2q gates: Trotter = 3775, AQC+Trotter = 1913\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/simulate-neutron-scattering/extracted-outputs/fbe57e1a-5.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Backend: ibm_fez\n",
            "Transpiled 2q depth (deepest, ISA on ibm_fez): 105 (2q gates: 2574)\n",
            "Job ID: d8v39vhropqc738biotg\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/simulate-neutron-scattering/extracted-outputs/fbe57e1a-7.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        },
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/simulate-neutron-scattering/extracted-outputs/fbe57e1a-8.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "# ── Parameters ──────────────────────────────────────────────────────────────\n",
        "n_qubits = 50  # 10 → 50\n",
        "interaction = 1.0  # J\n",
        "anisotropy = 1.0  # ε (isotropic point)\n",
        "time_step = 0.6\n",
        "n_steps = 20  # 10 → 20\n",
        "aqc_n_steps_1 = 6  # 3 → 6\n",
        "aqc_n_steps_2 = 4  # 2 → 4\n",
        "aqc_n_steps_total = aqc_n_steps_1 + aqc_n_steps_2\n",
        "gs_n_layers = 5  # 3 → 5\n",
        "mps_max_bond = 128  # 32 -> 128\n",
        "mps_cutoff = 1e-8\n",
        "center = n_qubits // 2 - 1\n",
        "\n",
        "# ── Step 1: Map ──────────────────────────────────────────────────────────────\n",
        "ham_mpo = xxz_hamiltonian_mpo(n_qubits, interaction, anisotropy)\n",
        "\n",
        "dmrg = qtn.DMRG2(ham_mpo)\n",
        "dmrg.solve(tol=1e-8)\n",
        "print(f\"Ground-state energy (DMRG): {dmrg.energy:.6f}\")\n",
        "\n",
        "gs_ansatz = build_ground_state_ansatz(n_qubits, gs_n_layers)\n",
        "\n",
        "rng = np.random.default_rng(12345)\n",
        "x0 = np.tile([0.0, np.pi / 2], gs_n_layers) + rng.normal(\n",
        "    scale=0.1, size=2 * gs_n_layers\n",
        ")\n",
        "print(\"Optimizing ground state ansatz...\")\n",
        "t0 = timeit.default_timer()\n",
        "result = optimize_ground_state_ansatz(\n",
        "    gs_ansatz,\n",
        "    x0,\n",
        "    dmrg.state,\n",
        "    max_bond=mps_max_bond,\n",
        "    cutoff=mps_cutoff,\n",
        "    options=dict(maxiter=100),\n",
        ")\n",
        "t1 = timeit.default_timer()\n",
        "print(f\"Finished optimizing ground state ansatz in {t1 - t0} seconds.\")\n",
        "print(f\"Ground state ansatz fidelity: {1 - result.fun:.6f}\")\n",
        "\n",
        "gs_circuit = gs_ansatz.assign_parameters(result.x)\n",
        "gs_circuit_mps = quimb_circuit(\n",
        "    gs_circuit.decompose([\"PauliEvolution\"]),\n",
        "    quimb_circuit_class=qtn.CircuitMPS,\n",
        "    max_bond=mps_max_bond,\n",
        "    cutoff=mps_cutoff,\n",
        ")\n",
        "gs_ansatz_energy = qtn.expec_TN_1D(\n",
        "    gs_circuit_mps.psi.H, ham_mpo, gs_circuit_mps.psi\n",
        ")\n",
        "print(f\"Ground state ansatz energy: {gs_ansatz_energy:.6f}\")\n",
        "\n",
        "perturbed = gs_circuit.copy()\n",
        "perturbed.rz(np.pi / 2, center)\n",
        "\n",
        "circuits = []\n",
        "for t in range(1, n_steps + 1):\n",
        "    circuit = perturbed.copy()\n",
        "    for instr in trotter_evolution(\n",
        "        circuit.qubits, interaction, anisotropy, time_step, t\n",
        "    ):\n",
        "        circuit.append(instr)\n",
        "    circuits.append(circuit)\n",
        "print(\n",
        "    f\"Built {len(circuits)} circuits, deepest 2q depth (uniform basis) = \"\n",
        "    f\"{uniform_2q_depth(circuits[-1])}\"\n",
        ")\n",
        "\n",
        "observables = [\n",
        "    SparsePauliOp.from_sparse_list([(\"Z\", [i], 1)], num_qubits=n_qubits)\n",
        "    for i in range(n_qubits)\n",
        "]\n",
        "print(f\"Defined {len(observables)} Z observables.\")\n",
        "\n",
        "# ── Step 2: AQC ──────────────────────────────────────────────────────────────\n",
        "target_circuits = {\n",
        "    k: circuits[k - 1] for k in range(1, aqc_n_steps_total + 1)\n",
        "}\n",
        "\n",
        "aqc_sim = QuimbSimulator(\n",
        "    quimb_circuit_factory=partial(\n",
        "        qtn.CircuitMPS,\n",
        "        gate_opts=dict(max_bond=mps_max_bond, cutoff=mps_cutoff),\n",
        "    ),\n",
        "    autodiff_backend=\"jax\",\n",
        ")\n",
        "print(\"\\nStep 2b — target MPS:\")\n",
        "target_mps = {}\n",
        "for k in range(1, aqc_n_steps_total + 1):\n",
        "    target_mps[k] = tensornetwork_from_circuit(\n",
        "        target_circuits[k].decompose([\"PauliEvolution\"]), aqc_sim\n",
        "    )\n",
        "    print(f\"  k={k}: max bond = {target_mps[k].psi.max_bond()}\")\n",
        "\n",
        "ansatz_1, initial_params_1 = generate_ansatz_from_circuit(\n",
        "    target_circuits[1].decompose([\"PauliEvolution\"]),\n",
        "    qubits_initially_zero=True,\n",
        ")\n",
        "initial_params_1 = np.array(initial_params_1)\n",
        "ansatz_2, initial_params_2 = generate_ansatz_from_circuit(\n",
        "    target_circuits[2].decompose([\"PauliEvolution\"]),\n",
        "    qubits_initially_zero=True,\n",
        ")\n",
        "initial_params_2 = np.array(initial_params_2)\n",
        "print(\n",
        "    f\"\\nStep 2c — one-layer ansatz: {ansatz_1.num_parameters} parameters, \"\n",
        "    f\"2q depth (uniform basis) = {uniform_2q_depth(ansatz_1)}\"\n",
        ")\n",
        "print(\n",
        "    f\"Step 2c — two-layer ansatz: {ansatz_2.num_parameters} parameters, \"\n",
        "    f\"2q depth (uniform basis) = {uniform_2q_depth(ansatz_2)}\"\n",
        ")\n",
        "\n",
        "aqc_circuits = {}\n",
        "aqc_params = {}\n",
        "for k in range(1, aqc_n_steps_total + 1):\n",
        "    if k <= aqc_n_steps_1:\n",
        "        ansatz, base_params = ansatz_1, initial_params_1\n",
        "    else:\n",
        "        ansatz, base_params = ansatz_2, initial_params_2\n",
        "    # Warm-start from the previous step only within the same stage\n",
        "    same_stage = (k - 1 >= 1) and (\n",
        "        (k - 1 <= aqc_n_steps_1) == (k <= aqc_n_steps_1)\n",
        "    )\n",
        "    x0 = aqc_params[k - 1] if same_stage else base_params\n",
        "    obj = MaximizeStateFidelity(target_mps[k], ansatz, aqc_sim)\n",
        "    t0 = timeit.default_timer()\n",
        "    result = scipy.optimize.minimize(\n",
        "        obj.loss_function,\n",
        "        x0,\n",
        "        method=\"L-BFGS-B\",\n",
        "        jac=True,\n",
        "        options=dict(maxiter=100),\n",
        "    )\n",
        "    elapsed = timeit.default_timer() - t0\n",
        "    aqc_params[k] = result.x\n",
        "    aqc_circuits[k] = ansatz.assign_parameters(result.x)\n",
        "    print(\n",
        "        f\"  k={k}: fidelity = {1 - result.fun:.4f}, \"\n",
        "        f\"2q depth (uniform basis) = \"\n",
        "        f\"{uniform_2q_depth(aqc_circuits[k])}, \"\n",
        "        f\"{elapsed:.1f}s\"\n",
        "    )\n",
        "\n",
        "all_circuits = []\n",
        "for k in range(1, aqc_n_steps_total + 1):\n",
        "    all_circuits.append(aqc_circuits[k])\n",
        "base = aqc_circuits[aqc_n_steps_total]\n",
        "for k in range(1, n_steps - aqc_n_steps_total + 1):\n",
        "    circuit = base.copy()\n",
        "    for instr in trotter_evolution(\n",
        "        circuit.qubits, interaction, anisotropy, time_step, k\n",
        "    ):\n",
        "        circuit.append(instr)\n",
        "    all_circuits.append(circuit)\n",
        "\n",
        "full_depths = [\n",
        "    uniform_2q_depth(circuits[k - 1]) for k in range(1, n_steps + 1)\n",
        "]\n",
        "aqc_depths = [uniform_2q_depth(circuit) for circuit in all_circuits]\n",
        "full_2q = [\n",
        "    circuits[k - 1].decompose([\"PauliEvolution\"]).num_nonlocal_gates()\n",
        "    for k in range(1, n_steps + 1)\n",
        "]\n",
        "aqc_2q = [\n",
        "    circuit.decompose([\"PauliEvolution\"]).num_nonlocal_gates()\n",
        "    for circuit in all_circuits\n",
        "]\n",
        "print(f\"\\nStep 2d — assembled {len(all_circuits)} circuits\")\n",
        "print(\n",
        "    f\"  At step {n_steps} (uniform basis): \"\n",
        "    f\"Trotter 2q depth = {full_depths[-1]}, AQC+Trotter 2q depth = {aqc_depths[-1]}; \"\n",
        "    f\"2q gates: Trotter = {full_2q[-1]}, AQC+Trotter = {aqc_2q[-1]}\"\n",
        ")\n",
        "\n",
        "steps_axis = np.arange(1, n_steps + 1)\n",
        "fig, ax = plt.subplots(figsize=(7, 4))\n",
        "ax.plot(steps_axis, full_depths, \"-o\", color=\"black\", label=\"Full Trotter\")\n",
        "ax.plot(\n",
        "    steps_axis, aqc_depths, \"-o\", color=\"cadetblue\", label=\"AQC + Trotter\"\n",
        ")\n",
        "ax.set_xlabel(\"Trotter steps\", fontsize=13)\n",
        "ax.xaxis.set_major_locator(plt.matplotlib.ticker.MaxNLocator(integer=True))\n",
        "ax.set_ylabel(\"2q gate depth (uniform basis)\", fontsize=13)\n",
        "ax.set_title(f\"AQC circuit-depth reduction ({n_qubits} qubits)\", fontsize=13)\n",
        "ax.legend(fontsize=11)\n",
        "plt.tight_layout()\n",
        "plt.show()\n",
        "\n",
        "# ── Step 3: Execute on IBM Quantum hardware ───────────────────────────────────\n",
        "# (replaces StatevectorEstimator)\n",
        "service = QiskitRuntimeService()\n",
        "backend = service.least_busy(\n",
        "    min_num_qubits=n_qubits,\n",
        "    operational=True,\n",
        "    simulator=False,\n",
        "    filters=lambda x: x.configuration().processor_type[\"family\"] == \"Heron\",\n",
        ")\n",
        "print(f\"Backend: {backend.name}\")\n",
        "\n",
        "pm = generate_preset_pass_manager(optimization_level=3, backend=backend)\n",
        "isa_circuits = pm.run(all_circuits, num_processes=1)\n",
        "isa_2q_depths = [\n",
        "    isa_circuit.depth(lambda inst: inst.operation.num_qubits == 2)\n",
        "    for isa_circuit in isa_circuits\n",
        "]\n",
        "print(\n",
        "    f\"Transpiled 2q depth (deepest, ISA on {backend.name}): \"\n",
        "    f\"{max(isa_2q_depths)} \"\n",
        "    f\"(2q gates: {max(isa_circuit.num_nonlocal_gates() for isa_circuit in isa_circuits)})\"\n",
        ")\n",
        "\n",
        "estimator = Estimator(backend)\n",
        "estimator.options.environment.job_tags = [\"TUT_SNS\"]\n",
        "estimator.options.dynamical_decoupling.enable = True\n",
        "estimator.options.dynamical_decoupling.sequence_type = \"XY4\"\n",
        "estimator.options.twirling.enable_gates = True\n",
        "estimator.options.twirling.num_randomizations = 1000\n",
        "estimator.options.twirling.shots_per_randomization = 128\n",
        "estimator.options.resilience.measure_mitigation = True\n",
        "estimator.options.resilience.measure_noise_learning.num_randomizations = 32\n",
        "estimator.options.resilience.measure_noise_learning.shots_per_randomization = 100\n",
        "\n",
        "pubs = [\n",
        "    (\n",
        "        isa_circuit,\n",
        "        [obs.apply_layout(isa_circuit.layout) for obs in observables],\n",
        "    )\n",
        "    for isa_circuit in isa_circuits\n",
        "]\n",
        "job = estimator.run(pubs)\n",
        "print(f\"Job ID: {job.job_id()}\")\n",
        "\n",
        "result = job.result()\n",
        "\n",
        "# ── Step 4: Post-process ──────────────────────────────────────────────────────\n",
        "rgf_mat = np.stack([pub_result.data.evs for pub_result in result])\n",
        "\n",
        "n_points_momentum, n_points_frequency = 100, 100\n",
        "spectrum = get_dsf(\n",
        "    n_qubits,\n",
        "    rgf_mat,\n",
        "    time_step,\n",
        "    n_steps,\n",
        "    n_points_momentum,\n",
        "    n_points_frequency,\n",
        ")\n",
        "\n",
        "plot_rgf(\n",
        "    n_qubits,\n",
        "    rgf_mat,\n",
        "    time_step,\n",
        "    n_steps,\n",
        "    title=f\"Retarded Green's function — {n_qubits} qubits (QPU)\",\n",
        ")\n",
        "plot_dsf(\n",
        "    spectrum,\n",
        "    time_step,\n",
        "    n_points_momentum,\n",
        "    n_points_frequency,\n",
        "    title=rf\"KCuF$_3$ DSF — {n_qubits} qubits (QPU)\"\n",
        "    \"\\n(AQC + DD + Pauli twirling + TREX)\",\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "ls-discussion",
      "metadata": {},
      "source": [
        "Die Hardware-Ergebnisse geben die wesentlichen Merkmale des Zwei-Spinon-Kontinuums wieder: Die Streuintensität konzentriert sich bei niedrigen Energien in der Nähe des antiferromagnetischen Wellenvektors $q = \\pi$ und wird nach unten durch die sinusförmige Spinon-Dispersion begrenzt, wobei darüber ein breites Kontinuum an spektralem Gewicht liegt und nicht nur ein einzelner scharfer Modus. Dies ist dieselbe Struktur, die mittels unelastischer Neutronenstreuung an KCuF$_3$ gemessen wurde, und sie bestätigt den gesamten Arbeitsablauf – von der Grundzustandspräparation über die Störung bis hin zur AQC-komprimierten Trotter-Evolution und der fehlermindernden Messung – im 50-Qubit-Maßstab.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "ls-nextsteps",
      "metadata": {},
      "source": [
        "<span id=\"next-steps\" />\n",
        "\n",
        "## Nächste Schritte\n",
        "\n",
        "Wenn Sie diese Arbeit interessant fanden, könnten Sie auch an folgenden Materialien interessiert sein:\n",
        "\n",
        "<Admonition type=\"tip\" title=\"Empfehlungen\">\n",
        "  * [Simulation der Neutronenstreuung mit einem AQC- und Trotter-Dynamik-basierten serverlosen Workflow](/docs/tutorials/simulate-neutron-scattering-with-a-serverless-workflow) – ein begleitendes Tutorial, in dem dasselbe Experiment mithilfe einer bereitgestellten „ Qiskit Serverless “-Funktionsvorlage durchgeführt wird\n",
        "  * [Lee et al., „Benchmarking von Quantensimulationen anhand von Neutronenstreuungsexperimenten“ (arXiv:2603.15608 )](https://arxiv.org/abs/2603.15608) – die Referenzarbeit, auf der dieses Tutorial basiert\n",
        "  * [Techniken zur Fehlerminderung und -unterdrückung](/docs/guides/error-mitigation-and-suppression-techniques) – DD, Pauli-Twirling und TREX, die in den Hardware-Experimenten zum Einsatz kamen\n",
        "  * [Ungefähre Quantenkompilierung für Zeitentwicklungsschaltungen](/docs/tutorials/approximate-quantum-compilation-for-time-evolution) – Tutorial zu AQC-Tensor\n",
        "  * [AQC-Tensor-Dokumentation](https://qiskit.github.io/qiskit-addon-aqc-tensor/)\n",
        "</Admonition>\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "id": "a1b8767d",
      "source": "© IBM Corp., 2017-2026"
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3"
    },
    "hours": 1.5,
    "qpuSeconds": 780
  },
  "nbformat": 4,
  "nbformat_minor": 5
}