{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "578810ce",
      "metadata": {},
      "source": [
        "---\n",
        "title: Spin and charge dynamics in the Fermi-Hubbard model\n",
        "description: Simulate real-time dynamics of a 1D Fermi-Hubbard chain on a QPU and watch spin and charge wavefronts spread at different speeds\n",
        "---\n",
        "\n",
        "{/* cspell:ignore Trotterized Hubbard Wigner holon holons spinon spinons Néel Fock Nup Ndn qctrl pcolormesh polyfit linewidth nanmax sharey Hartnett Najafi Khindanov Biercuk Baum Essler Frahm Gohmann Klumper Korepin Schutzman arrowprops arrowstyle fractionalizes frameon markerfacecolor udud wavefronts */}\n",
        "\n",
        "# Spin and charge dynamics in the Fermi-Hubbard model\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "a19187d8",
      "metadata": {},
      "source": [
        "*Usage estimate: 5 minutes on a Heron r3 processor. (NOTE: This is an estimate only. Your runtime might vary.)*\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "c3ee27fd",
      "metadata": {},
      "source": [
        "## Background\n",
        "\n",
        "The Fermi-Hubbard model describes electrons hopping on a lattice with an energy penalty for double occupancy. On a one-dimensional chain of $L$ sites with open boundaries, the Hamiltonian is defined as follows:\n",
        "\n",
        "$$\n",
        "H = -\\sum_{i=0}^{L-2} \\sum_{\\sigma \\in \\{\\uparrow,\\downarrow\\}} t_{h,i} \\left( c^\\dagger_{i\\sigma} c_{i+1,\\sigma} + \\mathrm{h.c.} \\right) + \\sum_{i=0}^{L-1} U_i\\, n_{i\\uparrow} n_{i\\downarrow} - \\sum_{i=0}^{L-1} \\sum_\\sigma \\mu_i\\, n_{i\\sigma},\n",
        "$$\n",
        "\n",
        "where $c^\\dagger_{i\\sigma}$ creates a fermion of spin $\\sigma$ on site $i$ and $n_{i\\sigma} = c^\\dagger_{i\\sigma} c_{i\\sigma}$ counts them. The three terms are the hopping amplitude $t_{h,i}$ on the bond between sites $i$ and $i+1$, the on-site interaction $U_i$, and the chemical potential $\\mu_i$. The tutorial measures energy in units of $t_h$ and time in units of $1/t_h$.\n",
        "\n",
        "Real-time dynamics of this model are hard to compute classically. Exact diagonalization stops near $L \\lesssim 14$ sites, and after a quench, the entanglement entropy of a tensor-network solver grows linearly in time, so the bond dimension needed for fixed accuracy grows exponentially.\n",
        "\n",
        "In one dimension, the elementary excitations are spinons, which carry spin and no charge, and holons, which carry charge and no spin. The two propagate at different speeds, so a localized disturbance splits into a fast front and a slow one. This effect is called spin-charge separation, and this tutorial measures it [\\[1\\]](#references).\n",
        "\n",
        "This tutorial uses the Dynamics Simulator Qiskit Function, which takes the Hamiltonian, an initial state, a Trotter schedule, and a list of observables, and returns expectation values. It handles the Jordan-Wigner mapping, Trotterization, qubit layout, and error suppression internally, so there are no circuits to build or transpile in this tutorial.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "4da1a887",
      "metadata": {},
      "source": [
        "## Requirements\n",
        "\n",
        "Before you start this tutorial, make sure that you have the following packages installed:\n",
        "\n",
        "* Qiskit SDK v2.3.0 or later (`pip install qiskit`)\n",
        "* Qiskit Runtime v0.50.0 or later (`pip install qiskit-ibm-runtime`)\n",
        "* Qiskit Functions Catalog v0.21.0 or later (`pip install qiskit-ibm-catalog`)\n",
        "* NumPy v2.5.1 or later (`pip install numpy`)\n",
        "* Matplotlib v3.10.3 or later (`pip install matplotlib`)\n",
        "* Q-CTRL Visualizer v11.0.0 or later (`pip install qctrl-visualizer`)\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "782b6eaa",
      "metadata": {},
      "source": [
        "## Setup\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "9b4b38bb",
      "metadata": {},
      "source": [
        "First, authenticate using your [IBM Quantum API key](/). Then, select the Qiskit Function as follows. (This code assumes that you already [saved your account](/docs/guides/functions-get-started#install-qiskit-functions-catalog-client) to your local environment.)\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 27,
      "id": "4be3b0bd",
      "metadata": {},
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "from scipy import stats\n",
        "from qiskit_ibm_catalog import QiskitFunctionsCatalog\n",
        "from qiskit_ibm_runtime import QiskitRuntimeService\n",
        "import qctrlvisualizer as qv\n",
        "\n",
        "plt.style.use(qv.get_qctrl_style())"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "242e3d57",
      "metadata": {},
      "outputs": [],
      "source": [
        "catalog = QiskitFunctionsCatalog(channel=\"ibm_quantum_platform\")\n",
        "\n",
        "# Access Function\n",
        "sim_qf = catalog.load(\"q-ctrl/dynamics-simulator\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "0bc24110",
      "metadata": {},
      "source": [
        "## Step 1: Map classical inputs to a quantum problem\n",
        "\n",
        "### Define the Hamiltonian\n",
        "\n",
        "The function takes the model as a dictionary with three keys. Each one accepts either a single float, which applies uniformly along the chain, or a list that specifies the value site by site. Lists have $L-1$ entries for `hopping_amplitude`, because an open chain has one fewer bond than sites, and $L$ entries for `on_site_interaction` and `chemical_potential`.\n",
        "\n",
        "This tutorial uses a uniform chain of 31 sites with a repulsive interaction $U/t_h = 8$. A uniform chemical potential contributes only a global phase when the particle number is fixed, so this example sets it to zero. The chemical potential matters only when it varies from site to site, for example, to add disorder or a single impurity.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "id": "738e696a",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "31 sites, 62 qubits, U/t_h = 8\n"
          ]
        }
      ],
      "source": [
        "site_count = 31\n",
        "hopping_amplitude = 1.0\n",
        "on_site_interaction = 8.0\n",
        "chemical_potential = 0.0\n",
        "\n",
        "model = {\n",
        "    \"hopping_amplitude\": [hopping_amplitude] * (site_count - 1),\n",
        "    \"on_site_interaction\": [on_site_interaction] * site_count,\n",
        "    \"chemical_potential\": [chemical_potential] * site_count,\n",
        "}\n",
        "\n",
        "print(\n",
        "    f\"{site_count} sites, {2 * site_count} qubits, U/t_h = {on_site_interaction:.0f}\"\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "ddf42054",
      "metadata": {},
      "source": [
        "### Set the Trotter schedule\n",
        "\n",
        "The `simulation` dictionary holds the time step ($\\Delta t$) and the number of Trotter steps, and the total evolution time is their product. Larger steps reach later times at fixed circuit depth but incur more Trotter error, which grows with $U$ because the Hamiltonian terms commute less. The published hardware study used $\\Delta t = 0.15$ for $U/t_h = 8$ [\\[2\\]](#references), and this tutorial uses the same value. Circuit depth grows linearly in the step count and does not depend on $L$ at all.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "id": "6a3366eb",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Total evolution time: t = 4.50 (units of 1/t_h)\n"
          ]
        }
      ],
      "source": [
        "time_step = 0.15\n",
        "step_count = 30\n",
        "\n",
        "simulation = {\"time_step\": time_step, \"step_count\": step_count}\n",
        "\n",
        "print(\n",
        "    f\"Total evolution time: t = {time_step * step_count:.2f} (units of 1/t_h)\"\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "1381c77a",
      "metadata": {},
      "source": [
        "### Prepare the initial state\n",
        "\n",
        "Each site holds two spin orbitals, so an $L$-site chain uses $2L$ qubits. The initial state is a Fock occupation string with one character per site:\n",
        "\n",
        "$$\n",
        "\\begin{aligned}\n",
        "\\texttt{empty (e)} &\\longleftrightarrow \\left(n_{i,\\uparrow}, n_{i,\\downarrow}\\right) = (0,0), \\\\\n",
        "\\texttt{up (u)} &\\longleftrightarrow \\left(n_{i,\\uparrow}, n_{i,\\downarrow}\\right) = (1,0), \\\\\n",
        "\\texttt{down (d)} &\\longleftrightarrow \\left(n_{i,\\uparrow}, n_{i,\\downarrow}\\right) = (0,1), \\\\\n",
        "\\texttt{both (b)} &\\longleftrightarrow \\left(n_{i,\\uparrow}, n_{i,\\downarrow}\\right) = (1,1).\n",
        "\\end{aligned}\n",
        "$$\n",
        "\n",
        "The string is read least significant site first, so the rightmost character is site 0 and the physical site index of character $k$ counted from the left is $i = L - 1 - k$. The string `udud` therefore means site 0 is `d`, site 1 is `u`, site 2 is `d`, and site 3 is `u`.\n",
        "\n",
        "Next, prepare a half-filled Néel state, alternating `u` and `d`, with the central site replaced by a vacancy. The vacancy is the disturbance whose spreading is tracked.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "id": "bc518931",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Vacancy at site 15\n",
            "udududududududueudududududududu\n"
          ]
        }
      ],
      "source": [
        "vacancy_site = site_count // 2\n",
        "\n",
        "# Build in increasing site order, then reverse: the leftmost character of the\n",
        "# string is the highest site index.\n",
        "site_occupations = [\n",
        "    \"u\" if site % 2 == 0 else \"d\" for site in range(site_count)\n",
        "]\n",
        "site_occupations[vacancy_site] = \"e\"\n",
        "initial_state = \"\".join(reversed(site_occupations))\n",
        "\n",
        "print(f\"Vacancy at site {vacancy_site}\")\n",
        "print(initial_state)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "3d87b6ec",
      "metadata": {},
      "source": [
        "### Build the observables\n",
        "\n",
        "Each observable is a dictionary that represents a weighted sum of products of occupation operators:\n",
        "\n",
        "$$\n",
        "O = \\sum_{\\alpha} c_{\\alpha} \\prod_{(\\sigma, j) \\in T_{\\alpha}} n_{j,\\sigma}.\n",
        "$$\n",
        "\n",
        "Every entry of `operators` is one term $T_{\\alpha}$, written as a list of `[name, site]` pairs that multiply together, and the matching entry of `coefficients` is its weight. The two names are `\"Nup\"` for $n_{j,\\uparrow}$ and `\"Ndn\"` for $n_{j,\\downarrow}$. A single one-point observable $n_{3,\\uparrow}$ is written as follows:\n",
        "\n",
        "```python\n",
        "{\"operators\": [[[\"Nup\", 3]]], \"coefficients\": [1.0]}\n",
        "```\n",
        "\n",
        "The next cell prints the smallest complete set of inputs the function accepts, which is useful as a shape reference while reading the code that follows.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "f996ddab",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "model: {'hopping_amplitude': 1.0, 'on_site_interaction': [2.0, 2.5], 'chemical_potential': 0.5}\n",
            "simulation: {'step_count': 5, 'time_step': 0.1}\n",
            "observables: [{'operators': [[['Nup', 0]]], 'coefficients': [1.0]}]\n",
            "initial_state: ud\n",
            "shots: 4096\n"
          ]
        }
      ],
      "source": [
        "minimal_inputs = {\n",
        "    \"model\": {\n",
        "        \"hopping_amplitude\": 1.0,\n",
        "        \"on_site_interaction\": [2.0, 2.5],\n",
        "        \"chemical_potential\": 0.5,\n",
        "    },\n",
        "    \"simulation\": {\"step_count\": 5, \"time_step\": 0.1},\n",
        "    \"observables\": [{\"operators\": [[[\"Nup\", 0]]], \"coefficients\": [1.0]}],\n",
        "    \"initial_state\": \"ud\",\n",
        "    \"shots\": 4096,\n",
        "}\n",
        "\n",
        "for key, value in minimal_inputs.items():\n",
        "    print(f\"{key}: {value}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "3365e2ae",
      "metadata": {},
      "source": [
        "Three groups of observables are needed. The charge density $n_i = n_{i\\uparrow} + n_{i\\downarrow}$ and the local spin $S^z_i = (n_{i\\uparrow} - n_{i\\downarrow})/2$ are both two-term sums over the same operators, differing only in their coefficients. The spin correlator $S^z_i S^z_{i_*}$ between site $i$ and the vacancy site expands into four two-point terms with coefficients $\\pm 1/4$.\n",
        "\n",
        "Occupation operators commute, so the order of factors within a term is irrelevant, but the order of the `observables` list fixes the order of the returned values. The code records an index for each group while it builds the list. The correlator skips the vacancy site, where the product would repeat a factor.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "id": "cd777495",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "92 observables\n",
            "Example correlator: {'operators': [[['Nup', 0], ['Nup', 15]], [['Nup', 0], ['Ndn', 15]], [['Ndn', 0], ['Nup', 15]], [['Ndn', 0], ['Ndn', 15]]], 'coefficients': [0.25, -0.25, -0.25, 0.25]}\n"
          ]
        }
      ],
      "source": [
        "observables = []\n",
        "charge_indices = []\n",
        "spin_indices = []\n",
        "correlator_indices = []\n",
        "\n",
        "# Charge density: n_i = n_{i,up} + n_{i,down}\n",
        "for site in range(site_count):\n",
        "    charge_indices.append(len(observables))\n",
        "    observables.append(\n",
        "        {\n",
        "            \"operators\": [[[\"Nup\", site]], [[\"Ndn\", site]]],\n",
        "            \"coefficients\": [1.0, 1.0],\n",
        "        }\n",
        "    )\n",
        "\n",
        "# Local spin: S^z_i = (n_{i,up} - n_{i,down}) / 2\n",
        "for site in range(site_count):\n",
        "    spin_indices.append(len(observables))\n",
        "    observables.append(\n",
        "        {\n",
        "            \"operators\": [[[\"Nup\", site]], [[\"Ndn\", site]]],\n",
        "            \"coefficients\": [0.5, -0.5],\n",
        "        }\n",
        "    )\n",
        "\n",
        "# Spin correlator: S^z_i S^z_{i_*}, expanded into four two-point terms.\n",
        "correlator_sites = [\n",
        "    site for site in range(site_count) if site != vacancy_site\n",
        "]\n",
        "for site in correlator_sites:\n",
        "    correlator_indices.append(len(observables))\n",
        "    observables.append(\n",
        "        {\n",
        "            \"operators\": [\n",
        "                [[\"Nup\", site], [\"Nup\", vacancy_site]],\n",
        "                [[\"Nup\", site], [\"Ndn\", vacancy_site]],\n",
        "                [[\"Ndn\", site], [\"Nup\", vacancy_site]],\n",
        "                [[\"Ndn\", site], [\"Ndn\", vacancy_site]],\n",
        "            ],\n",
        "            \"coefficients\": [0.25, -0.25, -0.25, 0.25],\n",
        "        }\n",
        "    )\n",
        "\n",
        "print(f\"{len(observables)} observables\")\n",
        "print(\"Example correlator:\", observables[correlator_indices[0]])"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "c2f5e708",
      "metadata": {},
      "source": [
        "## Step 2: Optimize problem for quantum hardware execution\n",
        "\n",
        "This tutorial has no transpilation step. The function performs the following steps internally [\\[2\\]](#references), and you cannot configure them from the notebook:\n",
        "\n",
        "* Maps the fermionic modes to qubits under an ordering chosen to keep Pauli weights low on a heavy-hex lattice\n",
        "* Uses a fermionic SWAP network so every hopping term becomes a two-qubit rotation\n",
        "* Embeds the chain along a path that avoids weak qubits\n",
        "* Applies dynamical decoupling and Pauli twirling during compilation\n",
        "\n",
        "You choose the backend and the shot count. The function takes a backend name as a string, so this tutorial selects the least busy device large enough for the chain and passes its name.\n",
        "\n",
        "The smallest signal that you need to resolve sets the shot count. The charge tracer reaches amplitudes near $0.15$, but the connected spin tracer peaks around $0.04$, roughly four times smaller. The published study used 20,000 shots per circuit for this reason [\\[2\\]](#references).\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "id": "df8de555",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "ibm_boston, 20000 shots per circuit\n"
          ]
        }
      ],
      "source": [
        "service = QiskitRuntimeService()\n",
        "backend = service.backend(\"ibm_boston\")\n",
        "backend_name = backend.name\n",
        "shot_count = 20_000\n",
        "\n",
        "print(f\"{backend_name}, {shot_count} shots per circuit\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "6b49a90c",
      "metadata": {},
      "source": [
        "## Step 3: Execute by using Qiskit primitives\n",
        "\n",
        "Next, submit the simulation. `run` returns a job handle immediately rather than blocking.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "85c07067",
      "metadata": {},
      "outputs": [],
      "source": [
        "job = sim_qf.run(\n",
        "    model=model,\n",
        "    simulation=simulation,\n",
        "    observables=observables,\n",
        "    initial_state=initial_state,\n",
        "    shots=shot_count,\n",
        "    backend_name=backend_name,\n",
        ")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "cd448379",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "6261b073-be19-48a1-be60-2f369b25613c QUEUED\n"
          ]
        }
      ],
      "source": [
        "print(job.job_id, job.status())"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "36ed5122",
      "metadata": {},
      "source": [
        "Use `job.status()` to follow progress and `job.result()` to block until the job finishes. The result is a mapping whose `observable_expectation_values` entry is indexed first by observable, in the order they were supplied, and then by time. There are `step_count + 1` time points because the value at $t=0$ is included.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "0b54e632",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Result keys: ['observable_expectation_values', 'observables', 'provider_job_ids', 'warnings']\n"
          ]
        }
      ],
      "source": [
        "result = job.result()\n",
        "\n",
        "print(\"Result keys:\", sorted(result.keys()))"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "5c936d75",
      "metadata": {},
      "source": [
        "## Step 4: Post-process and return result in desired classical format\n",
        "\n",
        "### Reconstruct the correlators\n",
        "\n",
        "Split the flat array back into the three groups using the indices recorded earlier, then form the two tracer correlators. The charge tracer is the change in per-site density relative to $t=0$:\n",
        "\n",
        "$$\n",
        "C^c_i(t) = \\langle n_{i\\uparrow}(t) + n_{i\\downarrow}(t) \\rangle - \\langle n_{i\\uparrow}(0) + n_{i\\downarrow}(0) \\rangle.\n",
        "$$\n",
        "\n",
        "The spin tracer is the connected correlator between each site and the vacancy site $i_*$:\n",
        "\n",
        "$$\n",
        "C^s_i(t) = 4 \\left( \\langle S^z_i S^z_{i_*} \\rangle - \\langle S^z_i \\rangle \\langle S^z_{i_*} \\rangle \\right).\n",
        "$$\n",
        "\n",
        "Both are built from Fock-basis operators, so they come from the same measured bitstrings and the second one costs no extra hardware execution.\n",
        "\n",
        "When you constructed the observables, the code recorded the position of each charge (`charge_indices`), spin (`spin_indices`), and spin-correlation (`correlator_indices`) operator in the list. Because the simulator returns the expectation values in the same order, you can use these positions to extract the corresponding results. For the spin correlations, the code also uses `correlator_sites` to place each result at the correct lattice position, and leaves the vacancy as `NaN`.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 13,
      "id": "6d6bc02e",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "92 observables, 31 time points up to t = 4.50\n",
            "Peak amplitudes: charge 0.978, spin 0.168\n"
          ]
        }
      ],
      "source": [
        "expectation_values = np.asarray(\n",
        "    result[\"observable_expectation_values\"], dtype=float\n",
        ")\n",
        "times = np.arange(expectation_values.shape[1]) * time_step\n",
        "\n",
        "charge_density = expectation_values[charge_indices]\n",
        "local_spin = expectation_values[spin_indices]\n",
        "spin_correlator = np.full_like(charge_density, np.nan)\n",
        "spin_correlator[correlator_sites] = expectation_values[correlator_indices]\n",
        "\n",
        "# Charge tracer: change in per-site density relative to t = 0.\n",
        "charge_tracer = charge_density - charge_density[:, [0]]\n",
        "\n",
        "# Spin tracer: connected correlator with the vacancy site.\n",
        "spin_tracer = 4.0 * (spin_correlator - local_spin * local_spin[vacancy_site])\n",
        "spin_tracer[vacancy_site] = np.nan\n",
        "\n",
        "print(\n",
        "    f\"{expectation_values.shape[0]} observables, \"\n",
        "    f\"{expectation_values.shape[1]} time points up to t = {times[-1]:.2f}\"\n",
        ")\n",
        "print(\n",
        "    f\"Peak amplitudes: charge {np.nanmax(np.abs(charge_tracer)):.3f}, \"\n",
        "    f\"spin {np.nanmax(np.abs(spin_tracer)):.3f}\"\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "9195ba20",
      "metadata": {},
      "source": [
        "### Visualize the wavefronts\n",
        "\n",
        "If you plot each correlator as a heatmap of space and time, the light cone is directly visible. Time flows downward, and the horizontal axis is the position along the chain.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 34,
      "id": "c75549b7",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/fermi-hubbard-spin-charge-dynamics/extracted-outputs/c75549b7-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "def plot_tracer(axis, tracer, title):\n",
        "    \"\"\"Draw one tracer correlator as a space-time heatmap, time flowing down.\"\"\"\n",
        "    # Scale from the off-centre sites. The vacancy site carries an order-one\n",
        "    # change in charge density and is left to clip rather than set the scale.\n",
        "    off_center = np.ones(site_count, dtype=bool)\n",
        "    off_center[vacancy_site] = False\n",
        "    limit = np.nanpercentile(np.abs(tracer[off_center]), 98)\n",
        "\n",
        "    colormap = qv.QCTRL_DIVERGENT_COLORMAP.reversed().copy()\n",
        "    colormap.set_bad(\"0.82\")  # the spin tracer has no value at the vacancy\n",
        "\n",
        "    image = axis.imshow(\n",
        "        tracer.T,\n",
        "        aspect=\"auto\",\n",
        "        origin=\"upper\",\n",
        "        cmap=colormap,\n",
        "        vmin=-limit,\n",
        "        vmax=limit,\n",
        "        extent=(-0.5, site_count - 0.5, times[-1], times[0]),\n",
        "        interpolation=\"nearest\",\n",
        "    )\n",
        "\n",
        "    axis.annotate(\n",
        "        \"vacancy\",\n",
        "        xy=(vacancy_site, times[-1]),\n",
        "        xytext=(vacancy_site, times[-1] + 0.55),\n",
        "        ha=\"center\",\n",
        "        va=\"top\",\n",
        "        fontsize=8,\n",
        "        annotation_clip=False,\n",
        "        arrowprops=dict(arrowstyle=\"-\", lw=0.8, color=\"0.4\"),\n",
        "    )\n",
        "\n",
        "    axis.set_title(title)\n",
        "    axis.set_xlabel(\"Site\")\n",
        "    return image\n",
        "\n",
        "\n",
        "fig, axes = plt.subplots(1, 2, figsize=(11, 4.5), sharey=True)\n",
        "\n",
        "charge_image = plot_tracer(\n",
        "    axes[0], charge_tracer, r\"Charge tracer $C^{c}_{i}(t)$\"\n",
        ")\n",
        "spin_image = plot_tracer(axes[1], spin_tracer, r\"Spin tracer $C^{s}_{i}(t)$\")\n",
        "\n",
        "axes[0].set_ylabel(r\"Time $t$  $(1/t_h)$\")\n",
        "fig.colorbar(charge_image, ax=axes[0], fraction=0.046)\n",
        "fig.colorbar(spin_image, ax=axes[1], fraction=0.046)\n",
        "fig.suptitle(f\"$L$ = {site_count},  $U/t_h$ = {on_site_interaction:.0f}\")\n",
        "fig.tight_layout()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "bcef6a44",
      "metadata": {},
      "source": [
        "### Extract the wavefront velocities\n",
        "\n",
        "Turn the heatmaps into numbers by following the detection protocol from the published study [\\[2\\]](#references). At each time slice, symmetrize the spatial profile about the vacancy, which averages away the antisymmetric part of the hardware noise, normalize by the largest off-center amplitude, and record the outermost site above a threshold. Detections are then cut into runs of consecutive time steps, and only the longest run is fit, which keeps the early-time transient out of the result.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "72411831",
      "metadata": {},
      "outputs": [],
      "source": [
        "def front_distance(profile, center, threshold):\n",
        "    \"\"\"Distance from `center` to the outermost site above a relative threshold.\"\"\"\n",
        "    values = np.nan_to_num(np.asarray(profile, dtype=float))\n",
        "    length = values.size\n",
        "\n",
        "    symmetrized = values.copy()\n",
        "    for site in range(length):\n",
        "        mirror = 2 * center - site\n",
        "        if 0 <= mirror < length:\n",
        "            symmetrized[site] = 0.5 * (values[site] + values[mirror])\n",
        "\n",
        "    amplitude = np.abs(symmetrized)\n",
        "    off_center = np.ones(length, dtype=bool)\n",
        "    off_center[center] = False\n",
        "\n",
        "    peak = amplitude[off_center].max()\n",
        "    if peak <= 0.0:\n",
        "        return np.nan\n",
        "    above = np.flatnonzero((amplitude / peak > threshold) & off_center)\n",
        "    return float(np.abs(above - center).max()) if above.size else np.nan\n",
        "\n",
        "\n",
        "def longest_run(mask):\n",
        "    \"\"\"Indices of the longest run of consecutive True values.\"\"\"\n",
        "    best, current = [], []\n",
        "    for index, keep in enumerate(mask):\n",
        "        current = current + [index] if keep else []\n",
        "        if len(current) > len(best):\n",
        "            best = current\n",
        "    return np.array(best, dtype=int)\n",
        "\n",
        "\n",
        "def front_velocity(tracer, center, threshold=0.3, quantile=0.75, min_block=5):\n",
        "    \"\"\"Fit front position against time, using only the longest contiguous run.\"\"\"\n",
        "    distances = np.array(\n",
        "        [\n",
        "            front_distance(tracer[:, k], center, threshold)\n",
        "            for k in range(tracer.shape[1])\n",
        "        ]\n",
        "    )\n",
        "    edge = max(center, tracer.shape[0] - 1 - center)\n",
        "    usable = (\n",
        "        np.isfinite(distances) & (distances >= 2) & (distances <= edge - 2)\n",
        "    )\n",
        "\n",
        "    block = longest_run(usable)\n",
        "    if block.size < min_block:\n",
        "        return np.nan, distances, np.array([], dtype=int), np.nan\n",
        "\n",
        "    candidates = [\n",
        "        stats.theilslopes(distances[block[:k]], times[block[:k]]).slope\n",
        "        for k in range(min_block, block.size + 1)\n",
        "    ]\n",
        "    cutoff = np.quantile(candidates, quantile)\n",
        "    keep = max(i for i, value in enumerate(candidates) if value >= cutoff)\n",
        "\n",
        "    fit_block = block[: keep + min_block]\n",
        "    slope, intercept = np.polyfit(times[fit_block], distances[fit_block], 1)\n",
        "    return float(slope), distances, fit_block, float(intercept)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "e4ec745a",
      "metadata": {},
      "source": [
        "The fitted velocity depends on the detection threshold and the trimming quantile, and that sensitivity dominates the shot noise. Scan both and report the median with the full envelope, as the published analysis does.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "d127f869",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "v_c = 1.93  [1.81, 2.18]\n",
            "v_s = 1.61  [1.23, 2.12]\n",
            "v_s / v_c = 0.83\n"
          ]
        }
      ],
      "source": [
        "thresholds = [0.15, 0.225, 0.30, 0.375, 0.45]\n",
        "quantiles = [0.60, 0.675, 0.75, 0.825, 0.90]\n",
        "\n",
        "\n",
        "def velocity_envelope(tracer, center):\n",
        "    \"\"\"Median, minimum, and maximum velocity over the parameter scan.\"\"\"\n",
        "    values = [\n",
        "        front_velocity(tracer, center, threshold=p, quantile=q)[0]\n",
        "        for p in thresholds\n",
        "        for q in quantiles\n",
        "    ]\n",
        "    values = np.array([v for v in values if np.isfinite(v)])\n",
        "    if values.size == 0:\n",
        "        return np.nan, np.nan, np.nan\n",
        "    return float(np.median(values)), float(values.min()), float(values.max())\n",
        "\n",
        "\n",
        "charge_velocity, charge_low, charge_high = velocity_envelope(\n",
        "    charge_tracer, vacancy_site\n",
        ")\n",
        "spin_velocity, spin_low, spin_high = velocity_envelope(\n",
        "    spin_tracer, vacancy_site\n",
        ")\n",
        "\n",
        "_, charge_front, charge_block, charge_intercept = front_velocity(\n",
        "    charge_tracer, vacancy_site\n",
        ")\n",
        "_, spin_front, spin_block, spin_intercept = front_velocity(\n",
        "    spin_tracer, vacancy_site\n",
        ")\n",
        "\n",
        "print(f\"v_c = {charge_velocity:.2f}  [{charge_low:.2f}, {charge_high:.2f}]\")\n",
        "print(f\"v_s = {spin_velocity:.2f}  [{spin_low:.2f}, {spin_high:.2f}]\")\n",
        "print(f\"v_s / v_c = {spin_velocity / charge_velocity:.2f}\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 31,
      "id": "640db359",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/fermi-hubbard-spin-charge-dynamics/extracted-outputs/640db359-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "PUBLISHED = {\"charge\": 1.78, \"spin\": 1.46}\n",
        "\n",
        "edge = max(vacancy_site, site_count - 1 - vacancy_site)\n",
        "fig, axis = plt.subplots(figsize=(7, 4.6))\n",
        "\n",
        "axis.axhspan(0, 2, color=\"0.92\", zorder=0)\n",
        "axis.axhspan(edge - 2, edge + 2, color=\"0.92\", zorder=0)\n",
        "\n",
        "series = [\n",
        "    (\n",
        "        \"charge\",\n",
        "        \"$v_c$\",\n",
        "        charge_front,\n",
        "        charge_block,\n",
        "        charge_velocity,\n",
        "        charge_intercept,\n",
        "        charge_low,\n",
        "        charge_high,\n",
        "        \"C0\",\n",
        "        \"o\",\n",
        "    ),\n",
        "    (\n",
        "        \"spin\",\n",
        "        \"$v_s$\",\n",
        "        spin_front,\n",
        "        spin_block,\n",
        "        spin_velocity,\n",
        "        spin_intercept,\n",
        "        spin_low,\n",
        "        spin_high,\n",
        "        \"C1\",\n",
        "        \"s\",\n",
        "    ),\n",
        "]\n",
        "\n",
        "for name, sym, front, block, v, b, lo, hi, color, marker in series:\n",
        "    excluded = np.setdiff1d(np.arange(front.size), block)\n",
        "    axis.plot(\n",
        "        times[excluded],\n",
        "        front[excluded],\n",
        "        marker,\n",
        "        color=color,\n",
        "        markersize=5,\n",
        "        markerfacecolor=\"none\",\n",
        "        alpha=0.45,\n",
        "        zorder=2,\n",
        "    )\n",
        "    axis.plot(\n",
        "        times[block],\n",
        "        front[block],\n",
        "        marker,\n",
        "        color=color,\n",
        "        markersize=5,\n",
        "        label=f\"{name}, {sym} = {v:.2f} [{lo:.2f}, {hi:.2f}]\",\n",
        "        zorder=3,\n",
        "    )\n",
        "\n",
        "    span = times[block]\n",
        "    anchor = v * span[0] + b\n",
        "    axis.plot(\n",
        "        span,\n",
        "        v * (span - span[0]) + anchor,\n",
        "        \"-\",\n",
        "        color=color,\n",
        "        lw=1.6,\n",
        "        zorder=4,\n",
        "    )\n",
        "    axis.fill_between(\n",
        "        span,\n",
        "        lo * (span - span[0]) + anchor,\n",
        "        hi * (span - span[0]) + anchor,\n",
        "        color=color,\n",
        "        alpha=0.15,\n",
        "        zorder=1,\n",
        "    )\n",
        "    axis.plot(\n",
        "        span,\n",
        "        PUBLISHED[name] * (span - span[0]) + anchor,\n",
        "        \"--\",\n",
        "        color=color,\n",
        "        lw=1.2,\n",
        "        alpha=0.8,\n",
        "        zorder=4,\n",
        "    )\n",
        "\n",
        "axis.plot([], [], \"--\", color=\"0.4\", lw=1.2, label=\"published [2]\")\n",
        "axis.plot(\n",
        "    [],\n",
        "    [],\n",
        "    \"o\",\n",
        "    color=\"0.5\",\n",
        "    markerfacecolor=\"none\",\n",
        "    label=\"excluded from fit\",\n",
        ")\n",
        "\n",
        "axis.set_xlabel(r\"Time $t$  $(1/t_h)$\")\n",
        "axis.set_ylabel(\"Front distance from vacancy (sites)\")\n",
        "axis.set_xlim(0, times[-1])\n",
        "axis.set_ylim(\n",
        "    0, min(edge + 1, np.nanmax(np.r_[charge_front, spin_front]) + 3)\n",
        ")\n",
        "axis.legend(loc=\"upper left\", frameon=False, fontsize=9)\n",
        "axis.set_title(\n",
        "    f\"L = {site_count}, $U/t_h$ = {on_site_interaction:.0f}, front propagation\",\n",
        "    fontsize=11,\n",
        ")\n",
        "\n",
        "fig.tight_layout()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "acd487cb",
      "metadata": {},
      "source": [
        "Drawing the detections back onto the heatmaps shows whether the algorithm tracked\n",
        "the leading edge. Filled markers are the detections that entered\n",
        "the fit. Open markers were cut by the masking rules, either because they sit too\n",
        "close to the vacancy, too close to the chain end, or were late enough that the front has\n",
        "broadened and its position is no longer sharp. The dashed line is the fit over the\n",
        "times it used, and its slope is the velocity in each panel title.\n",
        "\n",
        "The vacancy column is marked below the axis. The charge clips there because its value\n",
        "is of order one, and spin has no value there, so it is drawn in gray.\n",
        "\n",
        "The charge cone opens wider than the spin cone. Both disturbances leave the same\n",
        "defect at the same moment and travel outward at different speeds.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 35,
      "id": "25192e92",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/fermi-hubbard-spin-charge-dynamics/extracted-outputs/25192e92-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "def plot_tracer_front(axis, tracer, front, block, velocity, intercept, title):\n",
        "    \"\"\"Overlay the detected wavefront and the fitted cone on a tracer heatmap.\"\"\"\n",
        "    image = plot_tracer(axis, tracer, title)\n",
        "\n",
        "    excluded = np.setdiff1d(np.arange(front.size), block)\n",
        "    for index, filled in ((block, True), (excluded, False)):\n",
        "        good = index[np.isfinite(front[index])]\n",
        "        if good.size == 0:\n",
        "            continue\n",
        "        for sign in (-1, 1):\n",
        "            axis.plot(\n",
        "                vacancy_site + sign * front[good],\n",
        "                times[good],\n",
        "                \"o\",\n",
        "                markersize=4,\n",
        "                color=\"k\",\n",
        "                markerfacecolor=\"k\" if filled else \"none\",\n",
        "                markeredgewidth=0.9,\n",
        "                alpha=1.0 if filled else 0.4,\n",
        "            )\n",
        "\n",
        "    span = times[block]\n",
        "    for sign in (-1, 1):\n",
        "        axis.plot(\n",
        "            vacancy_site + sign * (velocity * span + intercept),\n",
        "            span,\n",
        "            \"--\",\n",
        "            color=\"k\",\n",
        "            lw=1.5,\n",
        "        )\n",
        "\n",
        "    axis.set_title(f\"{title}    $v$ = {velocity:.2f}\")\n",
        "    return image\n",
        "\n",
        "\n",
        "fig, axes = plt.subplots(1, 2, figsize=(11, 4.5), sharey=True)\n",
        "\n",
        "charge_image = plot_tracer_front(\n",
        "    axes[0],\n",
        "    charge_tracer,\n",
        "    charge_front,\n",
        "    charge_block,\n",
        "    charge_velocity,\n",
        "    charge_intercept,\n",
        "    r\"Charge tracer $C^{c}_{i}(t)$\",\n",
        ")\n",
        "spin_image = plot_tracer_front(\n",
        "    axes[1],\n",
        "    spin_tracer,\n",
        "    spin_front,\n",
        "    spin_block,\n",
        "    spin_velocity,\n",
        "    spin_intercept,\n",
        "    r\"Spin tracer $C^{s}_{i}(t)$\",\n",
        ")\n",
        "\n",
        "axes[0].set_ylabel(r\"Time $t$  $(1/t_h)$\")\n",
        "fig.colorbar(charge_image, ax=axes[0], fraction=0.046)\n",
        "fig.colorbar(spin_image, ax=axes[1], fraction=0.046)\n",
        "fig.suptitle(f\"$L$ = {site_count},  $U/t_h$ = {on_site_interaction:.0f}\")\n",
        "fig.tight_layout()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "913810be",
      "metadata": {},
      "source": [
        "The charge front outruns the spin front. A hole in a one-dimensional Mott background does not travel as an electron. Instead, it fractionalizes into a holon carrying the charge and a spinon carrying the spin disturbance, and at $U/t_h = 8$ the spinon is the slower of the two.\n",
        "\n",
        "The ordering depends on the interaction. At $U = 0$, the model is free fermions with $v_c = 2$ and $v_s = 4$, so spin is faster. Repulsion suppresses spin transport while leaving charge transport nearly untouched, and by $U/t_h = 8$ the published hardware values are $v_c = 1.78$ and $v_s = 1.46$ [\\[2\\]](#references). Rerunning with several values of $U$ and plotting $v_s / v_c$ maps the crossover.\n",
        "\n",
        "Note these caveats on the numbers. Threshold-based detection is biased slightly low, so a charge velocity a few percent under the published value is expected. The spin envelope is wide because the connected correlator sits closer to the noise floor, which is why the median of the scan rather than a single fit is quoted.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "84f72b94",
      "metadata": {},
      "source": [
        "## References\n",
        "\n",
        "\\[1] Essler, F. H. L., Frahm, H., Gohmann, F., Klumper, A., & Korepin, V. E. (2005). *The One-Dimensional Hubbard Model*. Cambridge University Press.\n",
        "\n",
        "\\[2] Hartnett, G. S., Najafi, K. S., Khindanov, A., Liao, H., Schutzman, M., Hush, M. R., Biercuk, M. J., & Baum, Y. (2026). Fast, accurate, high-resolution simulation of large-scale Fermi-Hubbard models on a digital quantum processor. *arXiv:2605.04025*.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "id": "a1b8767d",
      "source": "© IBM Corp., 2017-2026"
    }
  ],
  "metadata": {
    "hours": 1,
    "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": 180
  },
  "nbformat": 4,
  "nbformat_minor": 5
}