{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "df27beb5-e42d-4fa8-a47b-cbd4b7fd61a4",
   "metadata": {},
   "source": [
    "---\n",
    "title: Simulate a 127-qubit kicked-Ising model\n",
    "description: Use the pauli-prop Qiskit addon to simulate time dynamics of a kicked-Ising model\n",
    "---\n",
    "\n",
    "\n",
    "# Simulate a 127-qubit kicked-Ising model\n",
    "\n",
    "In this guide we use the ``pauli-prop`` package to classically simulate the time dynamics of a 127-qubit kicked-Ising model on a heavy-hex lattice and approximate the single-site magnetization, $\\langle Z_{62} \\rangle$. The Hamiltonian considered is:\n",
    "\n",
    "$H = -J\\sum\\limits_{\\langle i,j \\rangle} Z_iZ_j + h\\sum\\limits_iX_i$\n",
    "\n",
    "where $J>0$ describes the coupling of nearest-neighbor spins, $i<j$, and $h$ is the global transverse field. A first-order Trotter decomposition of the time-evolved operator will be implemented as a quantum circuit, $U$, over $20$ Trotter steps. The coupling constant, $J$, will be fixed at $J=-\\frac{\\pi}{2}$ such that $U$ is Clifford any time $h\\mod\\frac{\\pi}{2}=0$.\n",
    "\n",
    "Workflow:\n",
    "\n",
    "- Create quantum circuits implementing the Trotterized Ising model\n",
    "    - Vary $h$ between Clifford points, $0.0$ and $\\frac{\\pi}{2}$\n",
    "- Choose a site in the middle of the lattice (qubit $62$) and estimate the magnetization after $20$ Trotter steps for each circuit\n",
    "- Observe the magnetization roll off from $1.0$ to $0.0$ as $h$ moves from one Clifford point to another"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab545dd3-7771-4226-921a-e797ad4749c1",
   "metadata": {},
   "source": [
    "## Create circuits\n",
    "\n",
    "First, we use ``FakeSherbrooke`` from ``qiskit-ibm-runtime`` to get the edge indices for the coupling map of the 127-qubit Eagle QPU. Once we have the edge indices, we create the Trotter circuits for the varying values of $h$ we want to simulate. The circuits will implement a first-order Trotterization of the Hamiltonian using 20 Trotter steps and 5420 gates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "92612eff-e570-45a0-be1c-16d22dfa9690",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-08-20T23:46:48.371100Z",
     "iopub.status.busy": "2026-08-20T23:46:48.371006Z",
     "iopub.status.idle": "2026-08-20T23:46:49.865257Z",
     "shell.execute_reply": "2026-08-20T23:46:49.864711Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from qiskit import QuantumCircuit\n",
    "from qiskit_ibm_runtime.fake_provider import FakeSherbrooke\n",
    "\n",
    "# Get the edges in the device coupling map\n",
    "backend = FakeSherbrooke()\n",
    "edges = backend.coupling_map.get_edges()\n",
    "\n",
    "# Num Trotter steps\n",
    "num_steps = 20\n",
    "\n",
    "# Create Trotter circuits with varying global fields\n",
    "hs = [i * np.pi / 32 for i in range(17)]\n",
    "circuits = []\n",
    "for h in hs:\n",
    "    circuit = QuantumCircuit(backend.num_qubits)\n",
    "    for _ in range(num_steps):\n",
    "        circuit.rx(h, [i for i in range(backend.num_qubits)])\n",
    "        for edge in edges:\n",
    "            circuit.rzz(-np.pi / 2, edge[0], edge[1])\n",
    "    circuits.append(circuit)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2dfebe6e-3956-40eb-b19c-97ea1c20516c",
   "metadata": {},
   "source": [
    "## Simulate the time dynamics with Pauli propagation\n",
    "\n",
    "Next, we specify the observable, $O$, we want to measure. In this case, we estimate $\\langle Z_{62} \\rangle$.\n",
    "\n",
    "To study the magnetization between Clifford points, we approximately propagate $O$ to the beginning of the circuit, $U$, resulting in a new observable, $O^\\prime \\approx U^{\\dagger}OU$.\n",
    "\n",
    "During propagation, the number of terms in the evolved observable grows exponentially in the number of gates, so it is usually necessary to limit how large it can grow. Here, we limit the number of terms in $O^\\prime$ to $10^5$ and truncate terms with coefficients smaller than $10^{-5}$. Truncating terms results in some error.\n",
    "Once we have $O^\\prime$, we can easily estimate the expectation value, $\\langle0|U^{\\dagger}OU|0\\rangle = \\langle0|O^\\prime|0\\rangle$, by summing the coefficients of the diagonal terms in $O^\\prime$ (terms consisting only of ``I`` or ``Z`` on all qubits).\n",
    "\n",
    "We set ``frame='h'`` to specify we want to evolve the observable in the Heisenberg framework ($U^{\\dagger}OU$). Use ``frame=s`` for Schrödinger evolution ($UOU^{\\dagger}$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f0cfea45-2e5a-4e9e-9dd4-4b5502128e3a",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-08-20T23:46:49.866789Z",
     "iopub.status.busy": "2026-08-20T23:46:49.866636Z",
     "iopub.status.idle": "2026-08-20T23:49:49.909638Z",
     "shell.execute_reply": "2026-08-20T23:49:49.909193Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Finshed 17 simulations in 180s\n"
     ]
    }
   ],
   "source": [
    "import time\n",
    "\n",
    "from pauli_prop import propagate_through_circuit\n",
    "from qiskit.quantum_info import Pauli, SparsePauliOp\n",
    "\n",
    "# Z_62 observable\n",
    "id_pauli = Pauli(\"I\" * backend.num_qubits)\n",
    "observable = SparsePauliOp(id_pauli.dot(Pauli(\"Z\"), [62]))\n",
    "\n",
    "# Estimate the single-site magnetization for varying magnetic fields\n",
    "times = []\n",
    "approx_evs = []\n",
    "for circuit in circuits:\n",
    "    st = time.perf_counter()\n",
    "    propagated_obs = propagate_through_circuit(observable, circuit, 100_000, 1e-5, frame=\"h\")[0]\n",
    "    times.append(time.perf_counter() - st)\n",
    "    approx_evs.append(float(propagated_obs.coeffs[~propagated_obs.paulis.x.any(axis=1)].sum()))\n",
    "print(f\"Finshed {len(circuits)} simulations in {sum(times):.0f}s\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d214da4c-d5cf-43f2-97a7-380858d4a953",
   "metadata": {},
   "source": [
    "## Visualize\n",
    "\n",
    "Finally, we plot the magnetization as a function of $h$. When $h=0$, the circuit does not affect the initial state in the computational basis, and we know $\\langle0|O|0\\rangle = 1.0$. As $h$ moves further from $0$, the magnetization decays monotonically until $h$ reaches the next Clifford angle, $\\frac{\\pi}{2}$, and the magnetization converges to $0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "01a3e294-f5c3-485c-b027-a14d9a6384c9",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-08-20T23:49:49.911044Z",
     "iopub.status.busy": "2026-08-20T23:49:49.910974Z",
     "iopub.status.idle": "2026-08-20T23:49:50.161262Z",
     "shell.execute_reply": "2026-08-20T23:49:50.160903Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1190fcdd0>]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjoAAAG0CAYAAAA7Go31AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAARDFJREFUeJzt3Qd8FGX+x/HfbjqQBEJIIYQOUkKRHsAOoiIHKCdWlBM9OfAUjjtARcSG/fyfIhycCp4gKAqicFFEUYEAUqV3pKUQIIWE9P2/ngcSE0ggZZOZnf28X69ld2ZnM8/Mkt1vnjY2h8PhEAAAAAuyG10AAACAqkLQAQAAlkXQAQAAlkXQAQAAlkXQAQAAlkXQAQAAlkXQAQAAluUpbi4/P19OnDgh/v7+YrPZjC4OAAAoAzUNYFpamtSvX1/s9tLrbdw+6KiQExkZWZZzCgAATObo0aPSoEGDUp93+6CjanIKTlRAQEA1vjUAAKCiUlNTdUVFwfd4adw+6BQ0V6mQQ9ABAMC1XKnbCZ2RAQCAZRF0AACAZRF0AACAZRF0AACAZRF0AACAZRF0AACAZRF0AACAZRF0AACAZRF0AACAZbn9zMhVIS/fIesPnZbEtEwJ8feVbk2CxMPOBUMBAHDrGp2ffvpJBgwYoK9EqqZ0Xrx48RVfs3LlSunUqZP4+PhI8+bNZfbs2WKkmO1x0vvV7+WeWWvliflb9L1aVuurMljFHjglX245ru/VMgAAMFmNTnp6unTo0EH+9Kc/yR133HHF7Q8dOiT9+/eXxx57TObOnSsrVqyQESNGSHh4uPTr10+qmwozIz/eJBfHjPiUTL1++v2d5JaocKfvc8pXOyUuJbNwXXigr0we0Mbp+wIAwNXYHA6HKf/8VzU6ixYtkkGDBpW6zfjx42Xp0qWyffv2wnV33323JCcnS0xMTJmvfhoYGCgpKSmVuqinqkVRNTdFA0dRquEqLNBXVo2/0WnNWKUFq4KfXhXBCgAAMyjr97epanTKKzY2Vvr06VNsnarJefLJJ0t9TVZWlr4VPVHOoPrklBZyFBVG1PM9pn4ndWv6SA1vD6np4yk1vT2lho+H1PLxlBrenlLLx0Pf1/Qp8nzBtnr5/GNvD7uuySkppTouhB31fN82YfQPAgC4LZcOOvHx8RIaGlpsnVpW4eXcuXPi5+d3yWumTp0qU6ZMcXpZVMfjsjiZlq1vVa0gWKkAFt2sbpXvDwAAM3LpoFMREydOlLFjxxYuq1AUGRlZ6Z+rRleVxfMD20qT4JqSnpUnGdm5kp6VK+nZeZKRlStnL6w7m5UrGdl5F57LlYysPH2ffuG+PI2NE774Va5tUU+iIgKkbf1AaRnqL96epuqDDgBAlXHpoBMWFiYJCQnF1qll1VZXUm2OokZnqZuzqSHkqhOw6njsuEwfnfu6N6pUU5LqUpWZky8/7k2Uxz7edMXtfzuVIf899VvhsmryuirMXwefqIhAiaofqJd9vTzKVQ6G0AMAXIFLB53o6GhZtmxZsXXLly/X66ubCi9qpJPqHKxiTNGwUxBr1POV7YisOmn7eXvovjdXClbB/j7y9G2tZWdcqmw/nqJvqZm5su14ir6JHC0se4uQWtJOBR99C5DW4QG6r1BJGOkFAHAVphp1dfbsWdm/f79+fPXVV8tbb70lN9xwgwQFBUnDhg11s9Px48flo48+KhxeHhUVJaNGjdJD0r///nv561//qkdilXV4ubNGXRkRAgpGXUkpweriUVfqrT525pwOOTr4nDgfgE6nX9pnSOWxZvVq6eDTtn6ADkFt6gfI6v1JjPQCABiurN/fpgo6avI/FWwu9uCDD+qJAB966CE5fPiw3q7oa8aMGSM7d+6UBg0ayKRJk/R2ZeXsoFPdzTqVDVbq7VevLRp81C0x7feRaUWp4yhtQsKqGEIPAIBlgo4RqiLoVLeqCFaJqZmy40RqYe2Penw8+VyZXvvJIz0Y6QUAqFJuMY8OzlOhxtlDyEMCfPXthlYhhevmrftNnlr0++SMlR1qDwBAVWOcMcqsSXAtpw61BwCgqhF0UO4h9FdqFPt530nJys3jzAIADEfQQbmH0CuXCzvvrTwgA95ZJb8eS+bsAgAMRdBBuaiRXGrYuhpdVZSq6Zlxfyd5775OUremt+xNOCuD31sjr8bslswcancAAMZg1JUFRl2ZbaSXmpfnuSU7ZMnWE3q5eUgteX1Ie7m6YR2DSw0AsAqGlzv5RKH8YrbHyzOLt0vS2Sw9AeEj1zSVMX1blvtyEwAAVPT7m6YrVJlbosJk+ZhrZfDVEaLmGPz3Twfltn/9LBt/O81ZBwBUC4IOqlSdmt7yz6EdZdawLhLi7yMHT6bLkBmx8uLXO+VcNn13AABVi6CDatG3TagsH3Od3Nmpgai5uP+z6pCu3fnlMLU7AICqQ9BBtQms4SVv3tVBPnioi4QG+MihpHS569+xMuWrHZKRncs7AQBwOoIOqt2NrULl2zHXyV1dztfufLj6sNz6fz/L2oOneDcAAE5F0IEhAv285LUhHWT28K56Dp7fTmXI3TPXyuQvt0t6FrU7AADnIOjAUNdfFSLfjLlW7ukWqZfnxP4mt/zfT7LmQNIl8/bEHjglX245ru/VMgAAV8KEgcyjYxo/7T0pE7/YJseTz+nl+3s0lAm3tpZV+07KlK92SlzK71dFV7VA6nIUaqZmAID7SS3j9zdBh6BjKmmZOfLK/3bL3HVH9HJQDW85nZF9yXYF19pSl6Mg7ACA+0llwkC4In9fL3lpcDuZN6K7RNT2LTHkKAUNV6qmh2YsAEBp6KMDU+rZPFheHNTustuosKOas9Q1twAAKAlBB6aVmplTpu3UhUUBACgJQQempa6K7sztAADuh6AD0+rWJEiPriroeFwS9bzaDgCAkhB0YFoedpseQq6UFnYe7t1EbwcAQEkIOjA1NXRcDSEPCyzePOXlcT7czPjxgOxPPGtQ6QAAZsc8Osyj4xLUEHI1ukp1PFZ9cq4K85cH3l8nO06k6guEfvrnaGlUt6bRxQQAVBMmDHTyiYL5nE7PlrtnxsrehLMSUdtPPn0sWt8DAKwvlQkDYXVBNb3l4xHdpUlwTX3ZiPtmrZXEVIaaAwB+Rx8duDTVjDV3RHdpUMdPDp/KkPv+s05Onc0yulgAAJMg6MDl1a/tJ/NG9JCwAF/Zl3hWHnh/vaRklG2yQQCAtRF0YAkN69aQuY90l+Ba3rIzLlWGfbheXyAUAODeCDqwjGb1auk+O7VreMnWo8ny8OwNkpGda3SxAAAGIujAUlqFBch//9Rd/H09Zf3h0/LoRxslMyfP6GIBAAxC0IHltGsQKLOHd5Ma3h6yan+S/GXuJsnOzTe6WAAAAxB0YEmdG9WR9x/sKj6edvl+d6I8uWCz5OYRdgDA3RB0YFnRzerKzGFdxNvDLsu2xcvfF/4q+fkOo4sFAKhGBB1Y2nUt68m7916tL/y5aPNxeXrxNnE4CDsA4C4IOrC8m9uGydtDO4q6yPkn64/KlK92EnYAwE0QdOAWBnSoL68N6aAfz15zWF6N2UPYAQA3QNCB2xjSuYG8OChKP57x4wF55/v9RhcJAFDFCDpwK/f3aCTP9G+tH7+1fK/M/OmA0UUCAFQhgg7czohrmsq4m1vqxy8v2y3/jT1sdJEAAFWEoAO3NPrGFjLqhmb68aQvd8invxw1ukgAgCpA0IHbGnfzVfKnXk304/Ff/CpfbjludJEAAE5G0IHbstlsMun21nJf94aiptYZ++lWidkeb3SxAABORNCBuHvYeWFglNzRKULy8h3y+Ceb5Ic9ifpx7IFTupZH3atlAIDr8TS6AIDR7HabvHZne8nKzZelv8bJI3M2SICfl5xOzy7cJjzQVyYPaCO3RIUbWlYAQPlQowOoxO9h17Mnt48IlNx8R7GQo8SnZMrIjzdJzPY4zhcAuBCCDlDwy2CzSWJaZonno6DhSl0+gmYsAHAdBB3ggvWHTkt8alap50OFnbiUTL0dAMA1EHSAC0qrzanodgAA4xF0gAtC/H2duh0AwHgEHeCCbk2C9Ogq22XOiKfdJvX8fThnAOAiCDrABR52mx5CrpQWdtSIrMHvrdZz7QAAzI+gAxSh5smZfn8nCQss3jylanpevaOddGpYW9Iyc+VPs3+R91buF4eaUhkAYFo2h5t/UqempkpgYKCkpKRIQECA0cWBSagh5Gp0lep4rPrkqGYtVeOTlZsnzy3ZKZ+sP6K3698uXF7/Y3up4c3cmwBgxu9vgg5BBxUwd91v8tySHZKT55BWYf4y84Eu0rBuDc4lAJjs+5umK6AC7uveSD55pIcE1/KR3fFpMuDdVfLzvpOcSwAwGYIOUEFdGgfJ14/3lg6RtSXlXI48+MF6mfnTAfrtAICJmC7oTJs2TRo3biy+vr7SvXt3Wb9+/WW3f/vtt+Wqq64SPz8/iYyMlDFjxkhmJhO6oXqoTssLHu0hf+zcQNQFzl9etluemL9FzmXn8RYAgAmYKugsWLBAxo4dK5MnT5ZNmzZJhw4dpF+/fpKYWPJQ3nnz5smECRP09rt27ZL3339f/4ynnnqq2ssO9+Xr5SGvDWkvzw9sq+fZWbL1hNw5fY0cO5NhdNEAwO2ZqjOyqsHp2rWrvPvuu3o5Pz9f19I8/vjjOtBcbPTo0TrgrFixonDd3/72N1m3bp2sWrWqxH1kZWXpW9HOTGofjLqCM6w9eEpGzd0kp9KzJaimt7x779XSs1kwJxcA3L0zcnZ2tmzcuFH69OlTuM5ut+vl2NjYEl/Ts2dP/ZqC5q2DBw/KsmXL5Lbbbit1P1OnTtUnpuCmQg7gLD2a1pUlj/eWdhGBcjo9Wx54f718sOoQ/XYAwCCmCTpJSUmSl5cnoaGhxdar5fj4+BJfc++998rzzz8vvXv3Fi8vL2nWrJlcf/31l226mjhxok5/BbejR486/Vjg3iJq+8lnj0XL4Ksj9Hw8z3+9U/722VbJzKHfDgC4bdCpiJUrV8rLL78s7733nu7T88UXX8jSpUvlhRdeKPU1Pj4+uoqr6A2oin47b93VQSbd3kZPNPjFpuNy179j5UTyOU42ALhj0AkODhYPDw9JSEgotl4th4WFlfiaSZMmyQMPPCAjRoyQdu3ayeDBg3XwUc1Tqn8PYCSbzSYP924iH/2pm9Sp4SW/HkuRP7y7StYdPMUbAwDuFnS8vb2lc+fOxToWq7CilqOjo0t8TUZGhu7HU5QKS4qJ+ljDzfVqHixLRveW1uEBknQ2W+77zzr5b+xh/o8CgDsFHUUNLZ81a5bMmTNHj6YaOXKkpKeny/Dhw/Xzw4YN031sCgwYMECmT58u8+fPl0OHDsny5ct1LY9aXxB4ADOIDKohX4zsKQM61NdXQJ/05Q6Z8Pk2fe0sRfXliT1wSr7cclzfq2UAQOWZ6kqEQ4cOlZMnT8qzzz6rOyB37NhRYmJiCjsoHzlypFgNzjPPPKObB9T98ePHpV69ejrkvPTSSwYeBVAyP28P+dfdHaVdRIC88r/dsmDDUdmTkCZDu0bKv1bsk7iUzGJXS588oI2+mjoAwCLz6BiBi3rCCD/tPSmPf7JZXzqiJLYL99Pv70TYAQArzKMDuJNrW9aTRX/pqWdSLknBXx9TvtpJMxYAVAJBBzBIQmqW7q9TGvWMas5af+h0tZYLAKyEoAMYJDEt06nbAQAuRdABDBLi7+vU7QAAlyLoAAbp1iRIj64quZfO+Q7J6nm1HQCgYgg6gEHUpSHUEHKltLCjnlfbAQAqhqADGEjNk6OGkIcFFm+equXjydByALDahIGAu4advm3C9Oiqr389IXPXHZH6gb7MnwMATkCNDmACqnkqulld+Xu/q/TcOnsTz8rhpHSjiwUALo+gA5hI7Rre0qNpXf34mx3xRhcHAFweQQcwmX5tz1/bjaADAJVH0AFMRvXXUTYdSZbEVCYLBIDKIOgAJqNGYHWMrK0ff7szwejiAIBLI+gAJtSv7flaHZqvAKByCDqAifvpxB44JSnncowuDgC4LIIOYEJN69WSFiG19NXNf9idaHRxAMBlEXQAk6L5CgAqj6ADmDzorNxzUjJz8owuDgC4JIIOYFJREQESUdtPzuXkyc/7kowuDgC4JIIOYFI2m036tmHyQACoDIIO4ALNVyt2JUhuXr7RxQEAl0PQAUysa+M6UqeGl5zJyJH1h08bXRwAcDkEHcDEPD3s0qf1+earb3cwSzIAlBdBB3CR5qtvd8SLw+EwujgA4FIIOoDJ9W4RLDW8PeRESqZsO55idHEAwKUQdACT8/XykOuvqqcfc+0rACgfgg7gUrMk008HAMqDoAO4gBtahYiXh032J56VAyfPGl0cAHAZBB3ABQT4ekl0s2D9mOYrACg7gg7gIvq1LZglmeYrACgrgg7gItTlIGw2ka1HkyU+JdPo4gCASyDoAC4ixN9XOjWsox9/uzPe6OIAgEsg6AAu2XxF0AGAsiDoAC44zHztwdOSnJFtdHEAwPQIOoALaVS3prQK85e8fIes2JVodHEAwPQIOoCLublw8kCarwDgSgg6gIv20/lp30k5l51ndHEAwNQIOoCLaRMeIA3q+ElmTr78uPek0cUBAFMj6AAuxmazFXZK/pbmKwC4LIIO4IIKgs53uxIkJy/f6OIAgGkRdAAX1LlRHalb01tSM3Nl3cHTRhcHAEyLoAO4IA+7TV8SQmH0FQCUjqADuKjCfjo74yU/32F0cQDAlAg6gIvq2byu1PLxlITULNl6LNno4gCAKRF0ABfl4+kh119VTz/+ZkeC0cUBAFMi6AAurOgwc4eD5isAuBhBB3BhqkbH28MuB5PSZX/iWaOLAwCmQ9ABXJi/r5f0al5XP2b0FQBciqADWKT5in46AHApgg7g4vq0CRW7TWTb8RQ5dibD6OIAgKkQdAAXF1zLR7o0CtKPv2X0FQAUQ9ABLODmtsySDAAlIegAFuqn88vh03LqbJbRxQEA0yDoABYQGVRD2oQHiLoSxIpdiUYXBwBMg6ADWG70VbzRRQEA0yDoABbRL+p8P52f9yfJ2axco4sDAKZguqAzbdo0ady4sfj6+kr37t1l/fr1l90+OTlZRo0aJeHh4eLj4yMtW7aUZcuWVVt5AbO4KtRfGtWtIdm5+fLjnpNGFwcATMFUQWfBggUyduxYmTx5smzatEk6dOgg/fr1k8TEkvscZGdnS9++feXw4cOycOFC2bNnj8yaNUsiIiKqveyA0Ww2G81XAHARm8NEVwJUNThdu3aVd999Vy/n5+dLZGSkPP744zJhwoRLtp8xY4a8/vrrsnv3bvHy8irTPrKysvStQGpqqt5HSkqKBAQEOPFogOq38bfTcuf0WPH38ZSNk/qKt6ep/pYBAKdR39+BgYFX/P42zaegqp3ZuHGj9OnTp3Cd3W7Xy7GxsSW+ZsmSJRIdHa2brkJDQyUqKkpefvllycvLK3U/U6dO1Sem4KZCDmAVV0fWkXr+PpKWlStrDiQZXRwAMJxpgk5SUpIOKCqwFKWW4+NLHkVy8OBB3WSlXqf65UyaNEnefPNNefHFF0vdz8SJE3X6K7gdPXrU6ccCGMVut0nfNgWTBybwRgBwe6YJOhWhmrZCQkJk5syZ0rlzZxk6dKg8/fTTukmrNKrDsqriKnoDrDjMfPnOBMlTE+sAgBszTdAJDg4WDw8PSUgo/leoWg4LO//BfTE10kqNslKvK9C6dWtdA6SawgB3FN20rvj7ekrS2SzZfOSM0cUBAEOZJuh4e3vrWpkVK1YUq7FRy6ofTkl69eol+/fv19sV2Lt3rw5A6ucB7kh1QL6xVYh+zOSBANydaYKOooaWq+Hhc+bMkV27dsnIkSMlPT1dhg8frp8fNmyY7mNTQD1/+vRpeeKJJ3TAWbp0qe6MrDonA+7s91mSE8REAysBoNp5iomoPjYnT56UZ599Vjc/dezYUWJiYgo7KB85ckSPxCqgRkx98803MmbMGGnfvr2eP0eFnvHjxxt4FIDxrmtZT9fsHDmdIbvj06R1OH3RALgnU82jY+Zx+ICrGTHnF/luV6I82aeFPNmnpdHFAQD3nkcHgHPdXKT5CgDcFUEHsKg+rUPFbhPZFZcqR09nGF0cADAEQQewqKCa3tKtSZB+zOgrAO6KoAO4xeirkmcXBwCrI+gAbtBPZ8NvZ+Rk2u8XswUAd0HQASwsoraftIsIFDW28rtddEoG4H4IOoDF9WtbcJFPmq8AuB+CDuAm/XTW7D8laZk5RhcHAKoVQQewuOYhtaRpcE3JzsuXH/acNLo4AFCtCDqAxdlstiKTB9J8BcC9EHQAN+qns3J3omTm5BldHACoNgQdwA10aFBbQgN8JD07T9YcSDK6OADgOkHn9OnTkp+f75zSAKgSdrtNbm5zoflqO8PMAbiPCgWdnTt3yiuvvCI9e/aUevXqSUhIiAwbNkw+//xzSU9Pd34pATht9JWaTycv38EZBeAWyhx09uzZI3/729+kRYsW0qNHD/nll1/ksccek4SEBFm2bJk0atRInn/+eQkODpZbb71Vpk+fXrUlB1Au3ZsGSaCfl5xKz5YNh09z9gC4hTIHnTVr1ujamn/961+SlJSka29ULY4KNt26dZMXXnhBtm7dKrt27ZJbbrlFPw/APLw87HJTqxD9+JsdNF8BcA82h0NNDu++UlNTJTAwUFJSUiQgIMDo4gBVKmZ7vDz28UZ9aYhV42/QQ88BwMrf35XujHzu3Dk5fvz4Jet37NhR2R8NwMmua1lPfL3scjz5nOw4kcr5BWB5lQo6Cxcu1H12+vfvL+3bt5d169YVPvfAAw84o3wAnMjP20OubVFPP/6WyQMBuIFKBZ0XX3xRNm7cKFu2bJEPP/xQHn74YZk3b55+zs1bxADTj76inw4Ad+BZmRfn5ORIaOj5GVc7d+4sP/30kwwePFj2799P2z9gUje1DhG7TWRPQpq8//NBaVM/ULo1CRIPtRIALKZSNTpq/pxff/21cDkoKEiWL1+uR14VXQ/APNYePCWe9vO/+i8s3SX3zForvV/9XmK2xxldNAAw16irY8eOiaenp4SFna8KL2r16tXSq1cvMTtGXcGdqDAz8uNNcvEvfUFdzvT7O8ktUeEGlAwATDjqqkGDBiWGHMUVQg7gTtRsyFO+2nlJyFEK1qnnmTUZgJXYKzu0vOh1rlasWCFvvvmmnikZgLmsP3Ra4lIyS31ehR31vNoOAKyiUkFHXQpCVR0pU6dOlZdeekmPtpoxY4aMGzfOWWUE4ASJaZlO3Q4ALD/qKjc3V2rXrq0fq0s+qMtEeHt7y9ixY6Vjx47OKiMAJwjx93XqdgBg+RodNbRchRulfv36+hpYSlpamnNKB8Bp1BDy8EDfwo7HF1Pr1fNqOwCwikrV6HzwwQfy0EMP6VqcmjVrSocOHaRLly6SmJgor7/+uvNKCaDS1Dw5kwe00aOuVKgpqVOyep75dAC47fDye+65R55++mmJiooqtn7nzp2yb98+3ZQVEREhXbt2FQ8PD3EFDC+HOw4xV6OrLu6Y/NebmsvYvlcZVi4AqIrv73IFHbvdLvXq1dOjqy4OO4oagZWRkSG1atUSV0HQgTtSQ8jV6CrV8firrSfku12JMqhjfXn77quNLhoAGDuPjupkfOONN8r27dsvee7kyZOFnZMBmJdqnopuVlcGdoyQUTc01+tidsTL2axco4sGAE5VrqBjs9lk9uzZOuio27Zt2y7Zpui8OgDMr2NkbWkaXFMyc/Jl2TYuAwHAjYOOauVSfW/UFcpvuummEsOOCkMAXIf6nb2zcwP9+ItNx4wuDgAYP7xc9dWZO3eu9OnTR4cdLuAJuLZBV0eI+htl7cHTcvR0htHFAQDjmq4uDjt9+/bVtTuEHcB1RdT2k+imdfXjxZuPG10cADCu6arYi+12+fjjjwvDzpYtW5xXMgDV6s5OF5qvNh+/5HcdANwi6CxdulQP5Sop7Nx8881y5513Ort8AKrJLVFhUsPbQw4lpcumI8mcdwDuF3RuvfVW8fHxufSH2O3y3//+VwYOHMhfgoCLqunjqcOO8jmdkgG4W9CJj4+XrKys0n/QhZqdtWvX6uWDBw86p4QAqr356uutJyQzJ48zD8B9gs7ChQslKChIBg8eLB9++KGeHPBi69evl8WLF0vbtm31da8AuBbVIbl+oK+kZubKil2JRhcHAKov6IwePVq2bt0q11xzjZ40sEGDBtK7d295+eWX5ZFHHpHw8HAZNGiQvqDnK6+8UmIQAmBudrtNDzVXaL4CYAXlutZVUadOnZKvv/5ali1bJo0bN9b9c6Kjo11uwkCudQUUtz/xrPR560d9mYi1E2+Sev6X9ssDAEte1NOKCDrApQZOWy1bjybLM/1by4hrmnKKALjPRT0BWN+QTuebr77YxOSBAFwbQQfAJW5vX1+8PGyyMy5VdsWlcoYAuCyCDoBL1KnpLTe1CtWPudAnAFdG0AFQooIrmi/afEJy8/I5SwBcEkEHQImua1lPgmp6S9LZLPl5fxJnCYBLIugAKJG3p13+0KG+fvz5xmOcJQAuiaAD4IqXhPh2Z4KknMvhTAFwOQQdAKWKigiQlqG1JDs3X5Zti+NMAXA5BB0ApVIznd9xoVaH5isAroigA+CyBl8dIXabyIbfzshvp9I5WwBcCkEHwGWFBvhKr+bB+vHnzJQMwMUQdABc0ZALc+qoyQPz89368ngAXIzpgs60adP01dB9fX2le/fusn79+jK9bv78+bo/waBBg6q8jIC7ublNmNTy8ZRjZ87JL4dPG10cAHDNoLNgwQIZO3asTJ48WTZt2iQdOnSQfv36SWJi4mVfd/jwYRk3bpxcc8011VZWwJ34eXvIbe3C9GMu9AnAlZgq6Lz11lvyyCOPyPDhw6VNmzYyY8YMqVGjhnzwwQelviYvL0/uu+8+mTJlijRt2rRaywu445w6S7fFybnsPKOLAwCuFXSys7Nl48aN0qdPn8J1drtdL8fGxpb6uueff15CQkLk4YcfLtN+srKyJDU1tdgNwJV1bRwkDer4ydmsXPl2ZzynDIBLME3QSUpK0rUzoaHnr5hcQC3Hx5f8obpq1Sp5//33ZdasWWXez9SpUyUwMLDwFhkZWemyA+7Abi8ypw6jrwC4CNMEnfJKS0uTBx54QIec4ODzQ1/LYuLEiZKSklJ4O3r0aJWWE7CSO66O0Per9p2UhNRMo4sDAFfkKSahwoqHh4ckJCQUW6+Ww8LOd4Is6sCBA7oT8oABAwrX5efn63tPT0/Zs2ePNGvW7JLX+fj46BuA8mscXFO6NKqjJw9cvPm4/Pm6S3/HAMBMTFOj4+3tLZ07d5YVK1YUCy5qOTo6+pLtW7VqJdu2bZMtW7YU3v7whz/IDTfcoB/TJAVUjd+br46Jw8GcOgDMzTQ1OooaWv7ggw9Kly5dpFu3bvL2229Lenq6HoWlDBs2TCIiInQ/GzXPTlRUVLHX165dW99fvB6A8/RvHy7PfbVD9iaclR0nUiUqIpDTC8C0TBV0hg4dKidPnpRnn31Wd0Du2LGjxMTEFHZQPnLkiB6JBcA4gX5e0rdNqCz9NU7X6hB0AJiZzeHmdc9qeLkafaU6JgcEBBhdHMAl/LA7UYbP/kXq1vSWtU/dJF4e/AECwJzf33w6ASi3a1oES3AtHzmVni0r95zkDAIwLYIOgHLz9LDLoI71Cy/0CQBmRdABUCF3Xrii+YpdiZKckc1ZBGBKBB0AFdI6PEDfsvPy5atf4ziLAEyJoAOgwu7sdH6m5M830nwFwJwIOgAqbGDHCPGw22TL0WQ5cPIsZxKA6RB0AFRYPX8fua5lPf2YTskAzIigA6BS7rjQfLVo03HJz3frabkAmBBBB0Cl9GkdKv6+nnIiJVPWHjzF2QRgKgQdAJXi6+Uht7c/P6fO55uOczYBmApBB0ClDel8vvnqf9vjJD0rlzMKwDQIOgAqrVPDOtK4bg3JyM6Tb3bEc0YBmAZBB0Cl2Ww2uaPT+ZmS1RXNAcAsCDoAnGLw1eebr9YcOCUnks9xVgGYAkEHgFNEBtWQ7k2CxOEQWbSZTskAzIGgA8DpF/pUzVcOlXgAwGAEHQBOc2tUmPh62eXgyXTZeiyFMwvAcAQdAE7j7+slt7QN04+50CcAMyDoAHCqgtFXX/16QrJy8zi7AAxF0AHgVL2aB0togI8kZ+TID7sTObsADEXQAeBUHnabDLow1JxLQgAwGkEHgNPdeaH5StXonDqbxRkGYBiCDgCnaxnqL+0iAiU33yFfbT3BGQZgGIIOgCpxZyearwAYj6ADoEoM6FBfPO022XY8RfYmpHGWARiCoAOgStSt5SM3tArRj7nQJwCjEHQAVHmn5MWbj0tePpeEAFD9CDoAqswNrepJ7RpekpCaJav3J3GmAVQ7gg6AKuPj6SF/6FBfP6b5CoARCDoAquWSEMu2xcmCX45I7IFTNGMBqDae1bcrAO4oLvmcni05J88h4z/fpteFB/rK5AFt5JaocKOLB8DiqNEBUGVitsfJX+ZuuqQGJz4lU0Z+vEk/DwBViaADoEqocDPlq51S0lirgnXqeUZjAahKBB0AVWL9odMSl5JZ6vMq7Kjn1XYAUFUIOgCqRGJaplO3A4CKIOgAqBIh/r5O3Q4AKoKgA6BKdGsSpEdX2S6zjbeHXdrWD+AdAFBlCDoAqoQaUq6GkCulhZ3svHwZPvsXSc3M4V0AUCUIOgCqjJonZ/r9nSQssHjzlKrpGX9LKwnw9ZSNv52R+2atkzPp2bwTAJzO5nA43PpKe6mpqRIYGCgpKSkSEEAVOlAV1BByNbpKdTxWfXJUs5aq8dlxIkUeeH+9nE7PllZh/vLxiO4SXMuHNwGA076/CToEHcBQ+xLS5N7/rJOTaVnSrF5NmTuixyU1QABQ0e9vmq4AGKpFqL98+udoqR/oKwdOpstd/46VY2cyeFcAOAVBB4DhmgTXlAV/jpaGQTXkyOkMuWtGrBxOSje6WAAsgKADwBQig2romp2m9WrKiZRMXbOjmrUAoDIIOgBMQ/XNWfBotO6YnJiWJUNnrtUdlgGgogg6AEylnr+PfPJID2kXEahHY90zc61sOZpsdLEAuCiCDgDTqVPTW+Y+0l06N6ojqZm5cv9/1skvh7n4J4DyI+gAMKUAXy/56E/dpEfTIDmblSvD3l8vq/cnGV0sAC6GoAPAtGr6eMrs4d3k2pb15FxOnr5cxA+7E40uFgAXQtABYGq+Xh4ya1hn6dsmVLJz8+XR/26QmO1xRhcLgIsg6AAwPR9PD3nvvk7Sv3245OQ5ZNS8zfLlluNGFwuACyDoAHAJXh52+dfdV8udnRroa2c9uWCLfPrLUaOLBcDkCDoAXIa6EOjrQ9rLfd0biroc8T8+/1U+ij1sdLEAmBhBB4BLsdtt8uKgKPlTryZ6+dkvd8jMnw4YXSwAJkXQAeBybDabTLq9tYy6oZlefnnZbvm/7/aJQ1XzAEARBB0ALht2/t6vlYy7uaVe/ud3e+W1b/YQdgCYO+hMmzZNGjduLL6+vtK9e3dZv359qdvOmjVLrrnmGqlTp46+9enT57LbA7Ce0Te2kGf6t9aPp688IFO+2qnDjuqwHHvglB6dpe7VMgD34ykmsmDBAhk7dqzMmDFDh5y3335b+vXrJ3v27JGQkJBLtl+5cqXcc8890rNnTx2MXn31Vbn55ptlx44dEhERYcgxAKh+I65pKj5eHjJp8XaZveaw7E9Mk/2J6RKfmlm4TXigr0we0EZuiQrnLQLciM1hokZtFW66du0q7777rl7Oz8+XyMhIefzxx2XChAlXfH1eXp6u2VGvHzZsWJn2mZqaKoGBgZKSkiIBAQGVPgYAxvlsw1H5x8JfpaQPNduF++n3dyLsABZQ1u9v0zRdZWdny8aNG3XzUwG73a6XY2Njy/QzMjIyJCcnR4KCgkrdJisrS5+cojcA1nBHpwYS6OdV4nMF4Uc1bdGMBbgP0wSdpKQkXSMTGhpabL1ajo+PL9PPGD9+vNSvX79YWLrY1KlTdQIsuKkaIwDWsP7QaUk+l1Pq8yrsxKVk6u0AuAfTBJ3KeuWVV2T+/PmyaNEi3V+nNBMnTtTVXAW3o0eZWRWwisS0TKduB8D1maYzcnBwsHh4eEhCQkKx9Wo5LCzssq994403dND57rvvpH379pfd1sfHR98AWE+Iv69TtwPg+kxTo+Pt7S2dO3eWFStWFK5TnZHVcnR0dKmve+211+SFF16QmJgY6dKlSzWVFoAZdWsSpEdXFXQ8vphar55X2wFwD6YJOooaWq7mxpkzZ47s2rVLRo4cKenp6TJ8+HD9vBpJpZqeCqjh5JMmTZIPPvhAz72j+vKo29mzZw08CgBGXgtLDSFXbKX00Xn29jZ6OwDuwVRBZ+jQoboZ6tlnn5WOHTvKli1bdE1NQQflI0eOSFxcXOH206dP16O1hgwZIuHh4YU39TMAuCc1T44aQh4WWHLz1MmzWdVeJgDGMdU8OkZgHh3AmtQQcjW6SnU8Vn1yth1P1tfE8vawy+cje0q7BoFGFxFANXx/m6YzMgA4k2qeim5Wt3C5R9Mg+eXwGVm+M0FGzdskX/+1twT4ljznDgDrMFXTFQBU5UVA3xjSQRrU8ZMjpzPkH5/9ygVAATdA0AHgNgJreMm793YSLw+bxOyI19fFAmBtBB0AbqVjZG2ZeOv5q52/vGyXbD2abHSRAFQhgg4AtzO8V2O5pW2Y5OQ5dH+dlIzSLxsBwLURdAC4ZX+dV4e0l8ggPzl25pyMW7iV/jqARRF0ALgldZXz9+7trIebq5FYH6ymvw5gRQQdAG5LzaXzdP/z/XWmLtslm4+cMbpIAJyMoAPArQ2LbiT924VLbr5DRs/bLMkZ2UYXCYATEXQAiLv315l6ZztpVLeGHE8+J+M+o78OYCUEHQBuT82QPO3eTrq/zne7EuU/Px9y+3MCWAVBBwBEJCoiUCZduPL5qzG7ZeNv9NcBrICgAwAX3N+9odze/nx/ncfnbZIz6fTXAVwdQQcAivbXuaOdNAmuKSdSMuVvn22V/HwH5wdwYQQdACjCv6C/jqddvt+dKDN/Psj5AVwYQQcALtKmfoA8N6Ctfvz6N3vkl8OnOUeAiyLoAEAJ7ukWKQM71pc83V9ns5w6m8V5AlwQQQcASumv8/LgdtK0Xk2JT82UsZ/SXwdwRQQdAChFTR9Pee++TuLjaZcf956U6T8e4FwBLoagAwCX0SosQJ4feL6/zpvf7pF1B09xvgAXQtABgCu4q0uk3HF1hKiR5o9/slmS6K8DuAyCDgCUob/Oi4OjpHlILUlMy5IxC7Ywvw7gIgg6AFAGNbzP99fx9bLLz/uSZNoP+zlvgAsg6ABAGbUM9ZcXBkbpx//8bq+sOZDEuQNMjqADAOXwxy6RMqRzA91f54n5W+RkGvPrAGZG0AGAclK1Oi1Da+mQ8+SCzXpSQQDmRNABgHLy8/bQ/XX8vDxk9f5T8s73+3TYiT1wSr7cclzfE34Ac7A5HA63/lMkNTVVAgMDJSUlRQICAowuDgAX8sWmY3rGZCWoprecTs8ufC480FcmD2gjt0SFG1hCwLrK+v1NjQ4AVNAdnRpIr2Z19eOiIUeJT8mUkR9vkpjtcZxfwEAEHQCoINU8deDk2RKfK6gqn/LVTpqxAAMRdACggtYfOi3xqaWPulJhJy4lU28HwBgEHQCooMS0TKduB8D5CDoAUEEh/r5O3Q6A8xF0AKCCujUJ0qOrbJfZpm5Nb70dAGMQdACggjzsNj2EXCkt7CSfy5ZPNxzlHAMGIegAQCWoeXKm399JwgKLN0+FBfhKp4a1JS9fZOIX2+SZxdskOzefcw1UMyYMZMJAAE4aaq5GV6mOx6pPjmqusttE3lt5QN74do+oqVnVOjWjcnAtH845UE3f3wQdgg6AKrZiV4K+AOjZrFypH+grM4d1kaiIQM47UAnMjAwAJnFT61BZPKqnNAmuKSdSMmXIjDWyZOsJo4sFuAX66ABANWge4i+LR/WS61rWk8ycfPnrJ5vl1ZjdzJoMVDGCDgBUk0A/L/ngoa7y2HXN9PL0lQfk4Tm/SMq5HN4DoIoQdACgmoekT7i1lfzf3R3Fx9MuK/eclMHTVsv+xJKvmQWgcgg6AGCAgR0j5PORPXXn5INJ6TrsqE7LAJyLoAMABlEjr5Y83lu6NQ6StKxcGfHRBpn2w35xqLHoAJyCoAMABlJz6nw8orvc36Ohnmvn9W/2yOhPNktGdi7vC+AEBB0AMJi3p11eHNROXh7cTrw8bLL01zi5c3qsHD2dYXTRAJdH0AEAk7i3e0OZ90gPCa7lLbviUmXgtNUSe+CU0cUCXBpBBwBMpGvjIFkyurdERQTI6fRsuf/9dfJR7GH67QAVRNABAJOpX9tPFj7WUwZ2rK8nFHz2yx0y4fNtkpWbZ3TRAJdD0AEAE/L18pC3h3aUp25rpS8OumDDUbln5lpJTM3Uz6sApJq1vtxyXN+rZQCX4qKeXNQTgMn9uPekPD5vk6Rm5kpYgK881LOxzIk9LHEp50OPEh7oK5MHtJFbosINLStQXbh6uZNPFAAY6VBSujzy0YZSZ1C2Xbiffn8nwg7cQmoZv79pugIAF6CufL7wsWh92YiSFDRcTflqJ81YQBEEHQBwEbvi0iQrN7/U51XYUc1Z6w+drtZyAWbmaXQBAABlk5j2e5+cyxn32VaJblZX2oQHSJv6AdI6PEBfOR1wRwQdAHARIf6+ZdruePI5WbjxWLF1EbX9dOgpCD/qvkEdP7HZCnr3XJ4a1aVqilTYUuXo1iRIX4m9qlT3/ozYJ/tz06Azbdo0ef311yU+Pl46dOgg77zzjnTr1q3U7T/77DOZNGmSHD58WFq0aCGvvvqq3HbbbdVaZgCoDuqLV42uik/JLOyTU5T6Sq7n7yNT/tBWdsenyc64VNl5IlUHn4Lb8p2/XyE9wNdT1/YUDUAtQvz1JSmKitkep/v+VNcor+renxH7ZH9uOrx8wYIFMmzYMJkxY4Z0795d3n77bR1k9uzZIyEhIZdsv2bNGrn22mtl6tSpcvvtt8u8efN00Nm0aZNERUWVaZ+MugLgStQX5MiPN+nHjjKOukrJyDkfei4EH3W/PzFNcvIu/fhX19pqHuJfGHzOZubK29/tvSRYVdUor4Ljq679GbFP9tfJKefTJYeXq3DTtWtXeffdd/Vyfn6+REZGyuOPPy4TJky4ZPuhQ4dKenq6fP3114XrevToIR07dtRhqSwIOgBcjTNqA7Jz8/VQ9d/DT4q+V3P1lIe/r6eMvL6ZeHvYxW6ziaeHTTf3eNjO36tlvd5uFw+7iIfdLp52m9jVc2q7IjeVNB797wZJOptd6v5C/H1kwaPRlzQplbEFrsTmoz/+O1ZOpmVddp+fPXbpPiu6vyEz3Hd/NhEJC/SVVeNvrPT+XC7oZGdnS40aNWThwoUyaNCgwvUPPvigJCcny5dffnnJaxo2bChjx46VJ598snDd5MmTZfHixbJ169YS95OVlaVvRU+UClPMowPAlVRF/w71daCatwpqfVbtS5INv51xWpmBAp880kN3mK+OoGOaPjpJSUmSl5cnoaGhxdar5d27d5f4GtWPp6Tt1frSqGauKVOmOKnUAGAMFWoq+0VxMdUxuUGdGvp2c9swPXdPWYJOt8ZBElHHT3LzHZKXn69DmLqdXy5+U+vyHQ7JzbtwX+S5tMwcOZORc8X9+XiomiJ7if2UClzuT3hHkVeq/ZbUhHexgtqnyio4B+6+v8QyjiB0BtMEneoyceJEXQt0cY0OAKBio7zG9G3plNClrtl1z6y1V9xu9p+6Oy3klXWf/33YOftkf+X7v2WpCQODg4PFw8NDEhJ+HxGgqOWwsLASX6PWl2d7xcfHR1dxFb0BAEof5VXa3/lqvXpebeeK+zNin+xPnP4eukzQ8fb2ls6dO8uKFSsK16nOyGo5Ojq6xNeo9UW3V5YvX17q9gCAslNNGaqDs3JxEChYVs87a66Z6t6fEftkf+L099Blgo6impRmzZolc+bMkV27dsnIkSP1qKrhw4fr59XQc9X0VOCJJ56QmJgYefPNN3U/nueee042bNggo0ePNvAoAMA61CguNbxajZQpSi1XxVDv6t6fEftkf+FSnUwz6qqAGlpeMGGgGib+r3/9Sw87V66//npp3LixzJ49u3B7Nc/OM888Uzhh4GuvvVauCQMZXg4AV2b1WXyN2Cf7E/caXm4Ugg4AANb9/jZV0xUAAIAzEXQAAIBlEXQAAIBlEXQAAIBlEXQAAIBlEXQAAIBlEXQAAIBlEXQAAIBlEXQAAIBleYqbK5gYWs2wCAAAXEPB9/aVLvDg9kEnLS1Nn4jIyMhqeWMAAIBzv8fVpSBK4/bXusrPz5cTJ06Iv7+/2Gw2pyZNFZ6OHj162WtwuDp3OE6O0Tp4L62D99I6Uiv4PaJqclTIqV+/vtjtpffEcfsaHXVyGjRoIFVFvWlWDQDudpwco3XwXloH76V7v5eBl6nJKUBnZAAAYFkEHQAAYFkEnSri4+MjkydP1vdW5g7HyTFaB++ldfBeWodPFX+PuH1nZAAAYF3U6AAAAMsi6AAAAMsi6AAAAMsi6AAAAMsi6FTCtGnTpHHjxuLr6yvdu3eX9evXX3b7zz77TFq1aqW3b9eunSxbtkysdpyzZs2Sa665RurUqaNvffr0ueJ5ccX3ssD8+fP1jNqDBg0Sqx1jcnKyjBo1SsLDw/VoiJYtW7rE/9nyHufbb78tV111lfj5+enZWceMGSOZmZliVj/99JMMGDBAzwar/u8tXrz4iq9ZuXKldOrUSb+PzZs3l9mzZ4uZlfcYv/jiC+nbt6/Uq1dPTzgXHR0t33zzjZhdRd7LAqtXrxZPT0/p2LGjWO0Ys7Ky5Omnn5ZGjRrp/7Pq9/mDDz6ocBkIOhW0YMECGTt2rB4St2nTJunQoYP069dPEhMTS9x+zZo1cs8998jDDz8smzdv1l+M6rZ9+3ax0nGqD1R1nD/88IPExsbqL46bb75Zjh8/LlY5xgKHDx+WcePG6WBnduU9xuzsbP3FoY5x4cKFsmfPHh1iIyIixErHOW/ePJkwYYLefteuXfL+++/rn/HUU0+JWaWnp+vjUoGuLA4dOiT9+/eXG264QbZs2SJPPvmkjBgxwtRBoLzHqL5M1f9XFcQ3btyoj1V9uarPWjMr73EW/SNk2LBhctNNN4nZpVfgGO+66y5ZsWKF/n1Unz2ffPKJ/mOkwhyokG7dujlGjRpVuJyXl+eoX7++Y+rUqSVuf9dddzn69+9fbF337t0df/7zny11nBfLzc11+Pv7O+bMmeOw0jGq4+rZs6fjP//5j+PBBx90DBw40GFm5T3G6dOnO5o2berIzs52uJLyHqfa9sYbbyy2buzYsY5evXo5XIH6CF+0aNFlt/nHP/7haNu2bbF1Q4cOdfTr189hlWMsSZs2bRxTpkxxuIryHKd6/5555hnH5MmTHR06dHBY6Rj/97//OQIDAx2nTp1y2n6p0akA9deu+qtBNcsUvWaWWla1GCVR64tur6i/NEvb3lWP82IZGRmSk5MjQUFBYqVjfP755yUkJETX0JldRY5xyZIluvpfNV2FhoZKVFSUvPzyy5KXlydWOs6ePXvq1xQ0bx08eFDXCtx2221iFa742eOMizWriz2a9XOnMj788EP9/1TVQlrRkiVLpEuXLvLaa6/pGmTVZK5qzs+dO1fhn+n2F/WsiKSkJP2Br74AilLLu3fvLvE18fHxJW6v1lvpOC82fvx43TZ78QetKx/jqlWrdJWqagZwBRU5RvVB+v3338t9992nv/j3798vf/nLX3RoNesHbEWO895779Wv6927t74Scm5urjz22GOmbroqr9I+e9QVo9WXh+qbZDVvvPGGnD17VjeBWMm+fft0U+vPP/+s++dY0cGDB/VnrOpjt2jRIv37qT57Tp06pUNeRVCjgyrzyiuv6M666j+r+k9rBeqvxAceeED3VwkODhYr/0WsaqxmzpwpnTt3lqFDh+rOgTNmzBArUX3KVE3Ve++9p/v0qE6tS5culRdeeMHooqGCVL+rKVOmyKeffqr/D1uFCvEqmKtjU7UcVv7ssdlsMnfuXOnWrZuuXX3rrbdkzpw5Fa7VsWYkrGLqC87Dw0MSEhKKrVfLYWFhJb5GrS/P9q56nEX/olJB57vvvpP27duLVY7xwIEDuoOu6uhY9BdTUX9hqY5zzZo1E1d/H9VIKy8vL/26Aq1bt9a1A6qJyNvbW8ymIsc5adIkHVxV51xFjYZUnScfffRRHexU05erK+2zR41OslptjvrDSr2XaoSrWWuRK/NH1oYNG3QH69GjRxd+9qiaSPXZ8+2338qNN94ori48PFw3WQUGBhb77FHHeezYMWnRokW5f6br/xYbQH3Iq79yVa/wAuo/nFpW/RpKotYX3V5Zvnx5qdu76nEqqm1V/UUcExOj21rNrLzHqKYH2LZtm262Krj94Q9/KBzRokaZWeF97NWrl26uKghxyt69e/WHkBlDTkWPU/UhuzjMFIS7830nXZ8rfvZUhBqZM3z4cH2vRplZjQqmF3/2qGZWNRpJPVZTKVhBr1695MSJE7rpsehnj/o9bdCgQcV+qNO6NbuZ+fPnO3x8fByzZ8927Ny50/Hoo486ateu7YiPj9fPP/DAA44JEyYUbr969WqHp6en44033nDs2rVL95b38vJybNu2zWGl43zllVcc3t7ejoULFzri4uIKb2lpaQ6rHOPFXGHUVXmP8ciRI3q03OjRox179uxxfP31146QkBDHiy++6LDScarfQ3Wcn3zyiePgwYOOb7/91tGsWTM9StKs1O/S5s2b9U19hL/11lv68W+//aafV8enjrOAOq4aNWo4/v73v+vPnmnTpjk8PDwcMTExDqsc49y5c/Xnqzq2op87ycnJDjMr73FezBVGXaWV8xjV9g0aNHAMGTLEsWPHDsePP/7oaNGihWPEiBEVLgNBpxLeeecdR8OGDfUXuxrWunbt2sLnrrvuOv0FWNSnn37qaNmypd5eDfdcunSpw2rH2ahRI/2f+eKb+oW00nvpakGnIse4Zs0aPQWCCg5qqPlLL72kh9Vb6ThzcnIczz33nA43vr6+jsjISMdf/vIXx5kzZxxm9cMPP5T4O1ZwXOpeHefFr+nYsaM+J+q9/PDDDx1mVt5jVI8vt72V3ktXCzo/VOAYVSDv06ePw8/PT4ceNeVDRkZGhctgU/84r9IJAADAPOijAwAALIugAwAALIugAwAALIugAwAALIugAwAALIugAwAALIugAwAALIugAwAALIugAwAALIugAwAALIugA8CtXX/99fLkk09W+ueoq+m89dZb0qRJE6lRo4YMGjRIUlJSnFJGABVH0AHgNNddd53YbDZ98/b2ltatW8u8efPc4gz//e9/l+nTp8ucOXPk559/lo0bN8pzzz1ndLEAt0fQAeAUqkZj8+bN8sYbb0hcXJzs2bNHbrnlFhk2bJgcOnTI0md53bp1ujZnwYIFcu2110rnzp3lkUcekWXLlhldNMDtEXQAOMW+ffskLS1Nh5uwsDDdhPPwww9LXl6eDj3lFRMTI71795batWtL3bp15fbbb5cDBw4Ua3L661//Kv/4xz8kKChI7/PiGhRVnvvuu09q1qwp4eHh8s9//vOKTVX5+fkydepUXX4/Pz/p0KGDLFy48LJlVeHupptukk6dOhWuCw0NlaSkpHIfNwDnIugAcArVVFOnTh1p06aNXj527Jg8/fTT4uPjI+3bty/3z0tPT5exY8fKhg0bZMWKFWK322Xw4ME6iBRQzUQqxKgalddee02ef/55Wb58eeHz6vWrV6+WJUuW6PWqSWnTpk2X3a8KOR999JHMmDFDduzYIWPGjJH7779ffvzxxxK3z8rKkqVLl+qyFZWZmSmBgYHlPm4ATuYAACcYN26cw263O2rWrOnw9fV1qI8XPz8/x4cffuiU83vy5En9M7dt26aXr7vuOkfv3r2LbdO1a1fH+PHj9ePU1FSHl5eX47PPPit8Pjk52VGjRg3HE088UbhO/ZyC5czMTP38mjVriv3chx9+2HHPPfeUWC61rSqXOmZ17AU3b29vR79+/fQ2ixcvLrZPANXH09nBCYB7UjUlo0aN0s1JycnJMm7cOOnVq5c89NBDFW4Ke/bZZ3VtjWoCKqjJOXLkiERFRenHF9cUqeapxMRE/fjgwYOSk5Mj3bp1K3xe1bBcddVVpe5z//79kpGRIX379i22Pjs7W66++uoSX7N3715dq7Rly5Zi6/v376+PX/n111+lY8eO5TwDAJyBoAPAaUFHdcBt3ry5Xn7vvfd0EFHrGjdurNepL37Vabd79+66/44KLKppqCQDBgyQRo0ayaxZs6R+/fo66KjtVego4OXlVew1arRX0aat8jp79qy+V01RERERxZ5TTXAlSU1NleDg4MLjVn777Tcd1O68887CoHPu3Dl9/KqjtmpKKwhrAKoWfXQAVJqqPVG1OEW/vFVfnWbNmhUbXj5p0iR55ZVXdNhRfW5KCzmnTp3SHZifeeYZ3clXDVM/c+ZMucrUtGlTHYR++eWXwnVqXhtVA1MaVWYVaFStkQouRW+RkZElvkaFHPVz1aizAi+99JLcdttthf2VVNBRYU/1F1I1Xl9++WW5jgVAxVGjA8ApHZFVqGjZsmWx9SqkLFq0SJ566im9rEZkqQ7KqsZEjaoqjerUrEZazZw5UzdHqeAxYcKEcpXJ399fHnzwQT2/jRqVFRISIpMnT9YBS9X8lPYa1eSmApiqGVKjvlSIUQElICBA/7yL3XjjjbrjsQpwd999t8ydO1e++uorWb9+vX5e1eSokWePPvqoXlbNaWokGYDqQY0OAKc0W7Vo0UJPElhUnz59dAhSI7AUVbty+vRp3Vfm4manYh9MdrvMnz9fv7ageev1118vd7lUzVF0dLQemq7KopqOVO2Qr69vqa954YUXdM2TGn2ltlXhTAUzNdy8JGoY+ezZs/VkgW3btpW1a9fKqlWrCmuAtm/fLl26dCncftu2bXo7ANXDpnokV9O+ALix48eP6+acxYsX674ragh3dfdTUUPWVd+bN998U/cRqg7vv/++xMfH65osRXVq/u6773SNFYCqR40OgCqnmm/++Mc/yjvvvKNrRiZOnKhrTqqamqn5k08+0RMNqlonNXmgMnDgQKkuqn9Ou3bt9OPc3Fzdl4mQA1QfanQAWJYKOiNGjNAdm1Wzmro0g2rOKggeAKyPoAMAACyLpisAAGBZBB0AAGBZBB0AAGBZBB0AAGBZBB0AAGBZBB0AAGBZBB0AAGBZBB0AAGBZBB0AAGBZBB0AACBW9f8zf/nV++sdzwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.xlabel(r\"$R_x$ angle $\\theta_h$\")\n",
    "plt.ylabel(r\"$\\langle Z_{62} \\rangle$\")\n",
    "plt.plot(hs, approx_evs, marker=\"o\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}