{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "aedb05f4-d7f8-4cb1-9aab-8e7081946c79",
   "metadata": {},
   "source": [
    "# Entanglement Accelerates Quantum Simulation \n",
    "\n",
    "This notebook demonstrates how to reproduce the key numerical result from [*Entanglement Accelerates Quantum Simulation*](https://www.nature.com/articles/s41567-025-02945-2) using NVIDIA CUDA-Q. We follow the workflow presented in the paper for the quantum Ising model with mixed fields (QIMF) and show how entanglement growth can *reduce* the observed Trotter error of the first-order product formula (PF1).\n",
    "\n",
    "Notebook prepared by Jonah So (jonahso@ucsb.edu) and Jue Xu (xujue@connect.hku.hk)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "7dfdff13-d146-48fa-b997-1cc346e09487",
   "metadata": {
    "id": "7dfdff13-d146-48fa-b997-1cc346e09487"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CUDA-Q Version  (https://github.com/NVIDIA/cuda-quantum e50d28149779ca5a9c6757c8a93c670aa1febfc6)\n"
     ]
    }
   ],
   "source": [
    "# Relevant imports\n",
    "\n",
    "import cudaq\n",
    "import numpy as np\n",
    "import cupy as cp\n",
    "from cupyx.scipy.linalg import expm\n",
    "import time, sys\n",
    "from cudaq import spin\n",
    "from cudaq import State\n",
    "import matplotlib.pyplot as plt\n",
    "from typing import List\n",
    "\n",
    "print(cudaq.__version__)\n",
    "cudaq.set_target(\"nvidia\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e7a3747",
   "metadata": {},
   "source": [
    "## 2. Model Definition\n",
    "\n",
    "We focus on the $n$-qubit quantum Ising model with mixed fields (QIMF), the same Hamiltonian used throughout the paper. The subsections below walk through the exact construction, the PF1 Trotter step, and the helper routines needed for later analysis.\n",
    "\n",
    "### 2.1 Initial product state\n",
    "We start from the $|0000\\ldots\\rangle$ state, matching the initial condition explored in the paper. The helper kernel below prepares this state.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "184bb2a6-9e34-4aa4-ae59-4fafaac5a5e4",
   "metadata": {
    "id": "184bb2a6-9e34-4aa4-ae59-4fafaac5a5e4"
   },
   "outputs": [],
   "source": [
    "@cudaq.kernel\n",
    "def get_initial_state(n_spins: int):\n",
    "    '''\n",
    "    Create initial state |0000...> for n_spins\n",
    "    ''' \n",
    "    qubits = cudaq.qvector(n_spins)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "233965a0-62fa-4dbe-8ddc-c04287946eb3",
   "metadata": {
    "id": "233965a0-62fa-4dbe-8ddc-c04287946eb3"
   },
   "source": [
    "### 2.2 QIMF Hamiltonian\n",
    "\n",
    "The paper studies rapid entanglement growth in the one-dimensional QIMF Hamiltonian\n",
    "\n",
    "$$\n",
    "H = h_x \\sum_{j=1}^N X_j + h_y \\sum_{j=1}^N Y_j + J \\sum_{j=1}^{N-1} X_j X_{j+1},\n",
    "$$\n",
    "\n",
    "where $h_x$ and $h_y$ set the mixed fields and $J$ couples nearest neighbours. Every term is expressed in the Pauli basis, so we can assemble $H$ with CUDA-Q’s `spin` helpers. For later Trotterization we split the Hamiltonian into two commuting groups where the terms in the same group commute with each other:\n",
    "\n",
    "- $A = h_x \\sum_j X_j + J \\sum_j X_j X_{j+1}$\n",
    "- $B = h_y \\sum_j Y_j$\n",
    "\n",
    "The function below returns $(A, B, H)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "b0fed87f-0ea8-4ea3-8d65-a98252f00553",
   "metadata": {
    "id": "b0fed87f-0ea8-4ea3-8d65-a98252f00553"
   },
   "outputs": [],
   "source": [
    "def create_hamiltonian_qimf(n_spins: int, J: float, h_x: float, h_y: float):\n",
    "\n",
    "    A = 0\n",
    "    B = 0\n",
    "    \n",
    "    for j in range(n_spins - 1):\n",
    "        A += J * spin.x(j) * spin.x(j + 1)\n",
    "\n",
    "    for j in range(n_spins):\n",
    "        A += h_x * spin.x(j)\n",
    "        B += h_y * spin.y(j)\n",
    "        \n",
    "    H = A + B\n",
    "\n",
    "    return A, B, H"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cad3eef1",
   "metadata": {},
   "source": [
    "### 2.3 First-Order Trotter Formula (PF1)\n",
    "\n",
    "The paper benchmarks the first-order product formula as a baseline algorithm. For a Hamiltonian decomposition $H = \\sum_l H_l$, the ideal evolution $U(\\delta t) = e^{-i H \\delta t}$ is approximated by the ordered product of exponentials\n",
    "$$\n",
    "\\mathscr{U}_1(\\delta t) = \\overrightarrow{\\prod_{l}} e^{-i H_l \\, \\delta t},\n",
    "$$\n",
    "which incurs an error of $\\mathcal{O}(\\delta t^2)$ whenever the terms fail to commute. \n",
    "The authors show that once the simulated state has significant entanglement, the empirical PF1 error drops sharply to the average-case performance.\n",
    "\n",
    "The CUDA-Q kernel below implements a single PF1 step for two commuting groups (the $A$ and $B$ partitions defined above). "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e73881a-d0f0-4818-9ded-684320ce6e66",
   "metadata": {
    "id": "4e73881a-d0f0-4818-9ded-684320ce6e66"
   },
   "source": [
    "### 2.4 PF1 step for the QIMF partition\n",
    "\n",
    "Because each $A$ term commutes with every other $A$ term (and similarly for $B$), the PF1 update factorises neatly into\n",
    "$$\n",
    "\\mathscr{U}_1(\\delta t) = e^{-i A \\delta t} \\, e^{-i B \\delta t}.\n",
    "$$\n",
    "To simulate for a total time $t$, we apply this kernel $r$ times with $\\delta t = t/r$. The paper sets $t=10$ and $r=100$, which we adopt below so the numerical values line up with Figure 1a."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "89e43405-4d97-4c6c-8396-3e80ba3b7659",
   "metadata": {
    "id": "89e43405-4d97-4c6c-8396-3e80ba3b7659"
   },
   "outputs": [],
   "source": [
    "@cudaq.kernel\n",
    "def trotter_step_pf1 (state: cudaq.State, dt: float, c_a: List[complex], \n",
    "                      c_b: List[complex], words_a: List[cudaq.pauli_word], words_b: List[cudaq.pauli_word]):\n",
    "    '''\n",
    "    The following kernel performs one trotter step for a Hamiltonian, split into two non-commuting parts.  \n",
    "    Each represented by: \n",
    "    c_a, c_b: List[complex] \n",
    "        a set of coefficients corresponding to the Pauli terms \n",
    "    words_a, words_b: List[complex]\n",
    "        a list of Pauli words specifying the operators\n",
    "    '''\n",
    "    \n",
    "    qubits = cudaq.qvector(state)\n",
    "\n",
    "    for i in range(len(c_a)):\n",
    "        exp_pauli(-c_a[i].real * dt, qubits, words_a[i])\n",
    "\n",
    "    for i in range(len(c_b)):\n",
    "        exp_pauli(-c_b[i].real * dt, qubits, words_b[i])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ce68136",
   "metadata": {},
   "source": [
    "### 2.5 Hamiltonian helpers\n",
    "\n",
    "The simulation loop expects explicit coefficient lists and Pauli words. The utility functions below wrap CUDA-Q’s iterator API so that kernels can consume plain Python lists. This also keeps the PF1 kernel agnostic to the specific model — any Hamiltonian decomposed into commuting partitions can reuse the same workflow.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "f338169f-0064-405c-9770-0a67ce669d52",
   "metadata": {
    "id": "f338169f-0064-405c-9770-0a67ce669d52"
   },
   "outputs": [],
   "source": [
    "'''\n",
    "Here, we extract coefficients and words from Hamiltonian for decomposition.\n",
    "'''\n",
    "def extractCoefficients(hamiltonian: cudaq.SpinOperator):\n",
    "    result = []\n",
    "    for term in hamiltonian:\n",
    "        result.append(term.evaluate_coefficient())\n",
    "        # print(term.evaluate_coefficient)\n",
    "    return result\n",
    "\n",
    "def extractWords(hamiltonian: cudaq.SpinOperator):\n",
    "    n_spins = hamiltonian.qubit_count\n",
    "    result = []\n",
    "    for term in hamiltonian:\n",
    "        result.append(term.get_pauli_word(n_spins))\n",
    "    return result"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ddafbee3-4a29-44bb-b193-33b4440ec770",
   "metadata": {
    "id": "ddafbee3-4a29-44bb-b193-33b4440ec770"
   },
   "source": [
    "## 3. Entanglement metrics\n",
    "\n",
    "To match the paper we track how much entanglement each subsystem grows during time evolution. For a pure state $|\\psi\\rangle$, the reduced density matrix of a kept subsystem $A$ is\n",
    "$$\n",
    "\\rho_A = \\operatorname{Tr}_B (|\\psi\\rangle\\langle\\psi|),\n",
    "$$\n",
    "where $B$ denotes the qubits we trace out. The von Neumann entropy\n",
    "$$\n",
    "S(\\rho_A) = - \\sum_i \\lambda_i \\log_2 \\lambda_i\n",
    "$$\n",
    "with eigenvalues $\\lambda_i$ provides a convenient single-number entanglement measure. When $S(\\rho_A)$ approaches its maximum the subsystem is almost maximally mixed, signalling strong entanglement with the rest of the chain. The helpers below implement the partial trace and entropy evaluation using CuPy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "740ba553-1caa-4407-951a-76bc19e0c238",
   "metadata": {
    "id": "740ba553-1caa-4407-951a-76bc19e0c238"
   },
   "outputs": [],
   "source": [
    "def partial_trace (psi: cp.ndarray, keep: List[int], n_qubits: int) -> cp.ndarray:\n",
    "    '''\n",
    "    The following function computes the reduced density matrix for a pure state of an n-qubit system. \n",
    "    psi: cp.ndarray \n",
    "        The full statevector of the system \n",
    "    keep: List[int] \n",
    "        The list of qubit indices to \"keep\" (subsystem A). The rest of the qubits will be traced out (sybsystem B). \n",
    "    n_qubits: int \n",
    "        The number of qubits \n",
    "    ''' \n",
    "    psi = psi / cp.linalg.norm(psi)\n",
    "\n",
    "    traced_out = [i for i in range(n_qubits) if i not in keep]\n",
    "    dim_keep = 2 ** len(keep)\n",
    "    dim_traced = 2 ** len(traced_out)\n",
    "    \n",
    "    psi = psi.reshape([2] * n_qubits)\n",
    "    perm = keep + traced_out\n",
    "    psi = psi.transpose(perm)\n",
    "    psi = psi.reshape(dim_keep, dim_traced)\n",
    "\n",
    "    rho_A = psi @ psi.conj().T\n",
    "    \n",
    "    return rho_A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "8d77d00c-012a-4436-95a4-1d7eab1b7bcb",
   "metadata": {
    "id": "8d77d00c-012a-4436-95a4-1d7eab1b7bcb"
   },
   "outputs": [],
   "source": [
    "def entropy(rho: cp.ndarray) -> float:\n",
    "    '''\n",
    "    The following function calculates the entropy of a given density matrix: \n",
    "    rho: cp.ndarray\n",
    "    ''' \n",
    "    eigvals = cp.linalg.eigvalsh(rho)\n",
    "    eigvals = eigvals[eigvals > 0]\n",
    "    entropy = -cp.sum(eigvals * cp.log2(eigvals))\n",
    "    return entropy.get() "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d756259",
   "metadata": {},
   "source": [
    "## 4. Simulation workflow\n",
    "\n",
    "With the building blocks in place, we compare the ideal evolution with the PF1 approximation for a 12-qubit chain. The helpers below construct the exact unitary $U = e^{-iH\\delta t}$ and step both trajectories forward so that we can track their divergence and the associated entanglement.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "f3e5fe13-55ed-4ead-a1fa-bc42cd64d397",
   "metadata": {
    "id": "f3e5fe13-55ed-4ead-a1fa-bc42cd64d397"
   },
   "outputs": [],
   "source": [
    "def ideal_U (H: cudaq.SpinOperator, dt: float) -> cp.ndarray:\n",
    "    '''\n",
    "    This function constructs the exact unitary time-evolution operator U = exp(-i*H*dt), given\n",
    "        H: cudaq.SpinOperator \n",
    "            The Hamiltonian operator \n",
    "        dt: float\n",
    "    '''\n",
    "    Hm = cp.asarray(H.to_matrix())\n",
    "    return expm(-1j * Hm * dt)\n",
    "\n",
    "def ideal_next_state(U: cp.ndarray, state: cp.ndarray) -> cp.ndarray:\n",
    "    '''\n",
    "    The following function apply the exact unitary time-evolution operator U to the current state. \n",
    "        U: cp.ndarray\n",
    "            Unitary time-evolution operator\n",
    "        state: cp.ndarray \n",
    "            Current quantum state \n",
    "    ''' \n",
    "    return U @ state"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "adecb633",
   "metadata": {},
   "source": [
    "### 4.1 Single-step Trotter error\n",
    "Following Eq. (3) of the paper we evaluate the norm difference between the ideal state and its PF1 counterpart after each segment:\n",
    "$$\n",
    "\\epsilon(t) = \\big\\|(U(\\delta t) - \\mathscr{U}_1(\\delta t))\\,|\\psi(t)\\rangle\\big\\|.\n",
    "$$\n",
    "Tracking $\\epsilon(t)$ alongside entanglement lets us reproduce the result highlighted in Figure 1a."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "492b5a9b-7ef8-4b3d-bef6-5b04acc8547c",
   "metadata": {
    "id": "492b5a9b-7ef8-4b3d-bef6-5b04acc8547c"
   },
   "outputs": [],
   "source": [
    "def trotter_error(psi_ideal, psi_trot):\n",
    "    return ((cp.linalg.norm(psi_ideal - psi_trot)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f094fd9",
   "metadata": {},
   "source": [
    "### 4.2 Dual trajectory update\n",
    "\n",
    "The driver below advances both trajectories for $r$ steps:\n",
    "\n",
    "- **Ideal evolution** multiplies the state by $U = e^{-iH\\delta t}$ at each step.\n",
    "- **PF1 evolution** applies the Trotter evolution $e^{-iA\\delta t}e^{-iB\\delta t}$ to the CUDA-Q state.\n",
    "\n",
    "Storing the intermediate states allows us to probe entanglement and Trotter error at every segment boundary, exactly as in the reference study."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7421e3e9-00fb-412a-9b01-4c9b60103e55",
   "metadata": {
    "id": "7421e3e9-00fb-412a-9b01-4c9b60103e55"
   },
   "outputs": [],
   "source": [
    "def ideal_simulation (n_spins: int, initial_state: cudaq.State, r: int, h_x_val: float, h_y_val: float, J_val: float):   \n",
    "    ideal_state = initial_state\n",
    "    t_state = initial_state\n",
    "\n",
    "    out_list_ideal = []\n",
    "    out_list_trotter = []\n",
    "\n",
    "    time_list = []\n",
    "\n",
    "    # dt is set as 0.1, as indicated in the paper \n",
    "    dt = 0.1\n",
    "    \n",
    "    # prepare hamiltonian for evolution \n",
    "    A, B, H = create_hamiltonian_qimf(n_spins, J=J_val, h_x=h_x_val, h_y=h_y_val)\n",
    "    c_a = extractCoefficients(A)\n",
    "    c_b = extractCoefficients(B)\n",
    "    words_a = extractWords(A)\n",
    "    words_b = extractWords(B)\n",
    "    \n",
    "    start_time = time.time()\n",
    "\n",
    "    U = ideal_U(H, dt) \n",
    "    psi_ideal = cp.asarray(ideal_state, dtype=cp.complex128) \n",
    "    \n",
    "    for step in range(1, r + 1):\n",
    "        psi_ideal = ideal_next_state(U, psi_ideal)\n",
    "        t_state = cudaq.get_state(trotter_step_pf1, t_state, dt, c_a, c_b, words_a, words_b)\n",
    "        \n",
    "        psi_t = cp.asarray(t_state, dtype=cp.complex128) \n",
    "\n",
    "        out_list_ideal.append(psi_ideal)\n",
    "        out_list_trotter.append(psi_t)\n",
    "        time_list.append(dt * step)\n",
    "        \n",
    "        # Compare trotter at one step \n",
    "        t_state = cudaq.State.from_data(cp.asnumpy(psi_ideal))\n",
    "    \n",
    "    total_time = time.time() - start_time\n",
    "    print(f\"Simulation done in {total_time:.2f} seconds\")\n",
    "\n",
    "    return out_list_ideal, out_list_trotter "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ccb0dffb",
   "metadata": {},
   "source": [
    "## 5. Reproducing the paper’s Figure 1a\n",
    "\n",
    "We now evolve the $n=12$ QIMF chain for $t=10$ with $r=100$ PF1 steps, exactly matching the parameters in the paper. At every segment we store the PF1 error $\\epsilon(t)$ and the entanglement entropy for contiguous subsystems of size 1 through 4."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "00b574ce-9cd2-413c-a285-552d04305183",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Simulation done in 14.51 seconds\n"
     ]
    }
   ],
   "source": [
    "# Expect ~4 seconds on an NVIDIA A100 GPU (dense ideal evolution dominates runtime).\n",
    "cudaq.set_target(\"qpp-cpu\")\n",
    "\n",
    "n = 12  # system size, number of qubits\n",
    "r = 100  # number of PF1 Trotter steps\n",
    "\n",
    "# initial state \n",
    "initial_state = cudaq.get_state(get_initial_state, n)\n",
    "\n",
    "# field strength and coupling strength \n",
    "h_x_val, h_y_val, J_val = 0.8090, 0.9045, 1\n",
    "\n",
    "# the simulation returns a collection of states from computing ideal and trotter evolutions \n",
    "ideal, trotter = ideal_simulation(n, initial_state, r, h_x_val, h_y_val, J_val)\n",
    "\n",
    "tro_error = []\n",
    "# compute trotter error by comparing the two states at each step\n",
    "for i in range(r): \n",
    "    tro = trotter_error(ideal[i], trotter[i])\n",
    "    tro_error.append(tro.get()) \n",
    "\n",
    "# calculate entanglement entropy for the ideal trajectory at each step\n",
    "entropy_list = []\n",
    "for i in range(r):\n",
    "    rdm1 = partial_trace(ideal[i], list(range(1)), n_qubits=n)\n",
    "    rdm2 = partial_trace(ideal[i], list(range(2)), n_qubits=n)\n",
    "    rdm3 = partial_trace(ideal[i], list(range(3)), n_qubits=n)\n",
    "    rdm4 = partial_trace(ideal[i], list(range(4)), n_qubits=n)\n",
    "\n",
    "    temp = [entropy(rdm1), entropy(rdm2), entropy(rdm3), entropy(rdm4)]\n",
    "    entropy_list.append(temp)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "66aff8d4",
   "metadata": {},
   "source": [
    "### 5.1 Visualising the joint behaviour\n",
    "\n",
    "The plot below overlays the PF1 Trotter error with the entropy of increasingly large contiguous subsystems. This reproduces the qualitative behaviour of Figure 1a from [*Entanglement Accelerates Quantum Simulation*](https://www.nature.com/articles/s41567-025-02945-2):\n",
    "\n",
    "- the PF1 error (blue) drops within the first few units of simulated time,\n",
    "- each entropy curve (coloured) rises toward the maximally mixed value for its subsystem size, and\n",
    "- larger subsystems take longer to saturate, mirroring the light-cone spreading of entanglement.\n",
    "\n",
    "The key takeaway is that higher entanglement correlates with *lower* observed Trotter error, counter to the classical intuition that entanglement makes simulation harder."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "b0296a92-5b9e-48ae-88c6-798bc88b316f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 87604 characters]",
      "text/plain": [
       "<Figure size 800x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "times = np.linspace(0, 10, r)\n",
    "\n",
    "entropy_values = [\n",
    "    [step[k] for step in entropy_list]\n",
    "    for k in range(4)\n",
    "]\n",
    "labels = [r\"$S_{1}$\", r\"$S_{1,2}$\", r\"$S_{1,2,3}$\", r\"$S_{1,2,3,4}$\"]\n",
    "colors = [\"tab:orange\", \"tab:green\", \"tab:red\", \"tab:purple\"]\n",
    "\n",
    "fig, ax1 = plt.subplots(figsize=(8, 5))\n",
    "ax2 = ax1.twinx()\n",
    "\n",
    "# PF1 error trend (left axis)\n",
    "err_line = ax1.plot(times, tro_error, color=\"tab:blue\", label=\"PF1 Trotter error\", linewidth=2)\n",
    "ax1.set_ylabel(\"Trotter error\", color=\"tab:blue\")\n",
    "ax1.set_xlabel(\"Time $t$\")\n",
    "ax1.tick_params(axis=\"y\", labelcolor=\"tab:blue\")\n",
    "ax1.set_xlim(0, 10)\n",
    "ax1.grid(alpha=0.3)\n",
    "\n",
    "# Entropy growth (right axis)\n",
    "entropy_lines = []\n",
    "for vals, label, color in zip(entropy_values, labels, colors):\n",
    "    line = ax2.plot(times, vals, label=label, color=color, linewidth=2)\n",
    "    entropy_lines.extend(line)\n",
    "\n",
    "ax2.set_ylabel(\"Subsystem entanglement entropy\")\n",
    "ax2.set_ylim(0, 4.2)\n",
    "ax2.tick_params(axis=\"y\")\n",
    "\n",
    "# Combined legend\n",
    "lines = err_line + entropy_lines\n",
    "labels = [line.get_label() for line in lines]\n",
    "ax2.legend(lines, labels, loc=\"upper right\")\n",
    "\n",
    "ax1.set_title(\"PF1 error decreases as entanglement spreads\")\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c78dbd20",
   "metadata": {},
   "source": [
    "### 5.2 Interpreting the result\n",
    "\n",
    "Intuitively, the entangled state locally looks like an ensemble of random states, leading to average-case Trotter errors that shrink as entanglement increases. \n",
    "The tutorial demonstrates that CUDA-Q can capture this effect faithfully for realistic system sizes. Try experimenting with:\n",
    "\n",
    "- different subsystem cuts to probe entanglement propagation in more detail,\n",
    "- higher-order product formulas (PF2, PF4) to compare error scalings, or\n",
    "- alternative initial product states to see how quickly the system enters the high-entanglement regime.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d44b14d",
   "metadata": {},
   "source": [
    "## 6. References and further reading\n",
    "- [Entanglement Accelerates Quantum Simulation](https://arxiv.org/abs/2406.02379)\n",
    "- [Spin-Hamiltonian Simulation Using CUDA-Q](https://nvidia.github.io/cuda-quantum/latest/applications/python/hamiltonian_simulation.html)\n"
   ]
  }
 ],
 "metadata": {
  "accelerator": "GPU",
  "colab": {
   "gpuType": "A100",
   "machine_shape": "hm",
   "provenance": []
  },
  "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.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
