diff --git a/docs/tutorials/_toc.json b/docs/tutorials/_toc.json
index b99c2459c95..0afaaab4c38 100644
--- a/docs/tutorials/_toc.json
+++ b/docs/tutorials/_toc.json
@@ -53,6 +53,10 @@
         {
           "title": "Observable estimation",
           "children": [
+            {
+              "title": "Simulate neutron scattering in quantum materials with quantum circuits",
+              "url": "/docs/tutorials/simulate-neutron-scattering"
+            },
             {
               "title": "Krylov quantum diagonalization of lattice Hamiltonians",
               "url": "/docs/tutorials/krylov-quantum-diagonalization"
diff --git a/docs/tutorials/index.mdx b/docs/tutorials/index.mdx
index 38e9cb9a8f1..3d70cba8ea3 100644
--- a/docs/tutorials/index.mdx
+++ b/docs/tutorials/index.mdx
@@ -47,6 +47,8 @@ The tutorials highlight techniques where repeated sampling enables estimation of
 
 These tutorials focus on estimating physically meaningful quantities, such as energy or correlation values, by preparing quantum states and measuring observables. Techniques include both variational and Trotterized circuit approaches that balance circuit expressiveness with circuit-depth efficiency. Emphasis is placed on workflows that reduce quantum resource demands while maintaining accuracy, and enabling practical estimation of observables in chemical and physical systems.
 
+* [Simulate neutron scattering in quantum materials with quantum circuits](/docs/tutorials/simulate-neutron-scattering)
+
 * [Krylov quantum diagonalization of lattice Hamiltonians](/docs/tutorials/krylov-quantum-diagonalization)
 
 * [Nishimori phase transition](/docs/tutorials/nishimori-phase-transition)
diff --git a/docs/tutorials/simulate-neutron-scattering.ipynb b/docs/tutorials/simulate-neutron-scattering.ipynb
new file mode 100644
index 00000000000..60ae711d46c
--- /dev/null
+++ b/docs/tutorials/simulate-neutron-scattering.ipynb
@@ -0,0 +1,1340 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "8e82ead7",
+   "metadata": {},
+   "source": [
+    "---\n",
+    "title: Simulate neutron scattering in quantum materials with quantum circuits\n",
+    "description: Compute the dynamical structure factor of a quantum magnet using Trotter circuits and MPS simulation.\n",
+    "---\n",
+    "\n",
+    "{/* cspell:ignore DSFs spinon DMRG viridis fontsize vmax vmin cadetblue antiferromagnet COBYQA */}\n",
+    "\n",
+    "# Simulate neutron scattering in quantum materials with quantum circuits\n",
+    "\n",
+    "*Usage estimate: 13 minutes on a Heron r2 processor (NOTE: This is an estimate only. Your runtime might vary.)*"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "11033666",
+   "metadata": {},
+   "source": [
+    "## Learning outcomes\n",
+    "\n",
+    "After completing this tutorial, you can expect to understand the following information:\n",
+    "- How inelastic neutron scattering (INS) spectra connect to dynamical structure factors (DSFs) of quantum spin models.\n",
+    "- How to prepare a ground state, apply a local perturbation, and perform Trotter time evolution on a quantum circuit.\n",
+    "- How to use approximate quantum compiling (AQC) with `qiskit-addon-aqc-tensor` to compress deep Trotter circuits for hardware execution.\n",
+    "- How to extract the retarded Green's function (RGF) from qubit expectation values and Fourier-transform it into a DSF.\n",
+    "\n",
+    "## Prerequisites\n",
+    "It is recommended that you familiarize yourself with these topics:\n",
+    "- [Basics of quantum information](/learning/courses/basics-of-quantum-information)\n",
+    "- [Variational algorithm design](/learning/courses/variational-algorithm-design)\n",
+    "- [Introduction to Qiskit Primitives (Estimator and Sampler)](/docs/guides/qiskit-runtime-primitives)\n",
+    "\n",
+    "## Background\n",
+    "\n",
+    "In this tutorial, we reproduce the results of [Lee et al., arXiv:2603.15608](https://arxiv.org/abs/2603.15608).\n",
+    "\n",
+    "### Inelastic neutron scattering and the dynamical structure factor\n",
+    "\n",
+    "Inelastic neutron scattering (INS) is one of the most powerful experimental probes of magnetic excitations in quantum materials. When a beam of thermal or cold neutrons impinges on a crystal, individual neutrons exchange both momentum $\\mathbf{q}$ and energy $\\omega$ with the magnetic subsystem. The measured scattering intensity is proportional to the dynamical structure factor (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",
+    "which encodes the full space-time correlations of the spin degrees of freedom.\n",
+    "\n",
+    "### KCuF$_3$: a canonical Luttinger-liquid magnet\n",
+    "\n",
+    "Potassium copper fluoride (KCuF$_3$) is a quasi-one-dimensional antiferromagnet in which chains of spin-$\\frac{1}{2}$ Cu$^{2+}$ ions interact via a nearest-neighbor Heisenberg exchange $J$, while the interchain coupling is only $\\sim 2.7\\%$ of $J$. At $T = 6\\;\\mathrm{K}$ where INS data are available, the spectrum is dominated by fractionalized spinon excitations characteristic of a Tomonaga-Luttinger liquid. Because the intrachain dynamics are well described by the one-dimensional spin-$\\frac{1}{2}$ XXZ Hamiltonian at the isotropic point ($\\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$ serves as an ideal benchmark for quantum simulation: the Hamiltonian is simple enough to implement on a quantum processor, yet the ground state is strongly entangled and the excitation spectrum features a broad two-spinon continuum.\n",
+    "\n",
+    "Note: this tutorial sets $J = 1$ as the energy unit and adopts the normalization $H = J\\sum_i[\\ldots]$, which corresponds to the local Hamiltonian $H_\\mathrm{loc} = J(S^X S^X + S^Y S^Y + S^Z S^Z)$ used in the paper's circuit implementation (Fig. S3 of the supplementary). The paper's full Hamiltonian (Eq. 3) carries an additional overall factor of 2, so the paper's $J$ is twice the $J$ used here.\n",
+    "\n",
+    "### What we simulate and measure\n",
+    "\n",
+    "The physical quantity we compute is the **retarded Green's function** (RGF), defined as the time-dependent spin-spin correlation function\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",
+    "where $j_c$ is a reference site (the chain center) and $S_j^\\alpha(t) = e^{iHt}S_j^\\alpha e^{-iHt}$ is the Heisenberg-picture spin operator. In this tutorial we focus on the $zz$ component ($\\alpha = \\beta = z$). On a quantum computer the RGF is accessed by preparing the ground state, applying a local perturbation at $j_c$, time-evolving the perturbed state, and measuring the single-qubit expectation value $\\langle\\sigma_j^z\\rangle$ at every site $j$ for each time step. The key insight is that each $\\langle\\sigma_j^z\\rangle$ gives the difference from the ground-state magnetization; because the isotropic Heisenberg antiferromagnet has zero net per-site magnetization ($\\langle\\sigma_j^z\\rangle_\\mathrm{GS} = 0$), the raw measured value directly yields $G^R(j, j_c, t)$ without any explicit subtraction.\n",
+    "\n",
+    "By collecting $G^R(j, j_c, t)$ over all sites and time steps, we assemble a two-dimensional dataset that is then Fourier-transformed in both space and time to produce the **dynamical structure factor** $S(q,\\omega)$. The DSF is the quantity directly measured in an INS experiment: it tells us what magnetic excitations exist at each momentum $q$ and energy $\\omega$. For the isotropic Heisenberg chain the exact excitation spectrum is a **two-spinon continuum**, a broad band of scattering intensity whose shape serves as a stringent end-to-end benchmark for the quantum simulation: it validates the ground-state preparation, perturbation, Trotter time evolution, and measurement protocol all at once.\n",
+    "\n",
+    "### Quantum-simulation workflow\n",
+    "\n",
+    "The workflow mirrors the physics of an INS event. We (1) prepare the many-body ground state $|\\psi_{\\mathrm{GS}}\\rangle$ on $n$ qubits, (2) apply a local spin-flip perturbation $U_{j_c} = \\frac{1}{\\sqrt{2}}(I - i\\sigma^z_{j_c})$ at the chain center to mimic the neutron's spin transfer, (3) evolve under $H$ for discrete time steps using second-order Trotterization, and (4) measure $\\langle\\sigma_i^z\\rangle$ on every qubit at each step to obtain the RGF. A two-dimensional discrete Fourier transform then yields the DSF $S(q,\\omega)$.\n",
+    "\n",
+    "The observable we measure at every time step is $\\sigma_i^z$ on each qubit $i$. In Qiskit this is represented as a list of `SparsePauliOp` operators: one single-qubit $Z$ embedded in the $n$-qubit identity string for each site. These observables are constructed once per problem instance during the problem-mapping step (Step 1) and reused for every circuit at that scale.\n",
+    "\n",
+    "### Approximate quantum compiling (AQC)\n",
+    "\n",
+    "Deep Trotter circuits can be compressed with **approximate quantum compiling (AQC)**, which replaces the first several Trotter layers with a shorter parameterized ansatz whose parameters are classically optimized to maximize the MPS-level fidelity with the original deep circuit. The remaining Trotter steps are appended exactly, producing a \"mixed\" AQC + Trotter circuit with substantially fewer two-qubit gates.\n",
+    "\n",
+    "### MPS simulation\n",
+    "\n",
+    "For a one-dimensional system, matrix-product-state (MPS) methods can efficiently simulate both ground-state preparation (with the density matrix renormalization group, or DMRG) and circuit-level time evolution. By controlling the bond dimension $\\chi$, we trade off accuracy against computational cost. In this tutorial we use MPS simulation within `qiskit-addon-aqc-tensor` to compute high-fidelity AQC ansätze that compress the deep Trotter circuits for hardware execution."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "de1ad5a2",
+   "metadata": {},
+   "source": [
+    "## Requirements\n",
+    "Before starting this tutorial, be sure you have the following installed:\n",
+    "\n",
+    "- Qiskit SDK with [visualization](/docs/api/qiskit/visualization) support\n",
+    "- Qiskit Runtime (`pip install qiskit-ibm-runtime`)\n",
+    "- `qiskit-addon-aqc-tensor` with `quimb` and JAX extras (`pip install 'qiskit-addon-aqc-tensor[quimb-jax]'`)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "9938e4bd",
+   "metadata": {},
+   "source": [
+    "## Setup"
+   ]
+  },
+  {
+   "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": [
+    "## Small-scale simulator example\n",
+    "\n",
+    "We first demonstrate the full workflow on **10 qubits**, optimizing the HVA ground-state ansatz by using MPS simulation and using the Qiskit statevector simulator for time evolution. DMRG gives a reference ground-state energy and a reference MPS. This small-scale example lets us validate every step before scaling up."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "3818ed8f",
+   "metadata": {},
+   "source": [
+    "### Step 1: Map classical inputs to a quantum problem\n",
+    "\n",
+    "We begin by defining the physical model and constructing the quantum circuits.\n",
+    "\n",
+    "**Hamiltonian.** KCuF$_3$ is modeled by the 1D XXZ Hamiltonian at the isotropic point ($J = 1$, $\\epsilon = 1$, setting $J = 1$ as the energy unit).\n",
+    "\n",
+    "**Ground state.** We use a Hamiltonian variational ansatz (HVA) circuit, `build_ground_state_ansatz`, as the ground-state preparation circuit. The HVA parameters are optimized classically by maximizing the state fidelity $|\\langle\\psi_{\\mathrm{HVA}}(\\theta)|\\psi_{\\mathrm{DMRG}}\\rangle|^2$ with the DMRG ground-state MPS, where the HVA state is evaluated by simulating the circuit as a matrix product state with quimb's `CircuitMPS`. We use `scipy.optimize.minimize` for the optimization. The energy $\\langle H \\rangle$ of the optimized ansatz is also computed for reference.\n",
+    "\n",
+    "**Trotter gates.** Each nearest-neighbor interaction term $e^{-i\\Delta t\\, H_{\\mathrm{pair}}}$ with $H_{\\mathrm{pair}} = (J/4)(XX + YY + ZZ)$ is constructed with `PauliEvolutionGate(H_pair, time=time_step)`. Qiskit synthesizes this into the optimal three-CNOT decomposition during transpilation.\n",
+    "\n",
+    "**Perturbation.** An $R_z(\\pi/2)$ gate on the center qubit implements $U_{j_c} = \\frac{1}{\\sqrt{2}}(I - i\\sigma^z_{j_c})$, mimicking the local spin-flip produced by a scattered neutron.\n",
+    "\n",
+    "**Observables.** We construct a $\\sigma^z$ observable for each qubit site. These `SparsePauliOp` objects are passed to the Estimator primitive in Step 3 to extract $\\langle\\sigma_i^z\\rangle$ at every time step."
+   ]
+  },
+  {
+   "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": [
+    "### Step 2: Optimize problem for quantum hardware execution\n",
+    "\n",
+    "For real hardware, the Trotter circuits above would be too deep. **Approximate quantum compiling (AQC)** addresses this by replacing the first $k$ Trotter layers (including the ground-state circuit) with a shorter parameterized ansatz optimized to maximize the MPS-level fidelity with the original deep circuit. The remaining Trotter steps are appended exactly, producing a shallower \"AQC + Trotter\" circuit.\n",
+    "\n",
+    "To balance expressibility against circuit depth, two ansätze are used: a **one-layer ansatz** (generated from a single Trotter step) compresses the earliest time steps, and a deeper **two-layer ansatz** (generated from two Trotter steps) compresses the next few steps where higher fidelity is needed.\n",
+    "\n",
+    "The AQC workflow has four sub-steps:\n",
+    "\n",
+    "1. **Build target circuits** &mdash; the first $k_1 + k_2$ circuits from Step 1 serve directly as AQC targets.\n",
+    "2. **Compute target MPS** &mdash; simulate each target circuit as a matrix-product state using `quimb.tensor.CircuitMPS`.\n",
+    "3. **Generate and optimize ansätze** &mdash; `generate_ansatz_from_circuit` creates a one-layer and a two-layer parameterized ansatz; parameters are optimized by using L-BFGS-B with JAX-accelerated gradients to minimize $1 - |\\langle\\psi_{\\mathrm{ansatz}}|\\psi_{\\mathrm{target}}\\rangle|^2$. The first $k_1$ steps use the one-layer ansatz and the next $k_2$ steps use the two-layer ansatz; within each stage, each step warm-starts from the previous step's optimized parameters, and parameters are reset to the stage defaults at the stage boundary.\n",
+    "4. **Assemble mixed circuits** &mdash; for time steps beyond the AQC checkpoint, append exact Trotter layers to the optimized two-layer AQC circuit using `trotter_evolution`."
+   ]
+  },
+  {
+   "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": [
+    "### Step 3: Execute using Qiskit primitives\n",
+    "\n",
+    "We simulate each AQC-compiled circuit using the `StatevectorEstimator` primitive."
+   ]
+  },
+  {
+   "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": [
+    "### Step 4: Post-process and return result in desired classical format\n",
+    "\n",
+    "We now extract the expectation value $\\langle\\sigma_i^z\\rangle$ for every qubit $i$ at each time step $t$. These values form the retarded Green's function matrix $G^R(j, j_c, t)$. Then, we Fourier-transform the RGF into the dynamical structure factor $S(q, \\omega)$ and plot both the RGF and the DSF. `get_dsf` applies mirror symmetry and clips negative values internally."
+   ]
+  },
+  {
+   "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": [
+    "## Large-scale hardware execution\n",
+    "\n",
+    "We now scale up to **50 qubits**. At this scale, the optimizations reach lower fidelities than in the small-scale example: the ground-state ansatz fidelity drops to roughly 0.65, and the AQC fidelities decrease to roughly 0.7 at the later checkpoints. This is expected, and you can improve these fidelities by increasing the number of ground-state ansatz layers (`gs_n_layers`) or the number of optimization iterations (`maxiter`) at additional classical cost. Note also that the two-layer AQC fidelities come out lower than the one-layer ones. This is not a regression: later time steps generate more entanglement and are simply harder to compress, which is why the more expressive two-layer ansatz is used for them.\n",
+    "\n",
+    "The AQC optimization can also take several hours of classical compute time (roughly six hours in the run shown here, most of it spent on the two-layer checkpoints in Step 2c). To reduce the wall-clock time, consider running this notebook on more powerful classical hardware, such as a high-performance computing (HPC) system. Alternatively, you can scale down to a smaller problem instance (for example, fewer qubits or fewer time steps); your results will then differ from the ones shown here.\n",
+    "\n",
+    "The code below follows the same four-step structure as the small-scale example. The HVA ground-state parameters are again optimized by using MPS simulation. On the QPU, we enable dynamical decoupling (DD), Pauli twirling, and twirled readout error extinction (TREX) for error suppression and mitigation. The following table summarizes how the large-scale experiment differs from the small-scale one:\n",
+    "\n",
+    "| &nbsp; | Small scale | Large scale |\n",
+    "|---|---|---|\n",
+    "| Qubits | 10 | 50 |\n",
+    "| Time steps | 10 | 20 |\n",
+    "| AQC checkpoints (1-layer + 2-layer) | 3 + 2 = 5 | 6 + 4 = 10 |\n",
+    "| Ground state ansatz layers | 3 | 5 |\n",
+    "| MPS max bond dimension | 32 | 128 |\n",
+    "| Estimator | `StatevectorEstimator` | QPU with DD, Pauli twirling, and TREX |"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "9e85fbd1",
+   "metadata": {},
+   "source": [
+    "<Admonition type=\"note\" title=\"XLA slow-compile messages\">\n",
+    "\n",
+    "During the AQC optimization (Step 2c below) you may see stderr messages such as:\n",
+    "\n",
+    "```\n",
+    "[Compiling module jit_MakeArrayFn for CPU] Very slow compile? ...\n",
+    "The operation took 2m14s\n",
+    "```\n",
+    "\n",
+    "These are benign diagnostics from XLA, the compiler that backs `qiskit-addon-aqc-tensor`'s JAX autodiff. At 50 qubits with MPS bond dimension 128, XLA spends a couple of minutes compiling the gradient function the first time it is traced. The compilation succeeds, and the optimization results are unaffected.\n",
+    "\n",
+    "</Admonition>"
+   ]
+  },
+  {
+   "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": [
+    "The hardware results reproduce the key features of the two-spinon continuum: the scattering intensity is concentrated near the antiferromagnetic wave vector $q = \\pi$ at low energy and is bounded from below by the sinusoidal spinon dispersion, with a broad continuum of spectral weight above it rather than a single sharp mode. This is the same structure measured by inelastic neutron scattering on KCuF$_3$, and it validates the full workflow of ground-state preparation, perturbation, AQC-compressed Trotter evolution, and error-mitigated measurement at the 50-qubit scale."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "ls-nextsteps",
+   "metadata": {},
+   "source": [
+    "## Next steps\n",
+    "If you found this work interesting, you might be interested in the following material:\n",
+    "<Admonition type=\"tip\" title=\"Recommendations\">\n",
+    "\n",
+    "- [Lee et al., \"Benchmarking quantum simulation with neutron-scattering experiments\" (arXiv:2603.15608)](https://arxiv.org/abs/2603.15608) &mdash; the reference paper this tutorial is based on\n",
+    "- [Error mitigation and suppression techniques](/docs/guides/error-mitigation-and-suppression-techniques) &mdash; DD, Pauli twirling, and TREX used in the hardware experiments\n",
+    "- [Approximate quantum compilation for time evolution circuits](/docs/tutorials/approximate-quantum-compilation-for-time-evolution) &mdash; tutorial on AQC-Tensor\n",
+    "- [AQC-Tensor documentation](https://qiskit.github.io/qiskit-addon-aqc-tensor/)\n",
+    "\n",
+    "</Admonition>"
+   ]
+  }
+ ],
+ "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"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/public/docs/images/tutorials/simulate-neutron-scattering/extracted-outputs/4585d139-0.avif b/public/docs/images/tutorials/simulate-neutron-scattering/extracted-outputs/4585d139-0.avif
new file mode 100644
index 00000000000..185560dee37
Binary files /dev/null and b/public/docs/images/tutorials/simulate-neutron-scattering/extracted-outputs/4585d139-0.avif differ
diff --git a/public/docs/images/tutorials/simulate-neutron-scattering/extracted-outputs/4585d139-1.avif b/public/docs/images/tutorials/simulate-neutron-scattering/extracted-outputs/4585d139-1.avif
new file mode 100644
index 00000000000..20afe5af1b1
Binary files /dev/null and b/public/docs/images/tutorials/simulate-neutron-scattering/extracted-outputs/4585d139-1.avif differ
diff --git a/public/docs/images/tutorials/simulate-neutron-scattering/extracted-outputs/e78f8d9a-1.avif b/public/docs/images/tutorials/simulate-neutron-scattering/extracted-outputs/e78f8d9a-1.avif
new file mode 100644
index 00000000000..c94b9eba108
Binary files /dev/null and b/public/docs/images/tutorials/simulate-neutron-scattering/extracted-outputs/e78f8d9a-1.avif differ
diff --git a/public/docs/images/tutorials/simulate-neutron-scattering/extracted-outputs/fbe57e1a-5.avif b/public/docs/images/tutorials/simulate-neutron-scattering/extracted-outputs/fbe57e1a-5.avif
new file mode 100644
index 00000000000..980cd01f94f
Binary files /dev/null and b/public/docs/images/tutorials/simulate-neutron-scattering/extracted-outputs/fbe57e1a-5.avif differ
diff --git a/public/docs/images/tutorials/simulate-neutron-scattering/extracted-outputs/fbe57e1a-7.avif b/public/docs/images/tutorials/simulate-neutron-scattering/extracted-outputs/fbe57e1a-7.avif
new file mode 100644
index 00000000000..3d235ef8e31
Binary files /dev/null and b/public/docs/images/tutorials/simulate-neutron-scattering/extracted-outputs/fbe57e1a-7.avif differ
diff --git a/public/docs/images/tutorials/simulate-neutron-scattering/extracted-outputs/fbe57e1a-8.avif b/public/docs/images/tutorials/simulate-neutron-scattering/extracted-outputs/fbe57e1a-8.avif
new file mode 100644
index 00000000000..6c06f96c85c
Binary files /dev/null and b/public/docs/images/tutorials/simulate-neutron-scattering/extracted-outputs/fbe57e1a-8.avif differ
diff --git a/qiskit_bot.yaml b/qiskit_bot.yaml
index 60531ff8260..d58a69b1260 100644
--- a/qiskit_bot.yaml
+++ b/qiskit_bot.yaml
@@ -795,6 +795,9 @@ notifications:
     - "`@nathanearnestnoble`"
     - "@murilo-kipu"
     - "@HuangJunye"
+  "docs/tutorials/simulate-neutron-scattering":
+    - "`@nathanearnestnoble`"
+    - "@kevinsung"
   "docs/tutorials/readout-error-mitigation-sampler":
     - "`@nathanearnestnoble`"
     - "`@jlapeyre`"
diff --git a/scripts/config/notebook-testing.toml b/scripts/config/notebook-testing.toml
index 3ab287b0d5c..5add6d979c8 100644
--- a/scripts/config/notebook-testing.toml
+++ b/scripts/config/notebook-testing.toml
@@ -207,6 +207,7 @@ notebooks = [
   "docs/tutorials/edc-cut-bell-pair-benchmarking.ipynb",
   "docs/tutorials/compilation-methods-for-hamiltonian-simulation-circuits.ipynb",
   "docs/tutorials/solve-market-split-problem-with-iskay-quantum-optimizer.ipynb",
+  "docs/tutorials/simulate-neutron-scattering.ipynb",
 
   # Don't test any learning notebooks
   "learning/courses/quantum-computing-in-practice/introduction.ipynb",