{
  "cells": [
    {
      "attachments": {},
      "cell_type": "markdown",
      "id": "a94f5c51-348c-4fb0-9488-70ad4a19ace0",
      "metadata": {},
      "source": [
        "---\n",
        "title: Quantum simulation\n",
        "description: This course is about quantum simulation, including Trotterization for the transverse field Ising Hamiltonian\n",
        "---\n",
        "\n",
        "{/* cspell:ignore sharex */}\n",
        "\n",
        "# Quantum simulation\n",
        "\n",
        "<Admonition type=\"note\">\n",
        "  Yukio Kawashima (May 30, 2024)\n",
        "\n",
        "  [Download the pdf](https://ibm.ent.box.com/s/kzzsmxhw38vph1ioohczaet53euwi310) of the original lecture. Note that some code snippets might become deprecated since these are static images.\n",
        "\n",
        "  *Approximate QPU time to run this experiment is 7 seconds.*\n",
        "\n",
        "  (This notebook is mostly taken from a now-deprecated [tutorial notebook](https://github.com/qiskit-community/qiskit-algorithms/blob/main/docs/tutorials/13_trotterQRTE.ipynb) for Qiskit Algorithms.)\n",
        "</Admonition>\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "aef7b02a-2b09-4409-85a0-8cf8346df430",
      "metadata": {},
      "source": [
        "## 1. Introduction\n",
        "\n",
        "As a real time evolution technique, Trotterization consists in the successive application of a quantum gate or gates, chosen to approximate the time evolution of a system for a time slice. Following from the Schrödinger equation, the time evolution of a system initially in the state $\\vert\\psi(0)\\rangle$ takes the form:\n",
        "\n",
        "$$\n",
        "\\vert \\psi(t) \\rangle = e^{-i H t} \\vert \\psi(0) \\rangle \\text{,}\n",
        "$$\n",
        "\n",
        "where $H$ is the time-independent Hamiltonian governing the system. We consider a Hamiltonian that can be written as a weighted sum of Pauli terms $H=\\sum_j a_j P_j$, with $P_j$ representing a tensor product of Pauli terms acting on $n$ qubits. In particular, these Pauli terms might commute with one another, or they might not. Given a state at time $t=0$, how do we obtain the system's state at a later time $|\\psi(t)\\rangle$ using a quantum computer? The exponential of an operator can be most easily understood through its Taylor series:\n",
        "\n",
        "$$\n",
        "e^{-i H t} = 1-iHt-\\frac{1}{2}H^2t^2+...\n",
        "$$\n",
        "\n",
        "Some very basic exponentials, like $e^{iZ}$ can be implemented easily on quantum computers using a compact set of quantum gates. Most Hamiltonians of interest will not have just a single term, but will instead have many terms. Note what happens if $H = H_1+H_2$:\n",
        "\n",
        "$$\n",
        "e^{-i H t} = 1-i(H_1+H_2)t-\\frac{1}{2}(H_1+H_2)^2t^2+...\n",
        "$$\n",
        "\n",
        "When $H_1$ and $H_2$ commute, we have the familiar case (which is also true for numbers, and variables $a$ and $b$ below):\n",
        "\n",
        "$$\n",
        "e^{-i (a+b) t} = e^{-i a t}e^{-i b t}\n",
        "$$\n",
        "\n",
        "But when operators do not commute, terms cannot be rearranged in the Taylor series to simplify in this way. Thus, expressing complicated Hamiltonians in quantum gates is a challenge.\n",
        "\n",
        "One solution is to consider very small time $t$, such that the first-order term in the Taylor expansion dominates. Under that assumption:\n",
        "\n",
        "$$\n",
        "e^{-i (H_1+H_2) t} \\approx 1-i(H_1+H_2)t \\approx (1-i H_1 t)(1-i H_2 t) \\approx e^{-i H_1 t}e^{-i H_2 t}\n",
        "$$\n",
        "\n",
        "Of course, we might need to evolve our state for a longer time. That is accomplished by using many such small steps in time. This process is called Trotterization:\n",
        "\n",
        "$$\n",
        "\\vert \\psi(t) \\rangle \\approx \\left(\\prod_j e^{-i a_j P_j t/r} \\right)^r \\vert\\psi(0) \\rangle \\text{,}\n",
        "$$\n",
        "\n",
        "Here $t/r$ is the time slice (evolution step) that we are choosing. As a result, a gate to be applied $r$ times is created. A smaller timestep leads to a more accurate approximation. However, this also leads to deeper circuits which, in practice, leads to more error accumulation (a non-negligible concern on near-term quantum devices).\n",
        "\n",
        "Today, we will study the time evolution of the [Ising model](https://en.wikipedia.org/wiki/Ising_model) on linear lattices of $N=2$ and $N=6$ sites. These lattices consist of an array of spins $\\sigma_i$ that interact only with their nearest neighbors. These spins can have two orientations: $\\uparrow$ and $\\downarrow$, which correspond to a magnetization of $+1$ and $-1$ respectively.\n",
        "\n",
        "$$\n",
        "H = - J \\sum_{i=0}^{N-2} Z_i Z_{i+1} - h \\sum_{i=0}^{N-1} X_i  \\text{,}\n",
        "$$\n",
        "\n",
        "where $J$ describes the interaction energy, and $h$ the magnitude of an external field (in the x-direction above, but we will modify this). Let us write this expression using Pauli matrices, and considering that the external field has an angle $\\alpha$ with respect to the transversal direction,\n",
        "\n",
        "$$\n",
        "H = -J \\sum_{i=0}^{N-2} Z_i Z_{i+1} -h \\sum_{i=0}^{N-1} (\\sin\\alpha Z_i + \\cos\\alpha X_i) \\text{.}\n",
        "$$\n",
        "\n",
        "This Hamiltonian is useful in that it allows us to easily study the effects of an external field. In the computational basis, the system will be encoded as follows:\n",
        "\n",
        "|       Quantum state      |             Spin representation            |\n",
        "| :----------------------: | :----------------------------------------: |\n",
        "| $\\lvert 0 0 0 0 \\rangle$ |     $\\uparrow\\uparrow\\uparrow\\uparrow$     |\n",
        "| $\\lvert 1 0 0 0 \\rangle$ |    $\\downarrow\\uparrow\\uparrow\\uparrow$    |\n",
        "|         $\\ldots$         |                  $\\ldots$                  |\n",
        "| $\\lvert 1 1 1 1 \\rangle$ | $\\downarrow\\downarrow\\downarrow\\downarrow$ |\n",
        "\n",
        "We will start investigating the time evolution of such a quantum system. More specifically, we will visualize the time-evolution of certain properties of the system like magnetization.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "8d26b146-12b4-4882-a0c8-3f7b5b15f04d",
      "metadata": {},
      "source": [
        "### 1.1 Requirements\n",
        "\n",
        "Before starting this tutorial, be sure you have the following installed:\n",
        "\n",
        "* Qiskit SDK v1.2 or later ( `pip install qiskit` )\n",
        "* Qiskit Runtime v0.30 or later ( `pip install qiskit-ibm-runtime` )\n",
        "* Numpy v1.24.1 or later \\< 2 ( `pip install numpy` )\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "e73c8a4a-a540-4637-b179-670da6b1a6e6",
      "metadata": {},
      "source": [
        "### 1.2 Import the libraries\n",
        "\n",
        "Note that some libraries that might be useful (MatrixExponential, QDrift) are included even though they are not used in this current notebook. You can try them out if you have time!\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "id": "2edce0da-4b37-411f-ac13-191590d99e1b",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "'2.0.2'"
            ]
          },
          "execution_count": 1,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# Check the version of Qiskit\n",
        "import qiskit\n",
        "\n",
        "qiskit.__version__"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "id": "803f9fab-abab-482e-b8f5-90d817f2452f",
      "metadata": {
        "scrolled": true,
        "tags": []
      },
      "outputs": [],
      "source": [
        "# Import the qiskit library\n",
        "import numpy as np\n",
        "import matplotlib.pylab as plt\n",
        "import warnings\n",
        "\n",
        "from qiskit import QuantumCircuit\n",
        "from qiskit.circuit.library import PauliEvolutionGate\n",
        "from qiskit.primitives import StatevectorEstimator\n",
        "from qiskit.quantum_info import Statevector, SparsePauliOp\n",
        "from qiskit.synthesis import (\n",
        "    SuzukiTrotter,\n",
        "    LieTrotter,\n",
        ")\n",
        "from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager\n",
        "\n",
        "from qiskit_ibm_runtime import QiskitRuntimeService, SamplerV2\n",
        "\n",
        "warnings.filterwarnings(\"ignore\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "64dbeb76-bb61-469d-bdce-65f4b5cb65fb",
      "metadata": {},
      "source": [
        "## 2. Mapping your problem\n",
        "\n",
        "### 2.1 Defining the transverse-field Ising Hamiltonian\n",
        "\n",
        "We here consider the 1-D transverse-field Ising model.\n",
        "\n",
        "First, we will create a function that takes in the system parameters $N$, $J$, $h$ and $\\alpha$, and returns our Hamiltonian as a `SparsePauliOp`. A [SparsePauliOp](/docs/api/qiskit/qiskit.quantum_info.SparsePauliOp) is a sparse representation of an operator in terms of weighted [Pauli](/docs/api/qiskit/qiskit.quantum_info.Pauli) terms.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "id": "440e0b22-d221-485c-845f-9fefa1d10ae6",
      "metadata": {},
      "outputs": [],
      "source": [
        "def get_hamiltonian(nqubits, J, h, alpha):\n",
        "    # List of Hamiltonian terms as 3-tuples containing\n",
        "    # (1) the Pauli string,\n",
        "    # (2) the qubit indices corresponding to the Pauli string,\n",
        "    # (3) the coefficient.\n",
        "    ZZ_tuples = [(\"ZZ\", [i, i + 1], -J) for i in range(0, nqubits - 1)]\n",
        "    Z_tuples = [(\"Z\", [i], -h * np.sin(alpha)) for i in range(0, nqubits)]\n",
        "    X_tuples = [(\"X\", [i], -h * np.cos(alpha)) for i in range(0, nqubits)]\n",
        "\n",
        "    # We create the Hamiltonian as a SparsePauliOp, via the method\n",
        "    # `from_sparse_list`, and multiply by the interaction term.\n",
        "    hamiltonian = SparsePauliOp.from_sparse_list(\n",
        "        [*ZZ_tuples, *Z_tuples, *X_tuples], num_qubits=nqubits\n",
        "    )\n",
        "    return hamiltonian.simplify()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "b56952ca-f354-4719-8bd9-674f71f681a4",
      "metadata": {},
      "source": [
        "#### Define the Hamiltonian\n",
        "\n",
        "The system that we now consider has a size of $N=6$, $J=0.2$, $h=1.2$ and $\\alpha=\\frac{\\pi}{8.0}$ as an example.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "id": "0eca59f6-8aec-4c58-a22b-9a237b5787b8",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "SparsePauliOp(['IIIIZZ', 'IIIZZI', 'IIZZII', 'IZZIII', 'ZZIIII', 'IIIIIZ', 'IIIIZI', 'IIIZII', 'IIZIII', 'IZIIII', 'ZIIIII', 'IIIIIX', 'IIIIXI', 'IIIXII', 'IIXIII', 'IXIIII', 'XIIIII'],\n",
              "              coeffs=[-0.2       +0.j, -0.2       +0.j, -0.2       +0.j, -0.2       +0.j,\n",
              " -0.2       +0.j, -0.45922012+0.j, -0.45922012+0.j, -0.45922012+0.j,\n",
              " -0.45922012+0.j, -0.45922012+0.j, -0.45922012+0.j, -1.10865544+0.j,\n",
              " -1.10865544+0.j, -1.10865544+0.j, -1.10865544+0.j, -1.10865544+0.j,\n",
              " -1.10865544+0.j])"
            ]
          },
          "execution_count": 4,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "n_qubits = 6\n",
        "\n",
        "hamiltonian = get_hamiltonian(nqubits=n_qubits, J=0.2, h=1.2, alpha=np.pi / 8.0)\n",
        "hamiltonian"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "5a9fa555-e4ad-4e4c-acc1-bd7302b9361e",
      "metadata": {},
      "source": [
        "### 2.2 Set the parameters of the time-evolution simulation\n",
        "\n",
        "Here we will consider three different Trotterization techniques:\n",
        "\n",
        "* Lie–Trotter (first order)\n",
        "* second-order Suzuki–Trotter\n",
        "* fourth-order Suzuki–Trotter\n",
        "\n",
        "The latter two will be used in the exercise, and in the appendix.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "id": "9749e548-3238-4d0a-bcbf-c7b6699942d3",
      "metadata": {},
      "outputs": [],
      "source": [
        "num_timesteps = 60\n",
        "evolution_time = 30.0\n",
        "dt = evolution_time / num_timesteps\n",
        "product_formula_lt = LieTrotter()\n",
        "product_formula_st2 = SuzukiTrotter(order=2)\n",
        "product_formula_st4 = SuzukiTrotter(order=4)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "d910b0e0-5e95-4899-b433-d13710eb4ada",
      "metadata": {},
      "source": [
        "### 2.3 Prepare the quantum circuit 1 (Initial state)\n",
        "\n",
        "Create an initial state. Here we will start with a spin configuration of $\\uparrow\\uparrow\\downarrow\\downarrow\\uparrow\\uparrow$\n",
        ".\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "id": "c9d78393-1e0c-46bb-b909-ada0cadfd59d",
      "metadata": {
        "scrolled": true,
        "tags": []
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/learning/images/courses/utility-scale-quantum-computing/quantum-simulation/extracted-outputs/b338806d-67bf-47ae-9d03-944bde3e2c99-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "execution_count": 6,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "initial_circuit = QuantumCircuit(n_qubits)\n",
        "initial_circuit.prepare_state(\"001100\")\n",
        "# Change reps and see the difference when you decompose the circuit\n",
        "initial_circuit.decompose(reps=1).draw(\"mpl\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "2e31b0bd-ebbb-4e74-b95d-a3ae2e3fff62",
      "metadata": {},
      "source": [
        "### 2.4 Prepare the quantum circuit 2 (Single circuit for time evolution)\n",
        "\n",
        "We here construct a circuit for a single time step using Lie–Trotter.\n",
        "\n",
        "The Lie product formula (first-order) is implemented in the [LieTrotter](/docs/api/qiskit/qiskit.synthesis.LieTrotter) class. A first-order formula consists of the approximation stated in the introduction, where the matrix exponential of a sum is approximated by a product of matrix exponentials:\n",
        "\n",
        "$$\n",
        "e^{H_1+H_2} \\approx e^{H_1} e^{H_2}\n",
        "$$\n",
        "\n",
        "As previously mentioned, very deep circuits lead to the accumulation of errors, and cause problems for modern quantum computers. Because two-qubit gates have higher error rates than single-qubit gates, a quantity of particular interest is the two-qubit circuit depth. What really matters is the two-qubit circuit depth after transpilation (since that is the circuit the quantum computer actually executes). But let us get in the habit of counting the operations for this circuit, even now using the simulator.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "id": "e1ce17c2-02f2-4270-b271-55661ca14112",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "\n",
            "Trotter step with Lie-Trotter\n",
            "-----------------------------\n",
            "Depth: 17\n",
            "Gate count: 27\n",
            "Nonlocal gate count: 10\n",
            "Gate breakdown: U3: 12, CX: 10, U1: 5\n",
            "\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              "<Image src=\"/learning/images/courses/utility-scale-quantum-computing/quantum-simulation/extracted-outputs/252c2734-653d-4c8e-af6b-d3173845de88-1.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "execution_count": 8,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "single_step_evolution_gates_lt = PauliEvolutionGate(\n",
        "    hamiltonian, dt, synthesis=product_formula_lt\n",
        ")\n",
        "single_step_evolution_lt = QuantumCircuit(n_qubits)\n",
        "single_step_evolution_lt.append(\n",
        "    single_step_evolution_gates_lt, single_step_evolution_lt.qubits\n",
        ")\n",
        "\n",
        "print(\n",
        "    f\"\"\"\n",
        "Trotter step with Lie-Trotter\n",
        "-----------------------------\n",
        "Depth: {single_step_evolution_lt.decompose(reps=3).depth()}\n",
        "Gate count: {len(single_step_evolution_lt.decompose(reps=3))}\n",
        "Nonlocal gate count: {single_step_evolution_lt.decompose(reps=3).num_nonlocal_gates()}\n",
        "Gate breakdown: {\", \".join([f\"{k.upper()}: {v}\" for k, v in single_step_evolution_lt.decompose(reps=3).count_ops().items()])}\n",
        "\"\"\"\n",
        ")\n",
        "single_step_evolution_lt.decompose(reps=3).draw(\"mpl\", fold=-1)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "6710865f-8e47-485f-bab8-e281d37bfbcd",
      "metadata": {},
      "source": [
        "### 2.5 Set the operators to measure\n",
        "\n",
        "Let us define a *magnetization operator* $\\sum_i \\langle Z_i \\rangle / N$, and a *mean spin correlation operator* $\\sum_i \\langle Z_i Z_{i+1} \\rangle/ (N - 1)$.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "id": "008188d1-047a-40b9-aaea-cf6baa42c434",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "magnetization :  SparsePauliOp(['IIIIIZ', 'IIIIZI', 'IIIZII', 'IIZIII', 'IZIIII', 'ZIIIII'],\n",
            "              coeffs=[0.16666667+0.j, 0.16666667+0.j, 0.16666667+0.j, 0.16666667+0.j,\n",
            " 0.16666667+0.j, 0.16666667+0.j])\n",
            "correlation :  SparsePauliOp(['IIIIZZ', 'IIIZZI', 'IIZZII', 'IZZIII', 'ZZIIII'],\n",
            "              coeffs=[0.2+0.j, 0.2+0.j, 0.2+0.j, 0.2+0.j, 0.2+0.j])\n"
          ]
        }
      ],
      "source": [
        "magnetization = (\n",
        "    SparsePauliOp.from_sparse_list(\n",
        "        [(\"Z\", [i], 1.0) for i in range(0, n_qubits)], num_qubits=n_qubits\n",
        "    )\n",
        "    / n_qubits\n",
        ")\n",
        "correlation = SparsePauliOp.from_sparse_list(\n",
        "    [(\"ZZ\", [i, i + 1], 1.0) for i in range(0, n_qubits - 1)], num_qubits=n_qubits\n",
        ") / (n_qubits - 1)\n",
        "print(\"magnetization : \", magnetization)\n",
        "print(\"correlation : \", correlation)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "9a54c9da-d0e6-423d-aa4f-be0808159fa3",
      "metadata": {},
      "source": [
        "### 2.6 Perform time-evolution simulation\n",
        "\n",
        "We will monitor the energy (expectation value of the Hamiltonian), magnetization (expectation value of the magnetization operator), and mean spin correlation (expectation value of the mean spin correlation operator). Qiskit's `StatevectorEstimator` (EstimatorV2) primitive estimates expectation values of observables, $\\langle\\psi\\vert\\hat{O}\\vert\\psi\\rangle$.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "id": "0045114c-9f92-4e10-92d5-9cd478628781",
      "metadata": {},
      "outputs": [],
      "source": [
        "# Initiate the circuit\n",
        "evolved_state = QuantumCircuit(initial_circuit.num_qubits)\n",
        "# Start from the initial spin configuration\n",
        "evolved_state.append(initial_circuit, evolved_state.qubits)\n",
        "# Initiate Estimator (V2)\n",
        "estimator = StatevectorEstimator()\n",
        "# Set number of shots\n",
        "shots = 10000\n",
        "# Translate the precision required from the number of shots\n",
        "precision = np.sqrt(1 / shots)\n",
        "energy_list = []\n",
        "mag_list = []\n",
        "corr_list = []\n",
        "# Estimate expectation values for t=0.0\n",
        "job = estimator.run(\n",
        "    [(evolved_state, [hamiltonian, magnetization, correlation])], precision=precision\n",
        ")\n",
        "# Get estimated expectation values\n",
        "evs = job.result()[0].data.evs\n",
        "energy_list.append(evs[0])\n",
        "mag_list.append(evs[1])\n",
        "corr_list.append(evs[2])\n",
        "# Start time evolution\n",
        "for n in range(num_timesteps):\n",
        "    # Expand the circuit to describe delta-t\n",
        "    evolved_state.append(single_step_evolution_gates_lt, evolved_state.qubits)\n",
        "    # Estimate expectation values at delta-t\n",
        "    job = estimator.run(\n",
        "        [(evolved_state, [hamiltonian, magnetization, correlation])],\n",
        "        precision=precision,\n",
        "    )\n",
        "    # Retrieve results (expectation values)\n",
        "    evs = job.result()[0].data.evs\n",
        "    energy_list.append(evs[0])\n",
        "    mag_list.append(evs[1])\n",
        "    corr_list.append(evs[2])\n",
        "# Transform the list of expectation values (at each time step) to arrays\n",
        "energy_array = np.array(energy_list)\n",
        "mag_array = np.array(mag_list)\n",
        "corr_array = np.array(corr_list)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "91003bdc-0f5d-47cb-b664-2a483062a7ac",
      "metadata": {},
      "source": [
        "### 2.7 Plot the time evolution of the observables\n",
        "\n",
        "We plot the expectation values we measured against time.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "id": "3e3cd58b-0705-4cdf-a63b-01d1e3f2315b",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "Text(0.5, 0.98, 'Observable evolution')"
            ]
          },
          "execution_count": 11,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "text/plain": [
              "<Image src=\"/learning/images/courses/utility-scale-quantum-computing/quantum-simulation/extracted-outputs/9c0f4b39-801a-448c-9dce-8a6d7c83fe76-1.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "fig, axes = plt.subplots(3, sharex=True)\n",
        "times = np.linspace(0, evolution_time, num_timesteps + 1)  # includes initial state\n",
        "axes[0].plot(\n",
        "    times,\n",
        "    energy_array,\n",
        "    label=\"First order\",\n",
        "    marker=\"x\",\n",
        "    c=\"darkmagenta\",\n",
        "    ls=\"-\",\n",
        "    lw=0.8,\n",
        ")\n",
        "axes[1].plot(\n",
        "    times, mag_array, label=\"First order\", marker=\"x\", c=\"darkmagenta\", ls=\"-\", lw=0.8\n",
        ")\n",
        "axes[2].plot(\n",
        "    times, corr_array, label=\"First order\", marker=\"x\", c=\"darkmagenta\", ls=\"-\", lw=0.8\n",
        ")\n",
        "axes[0].set_ylabel(\"Energy\")\n",
        "axes[1].set_ylabel(\"Magnetization\")\n",
        "axes[2].set_ylabel(\"Mean spin correlation\")\n",
        "axes[2].set_xlabel(\"Time\")\n",
        "fig.suptitle(\"Observable evolution\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "481f1e48-0650-4f53-b3ce-7d0b8c1e171e",
      "metadata": {},
      "source": [
        "## 3. Exercise 1. Perform simulation using second-order Suzuki–Trotter\n",
        "\n",
        "Now let's try performing simulation with second-order Suzuki–Trotter following the example of Lie–Trotter shown above.\n",
        "\n",
        "The second-order Suzuki-Trotter can be used in Qiskit by means of the [SuzukiTrotter class](/docs/api/qiskit/qiskit.synthesis.SuzukiTrotter). Using this formula, a second order decomposition is:\n",
        "\n",
        "$$\n",
        "e^{H_1+H_2} \\approx e^{H_1/2}e^{H_2}e^{H_1/2}\n",
        "$$\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "1e47dd47-ce42-4c9b-88f3-7a44145edf31",
      "metadata": {},
      "source": [
        "### 3.1 Construct a circuit for a single time step\n",
        "\n",
        "Use product\\_formula\\_st2 (SuzukiTrotter(order=2)) and construct a circuit for a single time step using second-order Suzuki–Trotter. Also, count the number of gates and depth of the circuit and compare with Lie–Trotter.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 13,
      "id": "8b6389e3-9c44-4618-9d19-7854f122715a",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "\n",
            "Trotter step with second-order Suzuki-Trotter\n",
            "-----------------------------\n",
            "Depth: 34\n",
            "Gate count: 53\n",
            "Nonlocal gate count: 20\n",
            "Gate breakdown: U3: 23, CX: 20, U1: 10\n",
            "\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              "<Image src=\"/learning/images/courses/utility-scale-quantum-computing/quantum-simulation/extracted-outputs/9da4ff65-1989-4fc0-b63b-5e166a77f2d5-1.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "execution_count": 13,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# Modify the line below (Use PauliEvolutionGate)\n",
        "single_step_evolution_gates_st2 = PauliEvolutionGate(\n",
        "    hamiltonian, dt, synthesis=product_formula_st2\n",
        ")\n",
        "single_step_evolution_st2 = QuantumCircuit(n_qubits)\n",
        "single_step_evolution_st2.append(\n",
        "    single_step_evolution_gates_st2, single_step_evolution_st2.qubits\n",
        ")\n",
        "# Let us print some stats\n",
        "print(\n",
        "    f\"\"\"\n",
        "Trotter step with second-order Suzuki-Trotter\n",
        "-----------------------------\n",
        "Depth: {single_step_evolution_st2.decompose(reps=3).depth()}\n",
        "Gate count: {len(single_step_evolution_st2.decompose(reps=3))}\n",
        "Nonlocal gate count: {single_step_evolution_st2.decompose(reps=3).num_nonlocal_gates()}\n",
        "Gate breakdown: {\", \".join([f\"{k.upper()}: {v}\" for k, v in single_step_evolution_st2.decompose(reps=3).count_ops().items()])}\n",
        "\"\"\"\n",
        ")\n",
        "single_step_evolution_st2.decompose(reps=2).draw(\"mpl\", fold=-1)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "1748cc5d-360b-415f-8281-ab908d71c057",
      "metadata": {},
      "source": [
        "### 3.2 Perform time-evolution simulation\n",
        "\n",
        "Perform time evolution using second-order Suzuki–Trotter.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 14,
      "id": "59ceae10-5675-4af2-b88f-6b928df35606",
      "metadata": {},
      "outputs": [],
      "source": [
        "# Initiate the circuit\n",
        "evolved_state = QuantumCircuit(initial_circuit.num_qubits)\n",
        "# Start from the initial spin configuration\n",
        "evolved_state.append(initial_circuit, evolved_state.qubits)\n",
        "# Initiate Estimator (V2)\n",
        "estimator = StatevectorEstimator()\n",
        "# Set number of shots\n",
        "shots = 10000\n",
        "# Translate the precision required from the number of shots\n",
        "precision = np.sqrt(1 / shots)\n",
        "energy_list_st2 = []\n",
        "mag_list_st2 = []\n",
        "corr_list_st2 = []\n",
        "# Estimate expectation values for t=0.0\n",
        "job = estimator.run(\n",
        "    [(evolved_state, [hamiltonian, magnetization, correlation])], precision=precision\n",
        ")\n",
        "# Get estimated expectation values\n",
        "evs = job.result()[0].data.evs\n",
        "energy_list_st2.append(evs[0])\n",
        "mag_list_st2.append(evs[1])\n",
        "corr_list_st2.append(evs[2])\n",
        "# Start time evolution\n",
        "for n in range(num_timesteps):\n",
        "    # Expand the circuit to describe delta-t\n",
        "    evolved_state.append(single_step_evolution_gates_st2, evolved_state.qubits)\n",
        "    # Estimate expectation values at delta-t\n",
        "    job = estimator.run(\n",
        "        [(evolved_state, [hamiltonian, magnetization, correlation])],\n",
        "        precision=precision,\n",
        "    )\n",
        "    # Retrieve results (expectation values)\n",
        "    evs = job.result()[0].data.evs\n",
        "    energy_list_st2.append(evs[0])\n",
        "    mag_list_st2.append(evs[1])\n",
        "    corr_list_st2.append(evs[2])\n",
        "# Transform the list of expectation values (at each time step) to arrays\n",
        "energy_array_st2 = np.array(energy_list_st2)\n",
        "mag_array_st2 = np.array(mag_list_st2)\n",
        "corr_array_st2 = np.array(corr_list_st2)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "cd6468c8-bd78-4043-ab54-b6036178399f",
      "metadata": {},
      "source": [
        "### 3.3 Plot the second-order Suzuki–Trotter results\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 15,
      "id": "b333b560-a6ae-4ab0-a900-c098e545c9c8",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/learning/images/courses/utility-scale-quantum-computing/quantum-simulation/extracted-outputs/f859672f-cef3-4c2c-b087-9d36fa8162ff-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "execution_count": 15,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "axes[0].plot(\n",
        "    times,\n",
        "    energy_array_st2,\n",
        "    label=\"Second Order\",\n",
        "    marker=\"x\",\n",
        "    c=\"limegreen\",\n",
        "    ls=\"-\",\n",
        "    lw=0.8,\n",
        ")\n",
        "axes[1].plot(\n",
        "    times,\n",
        "    mag_array_st2,\n",
        "    label=\"Second Order\",\n",
        "    marker=\"x\",\n",
        "    c=\"limegreen\",\n",
        "    ls=\"-\",\n",
        "    lw=0.8,\n",
        ")\n",
        "axes[2].plot(\n",
        "    times,\n",
        "    corr_array_st2,\n",
        "    label=\"Second Order\",\n",
        "    marker=\"x\",\n",
        "    c=\"limegreen\",\n",
        "    ls=\"-\",\n",
        "    lw=0.8,\n",
        ")\n",
        "\n",
        "# Replace the legend\n",
        "# legend.remove()\n",
        "legend = fig.legend(\n",
        "    *axes[0].get_legend_handles_labels(),\n",
        "    bbox_to_anchor=(1.0, 0.5),\n",
        "    loc=\"center left\",\n",
        "    framealpha=0.5,\n",
        ")\n",
        "fig"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "e12ad5f0-eff1-4238-97fb-593e7caf67ba",
      "metadata": {},
      "source": [
        "### 3.4 Compare with exact results\n",
        "\n",
        "The data below is the precomputed exact results from the classical computer.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 16,
      "id": "912e69a9-fa10-4c34-ab8c-6d970a461f8f",
      "metadata": {},
      "outputs": [],
      "source": [
        "exact_times = np.array(\n",
        "    [\n",
        "        0.0,\n",
        "        0.3,\n",
        "        0.6,\n",
        "        0.8999999999999999,\n",
        "        1.2,\n",
        "        1.5,\n",
        "        1.7999999999999998,\n",
        "        2.1,\n",
        "        2.4,\n",
        "        2.6999999999999997,\n",
        "        3.0,\n",
        "        3.3,\n",
        "        3.5999999999999996,\n",
        "        3.9,\n",
        "        4.2,\n",
        "        4.5,\n",
        "        4.8,\n",
        "        5.1,\n",
        "        5.3999999999999995,\n",
        "        5.7,\n",
        "        6.0,\n",
        "        6.3,\n",
        "        6.6,\n",
        "        6.8999999999999995,\n",
        "        7.199999999999999,\n",
        "        7.5,\n",
        "        7.8,\n",
        "        8.1,\n",
        "        8.4,\n",
        "        8.7,\n",
        "        9.0,\n",
        "        9.299999999999999,\n",
        "        9.6,\n",
        "        9.9,\n",
        "        10.2,\n",
        "        10.5,\n",
        "        10.799999999999999,\n",
        "        11.1,\n",
        "        11.4,\n",
        "        11.7,\n",
        "        12.0,\n",
        "        12.299999999999999,\n",
        "        12.6,\n",
        "        12.9,\n",
        "        13.2,\n",
        "        13.5,\n",
        "        13.799999999999999,\n",
        "        14.1,\n",
        "        14.399999999999999,\n",
        "        14.7,\n",
        "        15.0,\n",
        "        15.299999999999999,\n",
        "        15.6,\n",
        "        15.899999999999999,\n",
        "        16.2,\n",
        "        16.5,\n",
        "        16.8,\n",
        "        17.099999999999998,\n",
        "        17.4,\n",
        "        17.7,\n",
        "        18.0,\n",
        "        18.3,\n",
        "        18.599999999999998,\n",
        "        18.9,\n",
        "        19.2,\n",
        "        19.5,\n",
        "        19.8,\n",
        "        20.099999999999998,\n",
        "        20.4,\n",
        "        20.7,\n",
        "        21.0,\n",
        "        21.3,\n",
        "        21.599999999999998,\n",
        "        21.9,\n",
        "        22.2,\n",
        "        22.5,\n",
        "        22.8,\n",
        "        23.099999999999998,\n",
        "        23.4,\n",
        "        23.7,\n",
        "        24.0,\n",
        "        24.3,\n",
        "        24.599999999999998,\n",
        "        24.9,\n",
        "        25.2,\n",
        "        25.5,\n",
        "        25.8,\n",
        "        26.099999999999998,\n",
        "        26.4,\n",
        "        26.7,\n",
        "        27.0,\n",
        "        27.3,\n",
        "        27.599999999999998,\n",
        "        27.9,\n",
        "        28.2,\n",
        "        28.5,\n",
        "        28.799999999999997,\n",
        "        29.099999999999998,\n",
        "        29.4,\n",
        "        29.7,\n",
        "        30.0,\n",
        "    ]\n",
        ")\n",
        "exact_energy = np.array(\n",
        "    [\n",
        "        -1.1184402376762155,\n",
        "        -1.1184402376762157,\n",
        "        -1.1184402376762157,\n",
        "        -1.1184402376762148,\n",
        "        -1.1184402376762153,\n",
        "        -1.1184402376762155,\n",
        "        -1.1184402376762148,\n",
        "        -1.118440237676216,\n",
        "        -1.118440237676216,\n",
        "        -1.1184402376762166,\n",
        "        -1.1184402376762148,\n",
        "        -1.118440237676216,\n",
        "        -1.1184402376762153,\n",
        "        -1.1184402376762148,\n",
        "        -1.118440237676217,\n",
        "        -1.118440237676215,\n",
        "        -1.1184402376762161,\n",
        "        -1.1184402376762157,\n",
        "        -1.118440237676217,\n",
        "        -1.1184402376762161,\n",
        "        -1.1184402376762137,\n",
        "        -1.1184402376762161,\n",
        "        -1.1184402376762161,\n",
        "        -1.118440237676218,\n",
        "        -1.1184402376762155,\n",
        "        -1.1184402376762166,\n",
        "        -1.1184402376762155,\n",
        "        -1.1184402376762137,\n",
        "        -1.1184402376762186,\n",
        "        -1.1184402376762215,\n",
        "        -1.1184402376762148,\n",
        "        -1.118440237676216,\n",
        "        -1.1184402376762166,\n",
        "        -1.1184402376762148,\n",
        "        -1.1184402376762121,\n",
        "        -1.1184402376762166,\n",
        "        -1.1184402376762181,\n",
        "        -1.1184402376762137,\n",
        "        -1.1184402376762148,\n",
        "        -1.1184402376762193,\n",
        "        -1.1184402376762108,\n",
        "        -1.1184402376762144,\n",
        "        -1.118440237676217,\n",
        "        -1.1184402376762197,\n",
        "        -1.1184402376762153,\n",
        "        -1.1184402376762161,\n",
        "        -1.1184402376762184,\n",
        "        -1.1184402376762126,\n",
        "        -1.118440237676214,\n",
        "        -1.118440237676214,\n",
        "        -1.1184402376762161,\n",
        "        -1.118440237676212,\n",
        "        -1.1184402376762164,\n",
        "        -1.118440237676217,\n",
        "        -1.1184402376762121,\n",
        "        -1.1184402376762157,\n",
        "        -1.1184402376762212,\n",
        "        -1.1184402376762217,\n",
        "        -1.1184402376762206,\n",
        "        -1.118440237676222,\n",
        "        -1.1184402376762166,\n",
        "        -1.118440237676212,\n",
        "        -1.1184402376762137,\n",
        "        -1.11844023767622,\n",
        "        -1.1184402376762206,\n",
        "        -1.118440237676219,\n",
        "        -1.1184402376762153,\n",
        "        -1.1184402376762164,\n",
        "        -1.118440237676209,\n",
        "        -1.1184402376762144,\n",
        "        -1.1184402376762161,\n",
        "        -1.118440237676216,\n",
        "        -1.1184402376762173,\n",
        "        -1.118440237676214,\n",
        "        -1.1184402376762093,\n",
        "        -1.1184402376762184,\n",
        "        -1.1184402376762126,\n",
        "        -1.118440237676213,\n",
        "        -1.1184402376762195,\n",
        "        -1.1184402376762095,\n",
        "        -1.1184402376762075,\n",
        "        -1.1184402376762197,\n",
        "        -1.1184402376762141,\n",
        "        -1.1184402376762146,\n",
        "        -1.1184402376762184,\n",
        "        -1.118440237676218,\n",
        "        -1.1184402376762224,\n",
        "        -1.118440237676219,\n",
        "        -1.118440237676218,\n",
        "        -1.1184402376762206,\n",
        "        -1.1184402376762168,\n",
        "        -1.118440237676221,\n",
        "        -1.118440237676218,\n",
        "        -1.1184402376762148,\n",
        "        -1.1184402376762106,\n",
        "        -1.1184402376762173,\n",
        "        -1.118440237676216,\n",
        "        -1.118440237676216,\n",
        "        -1.1184402376762113,\n",
        "        -1.1184402376762275,\n",
        "        -1.1184402376762195,\n",
        "    ]\n",
        ")\n",
        "exact_magnetization = np.array(\n",
        "    [\n",
        "        0.3333333333333333,\n",
        "        0.26316769633415005,\n",
        "        0.0912947227110664,\n",
        "        -0.09317712543141576,\n",
        "        -0.20391854332115245,\n",
        "        -0.19318196655046493,\n",
        "        -0.06411527074401464,\n",
        "        0.12558269854206197,\n",
        "        0.28252754464640606,\n",
        "        0.3264196194042506,\n",
        "        0.2361586169847769,\n",
        "        0.060894367906122224,\n",
        "        -0.10842387093076275,\n",
        "        -0.18636359582538073,\n",
        "        -0.1338364343947887,\n",
        "        0.020284606520827753,\n",
        "        0.19151142743926025,\n",
        "        0.2905341647678381,\n",
        "        0.2723014646745304,\n",
        "        0.15147481733047252,\n",
        "        -0.008179102877790292,\n",
        "        -0.1242999208732406,\n",
        "        -0.1372529247781061,\n",
        "        -0.04083616185958952,\n",
        "        0.11066094926716476,\n",
        "        0.23140661570567636,\n",
        "        0.2587109403786205,\n",
        "        0.1868237670027325,\n",
        "        0.061201779383143744,\n",
        "        -0.051391248969654205,\n",
        "        -0.09843899603365061,\n",
        "        -0.061297056158849166,\n",
        "        0.04199010081939773,\n",
        "        0.15861461430963147,\n",
        "        0.22336830674799552,\n",
        "        0.20179555623336537,\n",
        "        0.11407111438609417,\n",
        "        0.01609419104778282,\n",
        "        -0.04239611796730001,\n",
        "        -0.04249123521065924,\n",
        "        0.008850291714888112,\n",
        "        0.08780898151558082,\n",
        "        0.1561486776507056,\n",
        "        0.17627348772811832,\n",
        "        0.13870676179652253,\n",
        "        0.07205869195282538,\n",
        "        0.018300003064909465,\n",
        "        0.0001095640839572417,\n",
        "        0.015157929316037586,\n",
        "        0.05077755280969454,\n",
        "        0.09245534457650838,\n",
        "        0.12206907551110702,\n",
        "        0.12284950557969157,\n",
        "        0.09570215398601932,\n",
        "        0.06294378255078983,\n",
        "        0.045503313813986014,\n",
        "        0.043389819499542556,\n",
        "        0.046725117769796744,\n",
        "        0.054956411358382404,\n",
        "        0.0713814528253614,\n",
        "        0.08743689703248492,\n",
        "        0.08951216359166674,\n",
        "        0.07878386475305985,\n",
        "        0.06955669116405788,\n",
        "        0.06639892435963689,\n",
        "        0.05890378761746903,\n",
        "        0.04541796525844558,\n",
        "        0.0414221088331947,\n",
        "        0.05499634106912299,\n",
        "        0.07409418836014572,\n",
        "        0.08371859070160165,\n",
        "        0.08211623987959302,\n",
        "        0.07615055161378328,\n",
        "        0.06702584458783024,\n",
        "        0.051891407742740085,\n",
        "        0.038049378383635625,\n",
        "        0.03825614149768043,\n",
        "        0.054183218463525695,\n",
        "        0.0753534475741016,\n",
        "        0.08853147112587295,\n",
        "        0.08767917178542013,\n",
        "        0.07709383184439536,\n",
        "        0.06308595032042386,\n",
        "        0.0498812359204284,\n",
        "        0.04299040064096167,\n",
        "        0.04769159891460652,\n",
        "        0.06483569572288776,\n",
        "        0.08698035745435016,\n",
        "        0.10047391641776235,\n",
        "        0.09747255683203637,\n",
        "        0.08098863187287358,\n",
        "        0.05959496723987331,\n",
        "        0.04383882265040485,\n",
        "        0.04232138798062125,\n",
        "        0.05720514169944535,\n",
        "        0.08201306299870219,\n",
        "        0.10274898262000469,\n",
        "        0.10707552455080133,\n",
        "        0.09210856128265357,\n",
        "        0.06379922105742579,\n",
        "        0.03624325103307953,\n",
        "    ]\n",
        ")\n",
        "exact_correlation = np.array(\n",
        "    [\n",
        "        0.2,\n",
        "        0.1247704225763532,\n",
        "        0.01943938494098705,\n",
        "        0.03854917181332821,\n",
        "        0.11196616231067426,\n",
        "        0.0906546700356683,\n",
        "        0.01629373561896267,\n",
        "        0.011352652889791095,\n",
        "        0.0636185676540077,\n",
        "        0.09543834437789013,\n",
        "        0.10058518161011307,\n",
        "        0.11829217731417431,\n",
        "        0.1397812224038133,\n",
        "        0.12316460402216707,\n",
        "        0.08541383059335775,\n",
        "        0.06144846844403662,\n",
        "        0.020246372880505827,\n",
        "        -0.02693683090021662,\n",
        "        0.003919250903281282,\n",
        "        0.1117419430168554,\n",
        "        0.19676155181256794,\n",
        "        0.18594408880783336,\n",
        "        0.1002673802566004,\n",
        "        0.03821525827438024,\n",
        "        0.04485205090247377,\n",
        "        0.05348102743040269,\n",
        "        0.03160026140008638,\n",
        "        0.033437649060464834,\n",
        "        0.10486939975320728,\n",
        "        0.20249469538955758,\n",
        "        0.19735507621013149,\n",
        "        0.0553097261765083,\n",
        "        -0.04889114490131667,\n",
        "        0.011685690974970964,\n",
        "        0.11705971535823065,\n",
        "        0.11681165998194759,\n",
        "        0.06637091239560744,\n",
        "        0.10936684225958895,\n",
        "        0.20225454101061405,\n",
        "        0.16284420833341812,\n",
        "        -0.0025823294931362067,\n",
        "        -0.0763416631752919,\n",
        "        0.02985268630418397,\n",
        "        0.15234468006771007,\n",
        "        0.14606385406970995,\n",
        "        0.0935341856492092,\n",
        "        0.12325421854361143,\n",
        "        0.17130422930386324,\n",
        "        0.10383730044042278,\n",
        "        -0.031333159406547614,\n",
        "        -0.05241572078596815,\n",
        "        0.07722509925347705,\n",
        "        0.17642188574256007,\n",
        "        0.12765340239966838,\n",
        "        0.06309968945093776,\n",
        "        0.11574687130499339,\n",
        "        0.16978282647206913,\n",
        "        0.0736143632571229,\n",
        "        -0.05356602733119409,\n",
        "        -0.0009649396796768892,\n",
        "        0.15921620111869142,\n",
        "        0.17760366431811037,\n",
        "        0.04736297330213485,\n",
        "        0.012122870263181897,\n",
        "        0.13268065586830521,\n",
        "        0.1728473023503636,\n",
        "        0.03999259331072221,\n",
        "        -0.036997053070222885,\n",
        "        0.06951528580242439,\n",
        "        0.1769169993516561,\n",
        "        0.12290448295710298,\n",
        "        0.012897784654866427,\n",
        "        0.02859435620982225,\n",
        "        0.12895847695150875,\n",
        "        0.13629536955485938,\n",
        "        0.05394621059822597,\n",
        "        0.02298040588184324,\n",
        "        0.07036499900317271,\n",
        "        0.11706448623132719,\n",
        "        0.10435285842074606,\n",
        "        0.055721236329964965,\n",
        "        0.04676334743672697,\n",
        "        0.08417924910022263,\n",
        "        0.10611161955304965,\n",
        "        0.089304171047322,\n",
        "        0.06098589533081194,\n",
        "        0.06314519797488709,\n",
        "        0.09431492621892917,\n",
        "        0.09667836915967139,\n",
        "        0.0651298357290882,\n",
        "        0.05176966009147416,\n",
        "        0.06727229484222669,\n",
        "        0.08871788283607947,\n",
        "        0.09907054249093444,\n",
        "        0.09785167773502176,\n",
        "        0.09277216140054353,\n",
        "        0.07520999642062785,\n",
        "        0.05894392248382922,\n",
        "        0.07236135251622376,\n",
        "        0.08608284185200156,\n",
        "        0.07282922961856123,\n",
        "    ]\n",
        ")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 17,
      "id": "d626c298-49e5-4cba-9cf4-ba0935b13f95",
      "metadata": {
        "scrolled": true
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/learning/images/courses/utility-scale-quantum-computing/quantum-simulation/extracted-outputs/c95e61c9-229c-4eed-b240-e6540f80956a-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "execution_count": 17,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "axes[0].plot(exact_times, exact_energy, c=\"k\", ls=\":\", label=\"Exact\")\n",
        "axes[1].plot(exact_times, exact_magnetization, c=\"k\", ls=\":\", label=\"Exact\")\n",
        "axes[2].plot(exact_times, exact_correlation, c=\"k\", ls=\":\", label=\"Exact\")\n",
        "# Replace the legend\n",
        "legend.remove()\n",
        "# Select the labels of only the first axis\n",
        "legend = fig.legend(\n",
        "    *axes[0].get_legend_handles_labels(),\n",
        "    bbox_to_anchor=(1.0, 0.5),\n",
        "    loc=\"center left\",\n",
        "    framealpha=0.5,\n",
        ")\n",
        "fig.tight_layout()\n",
        "fig"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "6c72f95e-2259-4b76-9e06-ad296e0cfa7d",
      "metadata": {},
      "source": [
        "## 4. Executing on the quantum hardware\n",
        "\n",
        "We next run the time-evolution simulation on the quantum hardware. We will work on a smaller problem, lattice size N=2. We vary the $\\alpha$ parameter and see the difference in dynamics of the wavefunction.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "9a5bb8bf-9a96-4464-8bcc-eeaee193c23a",
      "metadata": {},
      "source": [
        "### 4.1  Step 1. Map classical inputs to a quantum problem\n",
        "\n",
        "Pick the initial setting of the simulation:\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 19,
      "id": "40ac8e04-7eab-4e92-8964-6987d707cc04",
      "metadata": {},
      "outputs": [],
      "source": [
        "n_qubits_2 = 2\n",
        "dt_2 = 1.6\n",
        "product_formula = LieTrotter(reps=1)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "1ad01369-146b-41c6-a557-5abae8eb0122",
      "metadata": {},
      "source": [
        "Then set the initial circuit:\n",
        "\n",
        "The initial spin configuration will be \"down-up\"\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 20,
      "id": "b7ed1249-1a69-4d62-97a3-102cf905a295",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/learning/images/courses/utility-scale-quantum-computing/quantum-simulation/extracted-outputs/f33796b1-7511-46b0-93ff-ee6e6188c412-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "execution_count": 20,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# We prepare an initial state ↓↑ (10).\n",
        "# Note that Statevector and SparsePauliOp interpret the qubits from right to left\n",
        "initial_circuit_2 = QuantumCircuit(n_qubits_2)\n",
        "initial_circuit_2.prepare_state(\"10\")\n",
        "# Change reps and see the difference when you decompose the circuit\n",
        "initial_circuit_2.decompose(reps=1).draw(\"mpl\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "b383a0e6-bebe-4d47-865d-50b0258a8fd9",
      "metadata": {},
      "source": [
        "Now compute the reference value using an ideal statevector simulator.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 21,
      "id": "d15f5cef-59a4-4946-9a3f-6bee81dd2587",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<matplotlib.legend.Legend at 0x11c816590>"
            ]
          },
          "execution_count": 21,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "text/plain": [
              "<Image src=\"/learning/images/courses/utility-scale-quantum-computing/quantum-simulation/extracted-outputs/8469258b-bb82-4e46-b3d6-43aee59d9474-1.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "bar_width = 0.1\n",
        "# initial_state = Statevector.from_label(\"10\")\n",
        "final_time = 1.6\n",
        "eps = 1e-5\n",
        "\n",
        "# We create the list of angles in radians, with a small epsilon\n",
        "# the exactly longitudinal field, which would present no dynamics at all\n",
        "alphas = np.linspace(-np.pi / 2 + eps, np.pi / 2 - eps, 5)\n",
        "\n",
        "for i, alpha in enumerate(alphas):\n",
        "    evolved_state_2 = QuantumCircuit(initial_circuit_2.num_qubits)\n",
        "    evolved_state_2.append(initial_circuit_2, evolved_state_2.qubits)\n",
        "    hamiltonian_2 = get_hamiltonian(nqubits=2, J=0.2, h=1.0, alpha=alpha)\n",
        "    single_step_evolution_gates_2 = PauliEvolutionGate(\n",
        "        hamiltonian_2, dt_2, synthesis=product_formula\n",
        "    )\n",
        "    evolved_state_2.append(single_step_evolution_gates_2, evolved_state_2.qubits)\n",
        "    evolved_state_2 = Statevector(evolved_state_2)\n",
        "    # Dictionary of probabilities\n",
        "    amplitudes_dict = evolved_state_2.probabilities_dict()\n",
        "    labels = list(amplitudes_dict.keys())\n",
        "    values = list(amplitudes_dict.values())\n",
        "    # Convert angle to degrees\n",
        "    alpha_str = f\"$\\\\alpha={int(np.round(alpha * 180 / np.pi))}^\\\\circ$\"\n",
        "    plt.bar(np.arange(4) + i * bar_width, values, bar_width, label=alpha_str, alpha=0.7)\n",
        "\n",
        "plt.xticks(np.arange(4) + 2 * bar_width, labels)\n",
        "plt.xlabel(\"Measurement\")\n",
        "plt.ylabel(\"Probability\")\n",
        "plt.suptitle(\n",
        "    f\"Measurement probabilities at $t={final_time}$, for various field angles $\\\\alpha$\\n\"\n",
        "    f\"Initial state: 10, Linear lattice of size $L=2$\"\n",
        ")\n",
        "plt.legend()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "0588757d-5dcf-4ee3-8d62-1b14070156eb",
      "metadata": {},
      "source": [
        "We have prepared a system initially with a sequence of spins $\\downarrow\\uparrow$, which corresponds to $\\vert\\psi(0)\\rangle = \\vert10\\rangle$. After letting it evolve for $t=1.6$ under a transversal field ($\\alpha=0^\\circ$), we are almost guaranteed to measure $\\uparrow\\downarrow$, that is, have a spin swap. (Note that the labels are interpreted from right to left). If the field is longitudinal ($\\alpha=\\pm90^\\circ$), we will have no evolution, therefore we will measure the system as it was initially prepared, $\\downarrow\\uparrow$. With intermediate angles, at $\\alpha=\\pm45^\\circ$, we will be able to measure all combinations will different probabilities, being a spin swap the most likely with a probability of 67%.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "36ef1a7e-1ab7-48ed-89ad-efbef819c871",
      "metadata": {},
      "source": [
        "#### Construct circuit for HW experiment\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 24,
      "id": "ef73cba1-3469-472d-a100-26fb1e103840",
      "metadata": {},
      "outputs": [],
      "source": [
        "circuit_list = []\n",
        "for i, alpha in enumerate(alphas):\n",
        "    evolved_state_2 = QuantumCircuit(initial_circuit_2.num_qubits)\n",
        "    evolved_state_2.append(initial_circuit_2, evolved_state_2.qubits)\n",
        "    hamiltonian_2 = get_hamiltonian(nqubits=2, J=0.2, h=1.0, alpha=alpha)\n",
        "    single_step_evolution_gates_2 = PauliEvolutionGate(\n",
        "        hamiltonian_2, dt_2, synthesis=product_formula\n",
        "    )\n",
        "    evolved_state_2.append(single_step_evolution_gates_2, evolved_state_2.qubits)\n",
        "    evolved_state_2.measure_all()\n",
        "    circuit_list.append(evolved_state_2)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "0cf2010d-2621-46fa-9828-f06979faad60",
      "metadata": {},
      "source": [
        "### 4.2 Step 2. Optimize for target hardware\n",
        "\n",
        "We specify a backend.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "0a7734a6-6ec3-48a8-91cb-8995a15ada88",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "'ibm_strasbourg'"
            ]
          },
          "execution_count": 25,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "service = QiskitRuntimeService()\n",
        "backend = service.least_busy(operational=True, simulator=False)\n",
        "backend.name"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "1d4be180-a241-4e84-a730-66e1a4fc9e9a",
      "metadata": {},
      "source": [
        "Then we transpile the circuit for the selected backend.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 26,
      "id": "b154db84-439c-40d0-bf14-d2a6b584b593",
      "metadata": {},
      "outputs": [],
      "source": [
        "pm = generate_preset_pass_manager(backend=backend, optimization_level=3)\n",
        "circuit_isa = pm.run(circuit_list)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "ce78ebe6-1ed5-470d-a000-490cfb6cc9e3",
      "metadata": {},
      "source": [
        "Check out the circuit.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 28,
      "id": "04de26ef-4854-4004-846b-3e2d7a2a3d0e",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/learning/images/courses/utility-scale-quantum-computing/quantum-simulation/extracted-outputs/64d9166e-296a-4a6d-a6d1-0b6fdd49399b-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "execution_count": 28,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "circuit_isa[1].draw(\"mpl\", idle_wires=False)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "d08cdec5-9127-434b-b435-7525d7379158",
      "metadata": {},
      "source": [
        "### 4.3 Step 3. Execute with Qiskit Runtime primitives\n",
        "\n",
        "Qiskit's `Sampler` (V2) primitive provides the counts of measured bitstrings.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 29,
      "id": "836b8fb2-6799-4c93-8c00-278c4dbb2744",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "job id: d13pswfmya70008ek070\n"
          ]
        }
      ],
      "source": [
        "sampler = SamplerV2(mode=backend)\n",
        "job = sampler.run(circuit_isa)\n",
        "job_id = job.job_id()\n",
        "print(\"job id:\", job_id)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "382d29d0-c9fc-452e-9584-d4276df6d93a",
      "metadata": {},
      "source": [
        "Save the results\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 30,
      "id": "8a2fc35c-8801-49bd-b715-b2906e3e422e",
      "metadata": {},
      "outputs": [],
      "source": [
        "results = job.result()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "b67e9e11-794d-4159-be82-d1935988f5af",
      "metadata": {},
      "source": [
        "### 4.4 Step 4. Post-process results\n",
        "\n",
        "Construct the histogram of the bitstrings, which corresponds to analyzing the wavefunction, and compare them with the ideal values shown above.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "336e6b34-deee-467f-bd8f-6d1a8ed1ccba",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<matplotlib.legend.Legend at 0x11d7af990>"
            ]
          },
          "execution_count": 32,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "text/plain": [
              "<Image src=\"/learning/images/courses/utility-scale-quantum-computing/quantum-simulation/extracted-outputs/1f4e5fbd-1994-4f4d-8d35-ddf1a37664b7-1.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "list_temp = [\"00\", \"01\", \"10\", \"11\"]\n",
        "\n",
        "for i, alpha in enumerate(alphas):\n",
        "    # Dictionary of probabilities\n",
        "    amplitudes_dict = results[i].data.meas.get_counts()\n",
        "    values = []\n",
        "    for str_temp in list_temp:\n",
        "        values.append(\n",
        "            amplitudes_dict[str_temp] / 4096.0\n",
        "        )  # divided by default number of shots\n",
        "    # Convert angle to degrees\n",
        "    alpha_str = f\"$\\\\alpha={int(np.round(alpha * 180 / np.pi))}^\\\\circ$\"\n",
        "    plt.bar(np.arange(4) + i * bar_width, values, bar_width, label=alpha_str, alpha=0.7)\n",
        "\n",
        "plt.xticks(np.arange(4) + 2 * bar_width, labels)\n",
        "plt.xlabel(\"Measurement\")\n",
        "plt.ylabel(\"Probabilities\")\n",
        "plt.suptitle(\n",
        "    f\"Measurement probabilities at $t={final_time}$, for various field angles $\\\\alpha$\\n\"\n",
        "    f\"Initial state: 10, Linear lattice of size $L=2$\"\n",
        ")\n",
        "plt.legend()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "8960b35a-ab79-48ad-a9a9-6024ad767699",
      "metadata": {},
      "source": [
        "We here show an example for constructing a circuit using higher-order (fourth-order) Suzuki–Trotter.\n",
        "Now let's try constructing a circuit simulation with fourth-order Suzuki–Trotter following the examples shown above.\n",
        "\n",
        "The fourth-order Suzuki–Trotter can be used in Qiskit by means of the [SuzukiTrotter class](/docs/api/qiskit/qiskit.synthesis.SuzukiTrotter). The fourth-order can be evaluated using the following recursion relation. Note that the order of Suzuki–Trotter is denoted as \"2k\" in the following equations.\n",
        "\n",
        "$$\n",
        "\\hat{U}_{ST(2k)}\\left(t\\right) = \\left[ \\hat{U}_{ST(2k-2)}\\left(p_k t\\right) \\right]^2 \\hat{U}_{ST(2k-2)}\\left( (1- 4 p_k) t\\right)\\left[ \\hat{U}_{ST(2k-2)}\\left(p_k t\\right) \\right]^2\n",
        "$$\n",
        "\n",
        "$$\n",
        "p_k = 1 / \\left(4-4^{\\frac{1}{2k-1}}\\right)\n",
        "$$\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "67041754-db9b-457a-b326-059e0e565d20",
      "metadata": {},
      "source": [
        "#### Construct a circuit for a single time step\n",
        "\n",
        "Use product\\_formula\\_st4 (SuzukiTrotter(order=4)) and construct a circuit for a single time step using fourth-order Suzuki–Trotter. Also, count the number of gates and depth of the circuit and compare with Lie–Trotter and second-order Suzuki–Trotter.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 33,
      "id": "be1a7a71-1f42-47c2-bee5-0f510f54d064",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "\n",
            "Trotter step with second-order Suzuki-Trotter\n",
            "-----------------------------\n",
            "Depth: 170\n",
            "Gate count: 265\n",
            "Nonlocal gate count: 100\n",
            "Gate breakdown: U3: 115, CX: 100, U1: 50\n",
            "\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              "<Image src=\"/learning/images/courses/utility-scale-quantum-computing/quantum-simulation/extracted-outputs/e20ea669-82d3-42a7-89f6-bab7d41689f1-1.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "execution_count": 33,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# Modify the line below (Use PauliEvolutionGate)\n",
        "single_step_evolution_gates_st4 = PauliEvolutionGate(\n",
        "    hamiltonian, dt, synthesis=product_formula_st4\n",
        ")\n",
        "single_step_evolution_st4 = QuantumCircuit(n_qubits)\n",
        "single_step_evolution_st4.append(\n",
        "    single_step_evolution_gates_st4, single_step_evolution_st4.qubits\n",
        ")\n",
        "# Let us print some stats\n",
        "print(\n",
        "    f\"\"\"\n",
        "Trotter step with second-order Suzuki-Trotter\n",
        "-----------------------------\n",
        "Depth: {single_step_evolution_st4.decompose(reps=3).depth()}\n",
        "Gate count: {len(single_step_evolution_st4.decompose(reps=3))}\n",
        "Nonlocal gate count: {single_step_evolution_st4.decompose(reps=3).num_nonlocal_gates()}\n",
        "Gate breakdown: {\", \".join([f\"{k.upper()}: {v}\" for k, v in single_step_evolution_st4.decompose(reps=3).count_ops().items()])}\n",
        "\"\"\"\n",
        ")\n",
        "single_step_evolution_st4.decompose(reps=2).draw(\"mpl\", fold=-1)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 34,
      "id": "87364170-60b0-4d7b-9dd1-a12676b1f80f",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "'2.0.2'"
            ]
          },
          "execution_count": 34,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# Check Qiskit version\n",
        "import qiskit\n",
        "\n",
        "qiskit.__version__"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "id": "a1b8767d",
      "source": "© IBM Corp., 2017-2026"
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3"
    },
    "widgets": {
      "application/vnd.jupyter.widget-state+json": {
        "state": {},
        "version_major": 2,
        "version_minor": 0
      }
    }
  },
  "nbformat": 4,
  "nbformat_minor": 4
}