{ "cells": [ { "cell_type": "markdown", "id": "frontmatter", "metadata": {}, "source": [ "---\n", "title: Simulate neutron scattering with an AQC + Trotter dynamics Serverless workflow\n", "description: Compute the dynamical structure factor S(q, w) of the quantum magnet KCuF3 by running a Trotter workflow with AQC compression as a deployed function template.\n", "---\n", "\n", "{/* cspell:ignore Trotter Trotterization spinon spinons KCuF quimb DMRG magnon antiferromagnetic Suzuki fidelities isa COBYQA Gjjc qpoints viridis fontsize vmax vmin */}\n", "\n", "# Simulate neutron scattering with an AQC + Trotter dynamics Serverless workflow\n", "\n", "*Usage estimate: 18 minutes on a Heron r3 processor (NOTE: This is an estimate only. Your runtime might vary.)*\n", "\n" ] }, { "cell_type": "markdown", "id": "learning-outcomes", "metadata": {}, "source": [ "## Learning outcomes\n", "\n", "* How an inelastic neutron-scattering spectrum maps to the dynamical structure factor $S(q, \\omega)$ of a 1D quantum magnet.\n", "* How to prepare the KCuF$_3$ (isotropic Heisenberg) ground state with the density matrix renormalization group (DMRG) and matrix product state (MPS) fidelity maximization.\n", "* How to run Trotter time-evolution, approximate quantum compilation (AQC) circuit compression, and mitigated execution as a single function call.\n", "* How to post-process the per-site $\\langle \\sigma_z \\rangle(t)$ time series into $S(q, \\omega)$ and identify the two-spinon continuum.\n", "\n" ] }, { "cell_type": "markdown", "id": "prerequisites", "metadata": {}, "source": [ "## Prerequisites\n", "\n", "* Familiarity with [Qiskit patterns](/docs/guides/intro-to-patterns), [`SparsePauliOp`](/docs/api/qiskit/qiskit.quantum_info.SparsePauliOp), and [Trotter time-evolution](/learning/courses/utility-scale-quantum-computing/quantum-simulation).\n", "* Basic exposure to tensor-network methods (DMRG and MPS) is helpful but not required, as is familiarity with the [`qiskit-addon-aqc-tensor`](https://github.com/Qiskit/qiskit-addon-aqc-tensor) library that the function uses to compress Trotter circuits.\n", "\n" ] }, { "cell_type": "markdown", "id": "background", "metadata": {}, "source": [ "## Background\n", "\n", "Inelastic neutron scattering measures the dynamical structure factor $S(q, \\omega)$, the space-and-time Fourier transform of the spin-spin correlation function, so reproducing $S(q, \\omega)$ from a microscopic spin model is a direct, falsifiable test of a quantum simulation. This tutorial studies KCuF$_3$, a spin-$\\frac{1}{2}$ antiferromagnetic Heisenberg chain whose excitations are not single spin flips but pairs of fractionalized spinons: instead of a sharp magnon dispersion, $S(q, \\omega)$ shows a broad *two-spinon continuum*, bounded below by $\\tfrac{\\pi}{2}|\\sin q|$ and above by $\\pi|\\sin(q/2)|$. Those are the dashed curves on the plots that follow. The physics in full, and the comparison against measured neutron data, are covered in the [original tutorial](/docs/tutorials/simulate-neutron-scattering) and in Lee et al., [arXiv:2603.15608](https://arxiv.org/abs/2603.15608).\n", "\n", "The quantum workflow mirrors the scattering experiment:\n", "\n", "1. Prepare the chain's ground state $|\\psi_0\\rangle$.\n", "2. Kick it with a local perturbation at the center site, a $\\pi/2$ $Z$-rotation, mimicking the momentum and energy transfer from the neutron.\n", "3. Time-evolve under the Heisenberg Hamiltonian, $e^{-iHt}$, with a Trotter product formula.\n", "4. Measure the per-site magnetization $\\langle \\sigma_z^j \\rangle(t)$. As a function of site $j$ and time $t$, this is exactly the retarded Green's function $G^R(j, j_c, t)$, so no conversion is needed before the Fourier transform in step 5.\n", "5. Fourier transform $G^R$ into $S(q, \\omega)$.\n", "\n", "Problems can arise in step 3, when exact Trotter circuits for long evolutions become too deep for hardware. AQC with tensor networks addresses this by compressing a block of Trotter steps into a fixed, shallow parameterized ansatz whose state fidelity to the exact evolution is maximized classically with an MPS simulator ([arXiv:2301.08609](https://arxiv.org/abs/2301.08609)). The [AQC Dynamics Template](/docs/guides/function-template-aqc-trotter) packages this whole quantum core (Trotter synthesis, AQC compression, and mitigated execution) behind one call:\n", "\n", "| PRE (this notebook) | FUNCTION (`aqc-dynamics-function`) | POST (this notebook) |\n", "| -------------------------------------------------------------------------------------------------------- | ------------------------------------------------------------------------------------------------------------------------------------------ | ---------------------------------------------- |\n", "| Ground state from DMRG plus MPS-fidelity maximization, with the neutron kick baked into the same circuit | Trotter synthesis → AQC compression → execution on `statevector`, `fake`, or `runtime`, returning $\\langle \\sigma_z^j \\rangle(t)$ per site | $S(q, \\omega)$, the dynamical structure factor |\n", "\n", "The experiment-specific work stays here in the notebook: ground-state preparation (PRE) and the $S(q, \\omega)$ post-processing (POST). The two quantum-heavy steps, compression and execution, run inside the function.\n", "\n", "This tutorial is a companion to [Simulate neutron scattering in quantum materials with quantum circuits](/docs/tutorials/simulate-neutron-scattering), which builds the same experiment inline: the same KCuF$_3$ model, ground-state preparation, neutron kick, and post-processing, with the Trotter synthesis, AQC compression, and mitigated execution written out step by step. Read that tutorial to learn how AQC compression works. Read this one to run the same experiment through a deployed function template: the quantum core becomes a single function call, and the multi-hour AQC compression runs inside the Serverless worker instead of on your machine, so you do not need an HPC system or an open kernel while it runs. The same call also drives other 1D dynamics experiments.\n", "\n" ] }, { "cell_type": "markdown", "id": "requirements", "metadata": {}, "source": [ "## Requirements\n", "\n", "Before starting this tutorial, be sure you have the following:\n", "\n", "* The function deployed to your Qiskit Serverless account. Run the companion function template first: [Deploy and run the AQC + Trotter dynamics function template](/docs/guides/function-template-aqc-trotter). That guide walks through getting the source files and uploading the function to your account. This tutorial only calls the deployed function.\n", "\n", "* IBM Quantum® credentials saved for `QiskitServerless` (see the function template). Both examples in this tutorial call the deployed function, so both need them.\n", "\n", "* Qiskit SDK v2.0 or later (`pip install qiskit`).\n", "\n", "* The Qiskit IBM Catalog client (`pip install qiskit-ibm-catalog`).\n", "\n", "* NumPy, SciPy, and Matplotlib (`pip install numpy scipy matplotlib`). SciPy 1.14 or later is needed for the COBYQA optimizer used in ground-state preparation.\n", "\n", "* The AQC tensor-network stack, because the ground-state preparation in Step 1 runs locally in this notebook: `pip install 'qiskit-addon-aqc-tensor[quimb-jax]==0.3.1'`.\n", "\n", "The first call to a newly deployed function waits while the Serverless worker installs its dependencies, so expect extra latency on that run.\n", "\n" ] }, { "cell_type": "markdown", "id": "setup-md", "metadata": {}, "source": [ "## Setup\n", "\n", "Import the libraries and define the experiment-specific helpers used later: `build_gs_ansatz` (the Hamiltonian variational ansatz, or HVA, for ground-state preparation), `prepare_ground_state` (DMRG plus MPS-fidelity maximization), and `get_spectrum`, `plot_green`, and `plot_spectrum` (the $S(q, \\omega)$ post-processing). These are adapted from the [original neutron-scattering tutorial](/docs/tutorials/simulate-neutron-scattering).\n", "\n" ] }, { "cell_type": "code", "execution_count": 1, "id": "setup-imports", "metadata": {}, "outputs": [], "source": [ "from functools import partial\n", "\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import scipy.optimize\n", "\n", "import quimb.tensor as qtn\n", "from qiskit import QuantumCircuit\n", "from qiskit.quantum_info import SparsePauliOp\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_catalog import QiskitServerless" ] }, { "cell_type": "code", "execution_count": 2, "id": "setup-helpers", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Setup complete - helpers defined.\n" ] } ], "source": [ "# Dynamical structure factor via discrete Fourier transform\n", "\n", "\n", "def get_spectrum(n, Gjjc, dt, time_steps, q_steps, w_steps):\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", " \"\"\"\n", " green = Gjjc / 4 # sigma -> S=1/2\n", " omega_max = np.pi / dt\n", " qpoints = np.arange(0, 2 * np.pi, 2 * np.pi / q_steps)\n", " omegas = np.arange(0, omega_max, omega_max / w_steps)\n", " green_map = np.zeros((omegas.shape[0], qpoints.shape[0]))\n", " center = n // 2 - 1\n", " for iw, w in enumerate(omegas):\n", " exponent = np.exp(1j * w * dt * np.arange(1, time_steps + 1))\n", " S_w = np.dot(green.T, exponent) * dt\n", " for iq, q in enumerate(qpoints):\n", " q_matrix = np.exp(-1j * q * np.arange(-center, center + 2, 1))\n", " green_map[iw, iq] = np.imag(np.dot(S_w, q_matrix))\n", " return green_map\n", "\n", "\n", "# Plotting helpers\n", "\n", "\n", "def plot_spectrum(\n", " dsf,\n", " dt,\n", " q_steps,\n", " w_steps,\n", " lower_bound=False,\n", " upper_bound=False,\n", " title=None,\n", "):\n", " \"\"\"Heat-map of the dynamical structure factor.\"\"\"\n", " omega_max = np.pi / dt\n", " qpoints = np.arange(0, 2 * np.pi, 2 * np.pi / q_steps)\n", " omegas = np.arange(0, omega_max, omega_max / w_steps)\n", " x, y = np.meshgrid(qpoints, omegas)\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", " if lower_bound:\n", " ax.plot(\n", " qpoints,\n", " np.pi * np.abs(np.sin(qpoints)) / 2,\n", " \"--\",\n", " color=\"white\",\n", " lw=1.5,\n", " label=\"Lower bound\",\n", " )\n", " if upper_bound:\n", " ax.plot(\n", " qpoints,\n", " np.pi * np.abs(np.sin(qpoints / 2)),\n", " \"--\",\n", " color=\"red\",\n", " lw=1.5,\n", " label=\"Upper bound\",\n", " )\n", " ax.set_ylim(0, 3.6)\n", " ax.set_xlim(0, 2 * np.pi - 2 * np.pi / q_steps)\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 lower_bound or upper_bound:\n", " ax.legend(loc=\"upper right\", fontsize=11)\n", " if title:\n", " ax.set_title(title, fontsize=14)\n", " plt.tight_layout()\n", " plt.show()\n", "\n", "\n", "def plot_green(n, Gjjc, time_steps, dt, title=None):\n", " \"\"\"Heat-map of the retarded Green's function in real space and time.\"\"\"\n", " fig, ax = plt.subplots(figsize=(8, 6))\n", " t_axis = np.arange(1, time_steps + 1) * dt\n", " site_axis = np.arange(n)\n", " x, y = np.meshgrid(t_axis, site_axis)\n", " c = ax.pcolormesh(\n", " x,\n", " y,\n", " np.real(Gjjc).T,\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(r\"Time ($t / J^{-1}$)\", fontsize=16)\n", " ax.set_ylabel(\"Site index $j$\", fontsize=16)\n", " if title:\n", " ax.set_title(title, fontsize=14)\n", " plt.tight_layout()\n", " plt.show()\n", "\n", "\n", "# Variational ground-state ansatz (HVA)\n", "\n", "\n", "def _apply_xxz_pair_gate(qc, q0, q1, theta):\n", " \"\"\"Apply the parameterized XXZ-type two-qubit gate used in the HVA.\"\"\"\n", " qc.cx(q0, q1)\n", " qc.rz(theta, q1)\n", " qc.h(q0)\n", " qc.rz(theta + np.pi / 2, q0)\n", " qc.cx(q0, q1)\n", " qc.rz(-theta, q1)\n", " qc.h(q1)\n", " qc.cx(q1, q0)\n", " qc.rz(np.pi / 2, q1)\n", " qc.rz(-np.pi / 2, q0)\n", " qc.h(q1)\n", " qc.h(q0)\n", "\n", "\n", "def build_gs_ansatz(n, params, layers):\n", " \"\"\"Build the Hamiltonian variational ansatz (HVA) circuit for\n", " ground-state preparation of the 1D Heisenberg model.\n", "\n", " Starts from a product of singlet pairs and applies alternating\n", " odd/even layers of parameterized XXZ gates. For layer r,\n", " params[2 * r] is the odd-layer (inter-pair) angle and\n", " params[2 * r + 1] is the even-layer (intra-pair) angle.\n", " \"\"\"\n", " qc = QuantumCircuit(n)\n", " # Initial singlet product state\n", " for i in range(n // 2):\n", " qc.x(2 * i)\n", " qc.x(2 * i + 1)\n", " qc.h(2 * i + 1)\n", " qc.cx(2 * i + 1, 2 * i)\n", " # Variational layers\n", " for r in range(layers):\n", " for i in range(1, (n + 1) // 2): # odd layer\n", " _apply_xxz_pair_gate(qc, 2 * i - 1, 2 * i, params[2 * r])\n", " for i in range(n // 2): # even layer\n", " _apply_xxz_pair_gate(qc, 2 * i, 2 * i + 1, params[2 * r + 1])\n", " return qc\n", "\n", "\n", "def prepare_ground_state(n, gs_layers=5, max_bond=128, cutoff=1e-8):\n", " \"\"\"Prepare the KCuF3 (isotropic Heisenberg) ground state as a QuantumCircuit.\n", "\n", " Runs DMRG (quimb MPO + DMRG2) to get the chain's ground state, then optimizes\n", " the HVA angles to maximize the MPS overlap ||^2. No exact\n", " diagonalization, so it scales to larger n.\n", " \"\"\"\n", " J = Jz = 1.0\n", " builder = qtn.SpinHam1D(S=1 / 2)\n", " builder += J * 0.5, \"+\", \"-\"\n", " builder += J * 0.5, \"-\", \"+\"\n", " builder += Jz, \"Z\", \"Z\"\n", " H_mpo = builder.build_mpo(L=n)\n", " dmrg = qtn.DMRG2(H_mpo)\n", " dmrg.solve(tol=1e-8, verbosity=0)\n", "\n", " gs_sim = QuimbSimulator(\n", " quimb_circuit_factory=partial(\n", " qtn.CircuitMPS, gate_opts=dict(cutoff=cutoff, max_bond=max_bond)\n", " ),\n", " autodiff_backend=\"jax\",\n", " )\n", "\n", " def gs_infidelity(params):\n", " psi = tensornetwork_from_circuit(\n", " build_gs_ansatz(n, params, gs_layers), gs_sim\n", " ).psi\n", " return 1 - abs(psi.H @ dmrg.state) ** 2\n", "\n", " # Seed and optimizer match the original tutorial. Each layer starts at\n", " # [0, pi/2]: an odd-layer angle of 0 makes the inter-pair gate the identity,\n", " # and an even-layer angle of pi/2 makes the intra-pair gate a SWAP (since\n", " # 0.5 * (XX + YY + ZZ) = SWAP - I/2). That puts the seed at the singlet-pair\n", " # product limit, which is already a decent approximation to the Heisenberg\n", " # ground state, so the optimizer only has to refine it. The small jitter\n", " # (fixed RNG seed, so runs are reproducible) breaks the exact symmetry\n", " # between layers; COBYQA then runs for up to 100 iterations.\n", " rng = np.random.default_rng(12345)\n", " x0 = np.tile([0.0, np.pi / 2], gs_layers) + rng.normal(\n", " scale=0.1, size=2 * gs_layers\n", " )\n", " result_gs = scipy.optimize.minimize(\n", " gs_infidelity, x0, method=\"COBYQA\", options={\"maxiter\": 100}\n", " )\n", " print(f\"DMRG ground-state energy: {dmrg.energy:.6f}\")\n", " print(f\"GS fidelity: {1 - result_gs.fun:.4f}\")\n", " return build_gs_ansatz(n, result_gs.x, gs_layers)\n", "\n", "\n", "print(\"Setup complete - helpers defined.\")" ] }, { "cell_type": "markdown", "id": "setup-load-md", "metadata": {}, "source": [ "### Load the function template\n", "\n", "Connect to Qiskit Serverless and load the deployed `aqc-dynamics-function`. Both examples in this tutorial call the same `fn` handle, so the function is loaded once, here.\n", "\n" ] }, { "cell_type": "code", "execution_count": 3, "id": "setup-connect", "metadata": {}, "outputs": [], "source": [ "# Credentials are read from the account saved once via QiskitServerless.save_account(...)\n", "serverless = QiskitServerless()\n", "fn = serverless.load(\"aqc-dynamics-function\")" ] }, { "cell_type": "markdown", "id": "small-md", "metadata": {}, "source": [ "## Small-scale simulator example\n", "\n", "We first run the full workflow on a small 10-site chain by using the exact `statevector` backend. This validates the PRE → FUNCTION → POST pipeline before spending any QPU time.\n", "\n" ] }, { "cell_type": "markdown", "id": "small-s1-md", "metadata": {}, "source": [ "### Step 1: Map classical inputs to a quantum problem\n", "\n", "Build the KCuF$_3$ Hamiltonian as a `SparsePauliOp` (isotropic Heisenberg: $XX + YY + ZZ$ at coupling $\\tfrac14$ on each nearest-neighbor bond; the strings are Pauli operators, so $\\tfrac14$ gives the spin-$\\frac{1}{2}$ coupling). Prepare the ground state with DMRG plus MPS-fidelity maximization, then bake in the neutron kick: a $\\pi/2$ $Z$-rotation at the center site. The prepared circuit is what we hand to the function as `initial_state`. We leave `observables` at its default (per-site $Z$), which is exactly the $\\langle \\sigma_z^j \\rangle(t)$ readout the neutron workflow needs.\n", "\n" ] }, { "cell_type": "code", "execution_count": 4, "id": "small-s1-code", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "DMRG ground-state energy: -4.258035\n", "GS fidelity: 0.9841\n", "Prepared 10-qubit ground state with the neutron kick at site 4.\n" ] } ], "source": [ "n = 10\n", "dt = 0.6 # physical time per Trotter step (also the omega-axis unit in POST)\n", "time_steps = 10\n", "center = n // 2 - 1\n", "\n", "# MPS-simulator settings, shared by the ground-state prep here and the AQC\n", "# compression inside the function (matches the original tutorial).\n", "mps_max_bond = 32\n", "mps_cutoff = 1e-8\n", "\n", "# 1D isotropic Heisenberg (KCuF3) Hamiltonian on n qubits\n", "H = SparsePauliOp.from_sparse_list(\n", " [(p, [i, i + 1], 0.25) for i in range(n - 1) for p in (\"XX\", \"YY\", \"ZZ\")],\n", " num_qubits=n,\n", ")\n", "\n", "# Ground state (DMRG + fidelity max) + neutron kick baked into the same circuit\n", "gs_circuit = prepare_ground_state(\n", " n, gs_layers=3, max_bond=mps_max_bond, cutoff=mps_cutoff\n", ")\n", "gs_circuit.rz(\n", " np.pi / 2, center\n", ") # exp(-i (pi/2)/2 Z_center): the neutron perturbation\n", "print(\n", " f\"Prepared {n}-qubit ground state with the neutron kick at site {center}.\"\n", ")" ] }, { "cell_type": "markdown", "id": "small-s3-md", "metadata": {}, "source": [ "### Steps 2 and 3: Compress and execute with the function template\n", "\n", "In a hand-written workflow these are two separate stages: optimize the circuits for hardware (Step 2) and execute them (Step 3). The function template collapses both into one call. It performs Trotter synthesis, AQC compression, and hardware transpilation, then runs the circuits (here on the exact simulator, later with built-in error mitigation on hardware). The two tuning parameters are `aqc_segments` (the compression plan) and `aqc_options` (the MPS and optimizer settings). Each segment `{\"n_steps\": k, \"ansatz_steps\": m}` compresses `k` consecutive Trotter steps into an ansatz built from an `m`-step Trotter target, and any steps beyond `sum(n_steps)` run as plain Trotter. Early, low-entanglement steps compress well into a shallow (`ansatz_steps=1`) ansatz, so here we compress the first three steps into a single-layer ansatz and the next two into a deeper two-layer ansatz; the remaining five of the 10 Trotter steps run as plain Trotter. For `aqc_options` we mirror the original tutorial: MPS bond dimension `max_bond=32`, `cutoff=1e-8`, and an L-BFGS-B optimizer capped at 100 iterations.\n", "\n", "Call the function loaded in Setup. `backend=\"statevector\"` runs the exact reference path: no QPU time, with the circuits running on an exact statevector simulator inside the serverless worker (a saved Qiskit Serverless account is still needed to call it). The `initial_state` carries the prepared ground state (including the kick); `observables` is omitted so the function measures the default per-site $Z$.\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "small-s3-run", "metadata": {}, "outputs": [], "source": [ "job = fn.run(\n", " t_steps=time_steps,\n", " aqc_segments=[\n", " {\n", " \"n_steps\": 3,\n", " \"ansatz_steps\": 1,\n", " }, # early steps -> shallow 1-layer ansatz\n", " {\n", " \"n_steps\": 2,\n", " \"ansatz_steps\": 2,\n", " }, # later steps -> deeper 2-layer ansatz\n", " ],\n", " aqc_options={\n", " \"max_bond\": mps_max_bond, # MPS bond dimension for AQC compression\n", " \"cutoff\": mps_cutoff,\n", " \"optimizer_settings\": {\n", " \"method\": \"L-BFGS-B\",\n", " \"jac\": True,\n", " \"options\": {\"maxiter\": 100},\n", " },\n", " },\n", " dt=dt,\n", " hamiltonian=H,\n", " initial_state=gs_circuit, # prepared ground state including the neutron kick\n", " # observables omitted -> default per-site Z (the neutron sigma_z readout)\n", " backend=\"statevector\",\n", ")" ] }, { "cell_type": "code", "execution_count": 20, "id": "small-status", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "DONE\n" ] } ], "source": [ "print(job.status()) # rerun this cell until status says DONE" ] }, { "cell_type": "code", "execution_count": 21, "id": "small-s3-result", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "AQC fidelities: {'1': 1.0, '2': 0.9999, '3': 0.9992, '4': 0.9998, '5': 0.9995}\n", "Green's function shape: (10, 10)\n" ] } ], "source": [ "# The per-site (t) the function returns is the retarded Green's function\n", "# G(j, j_c, t). The workflow samples t = 1..time_steps, so drop the t = 0 row (the\n", "# prepared+kicked state before any evolution) before post-processing.\n", "result = job.result()\n", "print(\n", " \"AQC fidelities:\",\n", " {k: round(v, 4) for k, v in result[\"metadata\"][\"aqc_fidelities\"].items()},\n", ")\n", "\n", "ev = np.array(result[\"expectation_values\"])\n", "Gjjc = ev[1:] # shape (time_steps, n)\n", "print(\"Green's function shape:\", Gjjc.shape)" ] }, { "cell_type": "markdown", "id": "small-s4-md", "metadata": {}, "source": [ "### Step 4: Post-process and return result in desired classical format\n", "\n", "Fourier-transform the Green's function into $S(q, \\omega)$, mirror-symmetrize, and clip negatives: the standard neutron post-processing. Mirroring is exact because $S(q, \\omega) = S(-q, \\omega)$ for this model, and the negative values that survive are artifacts of Fourier-transforming a finite, discretely sampled time series, so they are clipped to zero. On this small exact run the two-spinon continuum is only coarsely resolved, but the machinery is identical to the hardware run that follows.\n", "\n" ] }, { "cell_type": "code", "execution_count": 22, "id": "small-s4-code", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "q_res, w_res = 100, 100\n", "spectrum = get_spectrum(n, Gjjc, dt, time_steps, q_res, w_res)\n", "spectrum = -(spectrum + spectrum[:, ::-1]) / 2 # mirror symmetry\n", "spectrum = np.clip(spectrum, a_min=0, a_max=None) # clip negatives\n", "\n", "plot_green(\n", " n,\n", " Gjjc,\n", " time_steps,\n", " dt,\n", " title=f\"Retarded Green's function - {n} qubits (AQC, statevector)\",\n", ")\n", "plot_spectrum(\n", " spectrum,\n", " dt,\n", " q_res,\n", " w_res,\n", " lower_bound=True,\n", " upper_bound=True,\n", " title=f\"Dynamical structure factor - {n} qubits (AQC, statevector)\",\n", ")" ] }, { "cell_type": "markdown", "id": "large-md", "metadata": {}, "source": [ "## Large-scale hardware example\n", "\n", "The same workflow scales up without changing any of the science code: a 30-site chain, twice the Trotter depth (20 steps), a compression plan that varies the ansatz depth (a deeper ansatz for the later, more-entangled steps), and execution on an IBM Quantum processor with the function's built-in error mitigation (dynamical decoupling, Pauli twirling, and twirled readout error extinction (TREX)). We walk through the same four steps as the simulator example, reusing the `fn` handle from Setup.\n", "\n", "| | Small scale | Large scale |\n", "| ---------------------------------------- | ------------- | ------------------------------------- |\n", "| Qubits | 10 | 30 |\n", "| Trotter steps | 10 | 20 |\n", "| AQC-compressed steps (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", "| Backend | `statevector` | QPU with DD, Pauli twirling, and TREX |\n", "\n" ] }, { "cell_type": "markdown", "id": "large-s1-md", "metadata": {}, "source": [ "### Step 1: Map classical inputs to a quantum problem\n", "\n", "Build the same KCuF$_3$ Heisenberg `SparsePauliOp` and prepare the ground state, now with a deeper `gs_layers=5` ansatz for the longer chain, then bake in the $\\pi/2$ $Z$ neutron kick at the center site. This is identical to the small-scale mapping, but at $n = 30$.\n", "\n", "Expect a lower ground-state fidelity than the 10-site run: around 0.82 here against 0.98 for the smaller chain, because five HVA layers cannot fully capture a 30-site ground state. That is expected rather than a failure, and the original tutorial accepts roughly 0.65 at 50 sites for the same reason. Raising `gs_layers` or the COBYQA iteration cap improves it, at extra classical cost.\n", "\n" ] }, { "cell_type": "code", "execution_count": 34, "id": "large-code", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "DMRG ground-state energy: -13.111355\n", "GS fidelity: 0.8201\n", "Prepared 30-qubit ground state with the neutron kick at site 14.\n" ] } ], "source": [ "n = 30\n", "dt = 0.6\n", "time_steps = 20\n", "center = n // 2 - 1\n", "\n", "# Same MPS settings as the original large-scale run: a larger bond for the\n", "# longer, more-entangled chain (shared by GS prep and AQC compression).\n", "mps_max_bond = 128\n", "mps_cutoff = 1e-8\n", "\n", "# Same KCuF3 Hamiltonian and ground-state prep, on a larger chain\n", "H = SparsePauliOp.from_sparse_list(\n", " [(p, [i, i + 1], 0.25) for i in range(n - 1) for p in (\"XX\", \"YY\", \"ZZ\")],\n", " num_qubits=n,\n", ")\n", "gs_circuit = prepare_ground_state(\n", " n, gs_layers=5, max_bond=mps_max_bond, cutoff=mps_cutoff\n", ")\n", "gs_circuit.rz(np.pi / 2, center) # neutron kick at the center site\n", "print(\n", " f\"Prepared {n}-qubit ground state with the neutron kick at site {center}.\"\n", ")" ] }, { "cell_type": "markdown", "id": "large-s3-md", "metadata": {}, "source": [ "### Steps 2 and 3: Compress and execute with the function template\n", "\n", "The same single call as the simulator example, now with `backend_name` pointing at an IBM Quantum processor, so the function transpiles and executes there. The compression plan varies the ansatz depth: the first six (low-entanglement) Trotter steps compress into a shallow single-layer ansatz, the next four into a deeper two-layer ansatz, and the remaining 10 of the 20 steps run as plain Trotter. `aqc_options` raises the MPS bond dimension to `max_bond=128` for the longer, more-entangled chain (matching the original), keeping the same L-BFGS-B optimizer capped at 100 iterations. The `estimator_options` turn on the built-in error mitigation: dynamical decoupling (XY4), gate twirling, and TREX measurement mitigation. The function's defaults already match the original tutorial for all of these except the TREX learning budget (`measure_noise_learning`). The whole block is still written out because a caller-supplied `estimator_options` replaces the function's defaults wholesale instead of merging into them, so omitting a key would fall back to the IBM Quantum Compute default rather than the function's.\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "large-s3-run", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "job ID (save this to reconnect later): 43ed8d07-6d7d-4f33-b70a-7f31b765b310\n" ] } ], "source": [ "# Steps 2 + 3: the function compresses (varied ansatz) and executes on hardware.\n", "job = fn.run(\n", " t_steps=time_steps,\n", " aqc_segments=[\n", " {\n", " \"n_steps\": 6,\n", " \"ansatz_steps\": 1,\n", " }, # early steps -> shallow 1-layer ansatz\n", " {\n", " \"n_steps\": 4,\n", " \"ansatz_steps\": 2,\n", " }, # later steps -> deeper 2-layer ansatz\n", " ],\n", " aqc_options={\n", " \"max_bond\": mps_max_bond, # 128 for the longer chain\n", " \"cutoff\": mps_cutoff,\n", " \"optimizer_settings\": {\n", " \"method\": \"L-BFGS-B\",\n", " \"jac\": True,\n", " \"options\": {\"maxiter\": 100},\n", " },\n", " },\n", " dt=dt,\n", " hamiltonian=H,\n", " initial_state=gs_circuit,\n", " backend_name=\"ibm_pittsburgh\",\n", " # Mitigation settings from the original tutorial. Only the two\n", " # measure_noise_learning values differ from the function's defaults; the rest\n", " # restates them, because a caller-supplied estimator_options dict replaces the\n", " # function's defaults wholesale rather than merging into them.\n", " estimator_options={\n", " \"environment\": {\"job_tags\": [\"TUT-SNS\"]},\n", " \"dynamical_decoupling\": {\"enable\": True, \"sequence_type\": \"XY4\"},\n", " \"twirling\": {\n", " \"enable_gates\": True,\n", " \"num_randomizations\": 1000,\n", " \"shots_per_randomization\": 128,\n", " },\n", " \"resilience\": {\n", " \"measure_mitigation\": True,\n", " \"measure_noise_learning\": {\n", " \"num_randomizations\": 32,\n", " \"shots_per_randomization\": 100,\n", " },\n", " },\n", " },\n", ")\n", "print(\"job ID (save this to reconnect later):\", job.job_id)" ] }, { "cell_type": "markdown", "id": "large-reconnect-md", "metadata": {}, "source": [ "\n", " The large-scale run is not quick, and most of the time is classical rather than on the QPU. The AQC compression runs inside the function before anything reaches the QPU: at 30 sites with `max_bond=128` that took close to four hours in our run, against the roughly 18 minutes of QPU time quoted in the *Usage estimate* at the top of this tutorial. Queue wait is on top of both. You do not need to keep this notebook or kernel open while it runs.\n", "\n", " Copy the job ID printed by the preceding cell and save it. The next three cells let you pick the run back up later:\n", "\n", " 1. Reconnect, only needed in a new kernel session: re-run the [Setup](#setup) cells to recreate `serverless`, then rebuild the `job` handle from the ID you saved. Skip this cell if you are still in the session where you submitted, because the handle is already live.\n", " 2. Check status: re-run until it reports `DONE`.\n", " 3. Fetch the result: run only once the status is `DONE`.\n", "\n", " The following reconnect cell holds a placeholder. Replace it with your own `job_id`:\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "large-reconnect", "metadata": {}, "outputs": [], "source": [ "# Reconnect to a previously submitted job by its ID. Only needed in a NEW kernel\n", "# session; if you are still in the session where you submitted, the `job` handle\n", "# from the preceding cell is already live, so skip this cell. Replace the ID that follows with your own.\n", "job = serverless.get_job_by_id(\"\")" ] }, { "cell_type": "code", "execution_count": 27, "id": "large-status", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "DONE\n" ] } ], "source": [ "# Check where the job is. Re-run this until it reports DONE before fetching the\n", "# result in the following cell: QUEUED -> INITIALIZING -> RUNNING: OPTIMIZING_FOR_HARDWARE ->\n", "# RUNNING: WAITING_FOR_QPU -> RUNNING: EXECUTING_QPU -> RUNNING: POST_PROCESSING\n", "# -> DONE.\n", "print(job.status())" ] }, { "cell_type": "code", "execution_count": null, "id": "large-s3-result", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "AQC fidelities: {'1': 1.0, '2': 0.9994, '3': 0.9944, '4': 0.9853, '5': 0.9747, '6': 0.959, '7': 0.9495, '8': 0.9542, '9': 0.9533, '10': 0.9451}\n" ] } ], "source": [ "# Run this only once the preceding status cell reports DONE. result() blocks until\n", "# the job finishes, so calling it earlier just waits (possibly for hours).\n", "result = job.result()\n", "print(\n", " \"AQC fidelities:\",\n", " {k: round(v, 4) for k, v in result[\"metadata\"][\"aqc_fidelities\"].items()},\n", ")\n", "\n", "ev = np.array(result[\"expectation_values\"])\n", "Gjjc = ev[1:] # drop the t = 0 row -> shape (time_steps, n)" ] }, { "cell_type": "markdown", "id": "large-s4-md", "metadata": {}, "source": [ "### Step 4: Post-process and return result in desired classical format\n", "\n", "Identical post-processing to the simulator run: Fourier-transform the Green's function into $S(q, \\omega)$, mirror-symmetrize, and clip negatives. With the longer chain and evolution, the two-spinon continuum is far better resolved. It should fill the band between the dashed bounds, brightest near $q = \\pi$.\n", "\n" ] }, { "cell_type": "code", "execution_count": 33, "id": "large-result", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "n = result[\"metadata\"][\"n\"]\n", "q_res, w_res = 100, 100\n", "spectrum = get_spectrum(n, Gjjc, dt, time_steps, q_res, w_res)\n", "spectrum = -(spectrum + spectrum[:, ::-1]) / 2 # mirror symmetry\n", "spectrum = np.clip(spectrum, a_min=0, a_max=None) # clip negatives\n", "\n", "plot_green(\n", " n,\n", " Gjjc,\n", " time_steps,\n", " dt,\n", " title=f\"Retarded Green's function - {n} qubits (AQC, hardware)\",\n", ")\n", "plot_spectrum(\n", " spectrum,\n", " dt,\n", " q_res,\n", " w_res,\n", " lower_bound=True,\n", " upper_bound=True,\n", " title=f\"Dynamical structure factor - {n} qubits (AQC, hardware)\",\n", ")" ] }, { "cell_type": "markdown", "id": "appendix-md", "metadata": {}, "source": [ "## Appendix\n", "\n", "The preceding hardware example runs a single chain length. The three spectra that follow come from earlier hardware runs of this same workflow on `ibm_pittsburgh` at 10, 20, and 30 sites, with every other input held fixed: 20 Trotter steps at `dt = 0.6`, the compression plan of six one-layer plus four two-layer AQC-compressed steps, and `max_bond = 128`. These are recorded results, not output from the preceding cells.\n", "\n", "The same settings are used at all three sizes, so the spectra are directly comparable. Tuning them per chain length, with more ground-state ansatz layers or a larger `max_bond`, for example, can give better results than any shown here.\n", "\n", "![Dynamical structure factor at 10 sites, a single sharp bright peak at q = pi near the lower bound](https://eu-de.quantum.cloud.ibm.com/docs/images/tutorials/simulate-neutron-scattering-with-a-serverless-workflow/appendix-dsf-10.avif \"10 qubits\")\n", "\n", "![Dynamical structure factor at 20 sites, spectral weight filling the band between the two dashed two-spinon bounds](https://eu-de.quantum.cloud.ibm.com/docs/images/tutorials/simulate-neutron-scattering-with-a-serverless-workflow/appendix-dsf-20.avif \"20 qubits\")\n", "\n", "![Dynamical structure factor at 30 sites, the continuum resolved more finely with fainter contrast and some weight outside the bounds](https://eu-de.quantum.cloud.ibm.com/docs/images/tutorials/simulate-neutron-scattering-with-a-serverless-workflow/appendix-dsf-30.avif \"30 qubits\")\n", "\n" ] }, { "cell_type": "markdown", "id": "nextsteps", "metadata": {}, "source": [ "## Next steps\n", "\n", "\n", " * Adapt this workflow to your own system: the function accepts any 1D nearest-neighbor `SparsePauliOp`, so a different chain Hamiltonian, initial state, or set of observables runs the same PRE → FUNCTION → POST pipeline. See the full input/output contract in the [AQC Dynamics Template on GitHub](https://github.com/qiskit-community/qiskit-function-templates/tree/main/physics/aqc_trotter).\n", " * Read the paper this benchmark comes from: Lee et al., [*Benchmarking quantum simulation with neutron-scattering experiments*](https://arxiv.org/abs/2603.15608) (arXiv:2603.15608).\n", " * Compare with the [original \"Simulate neutron scattering\" tutorial](/docs/tutorials/simulate-neutron-scattering), the inline workflow this one ports onto a deployed function template.\n", " * Go deeper on the [error mitigation and suppression techniques](/docs/guides/error-mitigation-and-suppression-techniques) applied on the hardware run: dynamical decoupling, Pauli twirling, and TREX.\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "id": "a1b8767d", "source": "© IBM Corp., 2017-2026" } ], "metadata": { "hours": 4, "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" }, "qpuSeconds": 1080 }, "nbformat": 4, "nbformat_minor": 5 }