docs/tutorials/10_fermions_spin_charge_sep.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Spin-charge separation for fermions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook we take a look at the separation of spin and charge dynamics for fermionic atoms in a 1D system. It is based on the paper by \"Time-Resolved Observation of Spin-Charge Deconfinement in Fermionic Hubbard\n",
    "Chains\" the Munich group [[1]](#References).  In these experiments, fermions evolve in a one dimensional optical lattice. The experiment allows for site resolved observation of spin and charge, such that the famous spin charge separation is directly observable. It is therefore a very pedagogical example for microscopic observables in an extended system.\n",
    "\n",
    "Here, the physical dynamics of the 1-D tweezer array are governed by a Fermi-Hubbard Hamiltonian in second quantization: \n",
    "\n",
    "$$ \n",
    " H_{\\text{FH}}(\\boldsymbol{J},U,\\boldsymbol{\\mu}) =  -J \\underbrace{\\sum_{i=1,\\sigma}^{L-1} (f^\\dagger_{i,\\sigma} f_{i+1,\\sigma} + f^\\dagger_{i+1,\\sigma} f_{i,\\sigma} )}_{\\text{Tunneling/Hopping}} + \\underbrace{U \\sum_{i=1}^{L}  n_{i,\\uparrow}n_{i,\\downarrow}}_{\\text{interaction}} + \\underbrace{\\sum_{i=1,\\sigma}^{L} \\mu_i n_{i,\\sigma}}_{\\text{potential offset}} \n",
    "$$\n",
    "\n",
    "Here, $f_{i,\\sigma}, f^\\dagger_{i,\\sigma}$ are annihilation/creation operators for atoms in tweezers at site $i$ with spin $\\sigma$ and $n_{i,\\sigma} = f^\\dagger_{i,\\sigma} f_{i,\\sigma}$ is the number operator.  \n",
    "\n",
    "The dynamics depend on the parameters $\\{ \\boldsymbol{J}, U, \\boldsymbol{\\mu} \\}$ that determine the strength of the different contributions:\n",
    "\n",
    "- The first term describes hopping between the site $i$ and its neighboring site $i+1$ (for both spin species) through the tunnel effect. The parameter $J$ determines the strength of this hopping and can be tuned by adjusting the depth of the optical lattice.\n",
    "    \n",
    "- The second term describes an interaction between two atoms when they occupy the same site. It is controlled globally by the parameter $U$ set by an external magnetic field exploiting a Feshbach resonance.\n",
    "\n",
    "- The third term describes local offsets in the potential which can locally imprint a phase by tuning $\\mu_i$. Throughout this tutorial we will set $\\mu_i = 0$ for all sites $i$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initial state preparation\n",
    "\n",
    "As described in the paper, the initial state is a system of alternating spin-up and spin-down fermions in the lattice. To initialize the system only the site in the middle is left empty."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit_cold_atom.fermions import FermionSimulator, FfsimBackend\n",
    "\n",
    "from qiskit_cold_atom.fermions.fermion_gate_library import FermiHubbard\n",
    "from qiskit_nature.operators.second_quantization import FermionicOp\n",
    "from qiskit.circuit import Parameter\n",
    "\n",
    "backend = FermionSimulator()\n",
    "backend_ffsim = FfsimBackend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will attempt to obverve propagations of up to `d` sites. Let us initialize a system with `2d+1` sites. And in the middle of the lattice we will leave one site empty."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def prep_init_states(dmax: int, with_hole: bool) -> tuple[list, list]:\n",
    "    \"\"\"\n",
    "    Prepares a list of initial states for the Fermi-Hubbard model.\n",
    "\n",
    "    Args:\n",
    "        dmax: the maximum distance from the center\n",
    "        with_hole: whether to include a hole in the initial state\n",
    "\n",
    "    Returns:\n",
    "        two lists of initial states starting with spin-up and spin-down, respectively\n",
    "    \"\"\"\n",
    "    Nsites = 2 * dmax + 1\n",
    "\n",
    "    # the up state\n",
    "    state_up = [i % 2 for i in range(Nsites)]\n",
    "    state_down = [(i + 1) % 2 for i in range(Nsites)]\n",
    "\n",
    "    if with_hole:\n",
    "        state_up[dmax] = 0\n",
    "        state_down[dmax] = 0\n",
    "\n",
    "    init_state_1 = [state_up, state_down]\n",
    "    init_state_2 = [state_down, state_up]\n",
    "\n",
    "    return init_state_1, init_state_2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "dmax = 3\n",
    "Nsites = 2 * dmax + 1\n",
    "\n",
    "init_state_1, init_state_2 = prep_init_states(dmax, with_hole=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "hop_param = Parameter(\"hop_param\")\n",
    "int_param = Parameter(\"int_param\")\n",
    "# initialize one spin-up and one spin-down atom in the left tweezer\n",
    "qc_1 = backend.initialize_circuit(init_state_1)\n",
    "\n",
    "# apply a global Fermi-Hubbard gate with interaction\n",
    "qc_1.fhubbard(j=hop_param * np.ones(Nsites - 1), u=int_param, mu=np.zeros(Nsites), modes=range(2*Nsites))\n",
    "qc_1.measure_all()\n",
    "\n",
    "# qc_1.draw(output=\"mpl\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, this is not the only configuration that exists. Actually the total spin is on average zero. The only thing that is known is that we start out with an anti-ferromagnetic configuration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "hop_param_down = Parameter(\"hop_param\")\n",
    "int_param_down = Parameter(\"int_param\")\n",
    "\n",
    "# initialize one spin-up and one spin-down atom in the left tweezer\n",
    "qc_2 = backend.initialize_circuit(init_state_2)\n",
    "\n",
    "# apply a global Fermi-Hubbard gate with interaction\n",
    "qc_2.fhubbard(\n",
    "    j=hop_param_down * np.ones(Nsites - 1),\n",
    "    u=int_param_down,\n",
    "    mu=np.zeros(Nsites),\n",
    "    modes=range(2 * Nsites),\n",
    ")\n",
    "qc_2.measure_all()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "circuit_up_init = qc_1.bind_parameters({hop_param: 0, int_param: 0})\n",
    "job_up_init = backend.run(circuit_up_init, shots=5, num_species=2)\n",
    "\n",
    "circuit_down_init = qc_2.bind_parameters({hop_param_down: 0, int_param_down: 0})\n",
    "job_down_init = backend.run(circuit_down_init, shots=5, num_species=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let us plot up the density and charge density of the different states."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def nocc_from_counts(counts: dict[str, int], Nsites: int) -> np.ndarray:\n",
    "    \"\"\"\n",
    "    The calculates the occupation from all observed configurations.\n",
    "\n",
    "    Args:\n",
    "        counts: the counts of the observed configurations\n",
    "        Nsites: the number of sites in the system\n",
    "\n",
    "    Returns:\n",
    "        1D array The occupation for each site\n",
    "    \"\"\"\n",
    "    nocc_t = np.zeros(Nsites)\n",
    "    Nobs = 0\n",
    "    for occ_string, observations in counts.items():\n",
    "        # get the occupation of each site in the observed configuration\n",
    "        nocc = np.zeros(Nsites)\n",
    "\n",
    "        for site_index in range(Nsites):\n",
    "            nocc[site_index] = int(occ_string[site_index]) + int(occ_string[site_index + Nsites])\n",
    "\n",
    "        # now add the value to the correct time step and with it by the number of observations\n",
    "        Nobs = Nobs + observations\n",
    "        nocc_t = nocc_t + nocc * observations\n",
    "    return nocc_t / Nobs\n",
    "\n",
    "\n",
    "def squeezed_spin_corr_from_counts(counts: dict[str, int], Nsites: int) -> np.ndarray:\n",
    "    \"\"\"\n",
    "    The calculates the spin projection from all observed configurations.\n",
    "    However, we filter out situations with a hole and a doublon in the system.\n",
    "\n",
    "    Args:\n",
    "        counts: the counts of the observed configurations\n",
    "        Nsites: the number of sites in the system\n",
    "\n",
    "    Returns:\n",
    "        1D array The occupation for each site\n",
    "    \"\"\"\n",
    "    sz_squeeze = np.zeros(Nsites)\n",
    "    szvar = np.zeros(Nsites)\n",
    "    Nobs = 0\n",
    "    Nobs_var = np.zeros(Nsites)\n",
    "    Nobs = np.zeros(Nsites)\n",
    "\n",
    "    for occ_string, observations in counts.items():\n",
    "        occ_int = [int(occ) for occ in occ_string]\n",
    "        # get the spin of each site in the observed configuration\n",
    "        sz = np.zeros(Nsites)\n",
    "        nocc = np.zeros(Nsites)\n",
    "        for site_index in range(Nsites):\n",
    "            sz[site_index] = 1 / 2 * (occ_int[site_index] - occ_int[site_index + Nsites])\n",
    "            nocc[site_index] = occ_int[site_index] + occ_int[site_index + Nsites]\n",
    "\n",
    "        # now obtain the symmetrized variance\n",
    "        # first for the first site\n",
    "        site_index = 0\n",
    "        if nocc[site_index] == 1 and nocc[site_index + 1] == 1:\n",
    "                    szvar[site_index] = (\n",
    "                        szvar[site_index] + sz[site_index] * sz[site_index + 1] * observations\n",
    "                    )\n",
    "                    Nobs_var[site_index] = Nobs_var[site_index] + observations\n",
    "        \n",
    "        # then for the last site\n",
    "        site_index = Nsites - 1\n",
    "        if nocc[site_index] == 1 and nocc[site_index - 1] == 1:\n",
    "                    szvar[site_index] = (\n",
    "                        szvar[site_index] + sz[site_index] * sz[site_index - 1] * observations\n",
    "                    )\n",
    "                    Nobs_var[site_index] = Nobs_var[site_index] + observations\n",
    "        \n",
    "        # and finally for all the other sites\n",
    "        for site_index in range(1,Nsites-1):\n",
    "            if (\n",
    "                    nocc[site_index] == 1\n",
    "                    and nocc[site_index + 1] == 1\n",
    "                    and nocc[site_index - 1] == 1\n",
    "                ):\n",
    "                    szvar[site_index] = (\n",
    "                        szvar[site_index]\n",
    "                        + 1\n",
    "                        / 2\n",
    "                        * sz[site_index]\n",
    "                        * (sz[site_index + 1] + sz[site_index - 1])\n",
    "                        * observations\n",
    "                    )\n",
    "                    Nobs_var[site_index] = Nobs_var[site_index] + observations\n",
    "            elif (\n",
    "                    nocc[site_index] == 1\n",
    "                    and nocc[site_index + 1] == 1\n",
    "                    and nocc[site_index - 1] == 0\n",
    "                ):\n",
    "                    szvar[site_index] = (\n",
    "                        szvar[site_index] + sz[site_index] * sz[site_index + 1] * observations\n",
    "                    )\n",
    "                    Nobs_var[site_index] = Nobs_var[site_index] + observations\n",
    "            elif (\n",
    "                    nocc[site_index] == 1\n",
    "                    and nocc[site_index + 1] == 0\n",
    "                    and nocc[site_index - 1] == 1\n",
    "                ):\n",
    "                    szvar[site_index] = (\n",
    "                        szvar[site_index] + sz[site_index] * sz[site_index - 1] * observations\n",
    "                    )\n",
    "                    Nobs_var[site_index] = Nobs_var[site_index] + observations\n",
    "\n",
    "        sz_squeeze += sz * (nocc == 1) * observations\n",
    "        Nobs += (nocc == 1) * observations\n",
    "\n",
    "    # regularize the observations\n",
    "    var_weight = 1 / Nobs_var\n",
    "    var_weight[np.isinf(var_weight)] = 0\n",
    "\n",
    "    mean_weight = 1 / Nobs\n",
    "    mean_weight[np.isinf(mean_weight)] = 0\n",
    "    # now that we have gone through all observations, we can average\n",
    "    fluc = np.zeros(Nsites)\n",
    "    for site_index in range(Nsites):\n",
    "        if site_index == 0:\n",
    "            fluc[site_index] = 4 * (\n",
    "                szvar[site_index] * var_weight[site_index]\n",
    "                - sz_squeeze[site_index]\n",
    "                * mean_weight[site_index]\n",
    "                * sz_squeeze[site_index + 1]\n",
    "                * mean_weight[site_index + 1]\n",
    "            )\n",
    "        elif site_index == Nsites - 1:\n",
    "            fluc[site_index] = 4 * (\n",
    "                szvar[site_index] * var_weight[site_index]\n",
    "                - sz_squeeze[site_index]\n",
    "                * mean_weight[site_index]\n",
    "                * sz_squeeze[site_index - 1]\n",
    "                * mean_weight[site_index - 1]\n",
    "            )\n",
    "        else:\n",
    "            fluc[site_index] = 4 * (\n",
    "                szvar[site_index] * var_weight[site_index]\n",
    "                - sz_squeeze[site_index]\n",
    "                * mean_weight[site_index]\n",
    "                * (\n",
    "                    sz_squeeze[site_index + 1] * mean_weight[site_index + 1]\n",
    "                    + sz_squeeze[site_index - 1] * mean_weight[site_index - 1]\n",
    "                )\n",
    "                / 2\n",
    "            )\n",
    "    return fluc, Nobs_var, Nobs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us just visualize the initial state:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/_v/vw5sh5fd4lv4x0cnk8cxz__w0000gn/T/ipykernel_82435/3602395813.py:110: RuntimeWarning: divide by zero encountered in divide\n",
      "  var_weight = 1 / Nobs_var\n",
      "/var/folders/_v/vw5sh5fd4lv4x0cnk8cxz__w0000gn/T/ipykernel_82435/3602395813.py:113: RuntimeWarning: divide by zero encountered in divide\n",
      "  mean_weight = 1 / Nobs\n"
     ]
    }
   ],
   "source": [
    "noccs_init = np.zeros( Nsites)\n",
    "corrs_init = np.zeros(Nsites)\n",
    "\n",
    "counts_1 = job_up_init.result().get_counts()\n",
    "counts_2 = job_down_init.result().get_counts()\n",
    "total_counts = counts_1\n",
    "total_counts.update(counts_2)\n",
    "\n",
    "nocc = nocc_from_counts(total_counts, Nsites)\n",
    "corr, _, _ = squeezed_spin_corr_from_counts(total_counts, Nsites)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'site index')"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAGwCAYAAACkfh/eAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA5X0lEQVR4nO3de3wU9b3/8fcm5EICScBALrIQLgqkDZcGiFGsCCkJWEu8FTxoABGOFLAIItAKqFQjoBYRKmq5HqFQLRcvbSQnEvTUGDBAFQgoFAiXXMA0WRJKAtn5/ZGfq1sCbCDLbobX8/GYB+zMd77zma+r+3Z25rsWwzAMAQAAmIiPpwsAAABoaAQcAABgOgQcAABgOgQcAABgOgQcAABgOgQcAABgOgQcAABgOk08XYAn2O12nThxQs2bN5fFYvF0OQAAwAWGYej06dOKjo6Wj8+lr9FclwHnxIkTslqtni4DAABcgaNHj6pNmzaXbHNdBpzmzZtLqh2gkJAQD1cDAABcYbPZZLVaHZ/jl3JdBpzvvpYKCQkh4AAA0Mi4cnsJNxkDAADTIeAAAADTIeAAAADTIeAAAADTIeAAAADTIeAAAADTIeAAAADTIeAAAADTIeAAAADTuS5nMnaXmhrp00+lwkIpKkq6/XbJ19fTVXknxsp1jBXcgfeV6xir+vGa8TLcaOvWrcbPf/5zIyoqypBkbNiw4bL7bNmyxejZs6fh7+9vdOzY0Vi+fPkFbRYtWmS0a9fOCAgIMPr06WPk5ubWq67y8nJDklFeXl6v/S7lL38xjDZtDEP6fmnTpnY9nDFWrmOs4A68r1zHWNWPu8erPp/fbg04f/3rX43f/va3xvr1610KOP/85z+NoKAgY/LkycbevXuN1157zfD19TUyMjIcbdauXWv4+/sby5YtM/bs2WOMGTPGCAsLM4qLi12uq6EDzl/+YhgWi/M/UKl2ncXCvwg/xFi5jrGCO/C+ch1jVT/XYrzq8/ltMQzDuBZXiiwWizZs2KDU1NSLtpk2bZo+/PBD7d6927Fu2LBhKisrU0ZGhiQpISFBvXv31qJFiyRJdrtdVqtVEydO1PTp012qxWazKTQ0VOXl5Vf9Y5s1NVJMjHTsWN3bLRapTRvp0CEuaTJWrmOs4A68r1zHWNXPtRqv+nx+e9VNxjk5OUpKSnJal5ycrJycHElSdXW18vLynNr4+PgoKSnJ0aYuVVVVstlsTktD+fTTi/8DlWrz69Gjte2ud4yV6xgruAPvK9cxVvXjjePlVQGnqKhIERERTusiIiJks9n073//W6dOnVJNTU2dbYqKii7ab3p6ukJDQx2L1WptsJoLCxu2nZkxVq5jrOAOvK9cx1jVjzeOl1cFHHeZMWOGysvLHcvRo0cbrO+oqIZtZ2aMlesYK7gD7yvXMVb1443j5VUBJzIyUsXFxU7riouLFRISoqZNmyo8PFy+vr51tomMjLxovwEBAQoJCXFaGsrtt9d+r2ix1L3dYpGs1tp21zvGynWMFdyB95XrGKv68cbx8qqAk5iYqKysLKd1mZmZSkxMlCT5+/srPj7eqY3dbldWVpajzbXm6yu9+mrt3//zH+x3rxcs4CY0ibGqD8YK7sD7ynWMVf145Xhd/UNbF3f69Glj586dxs6dOw1JxiuvvGLs3LnTOHLkiGEYhjF9+nTj4YcfdrT/7jHxqVOnGvn5+cbixYvrfEw8ICDAWLFihbF3715j7NixRlhYmFFUVORyXddqHhyrlccI68JYuY6xgjvwvnIdY1U/7h4vr3lMPDs7W3feeecF60eMGKEVK1Zo5MiROnz4sLKzs532eeKJJ7R37161adNGM2fO1MiRI532X7RokebPn6+ioiL16NFDCxcuVEJCgst1NeRj4j/kNbM3NgKMlesYK7gD7yvXMVb1487xqs/n9zWbB8ebuCvgAAAA92m08+AAAAA0BAIOAAAwHQIOAAAwHQIOAAAwHQIOAAAwHQIOAAAwHQIOAAAwHQIOAAAwHQIOAAAwHQIOAAAwHQIOAAAwHQIOAAAwHQIOAAAwHQIOAAAwHQIOAAAwHQIOAAAwHQIOAAAwHQIOAAAwHQIOAAAwHQIOAAAwHQIOAAAwHQIOAAAwHQIOAAAwHQIOAAAwHQIOAAAwHQIOAAAwHQIOAAAwHQIOAAAwHQIOAAAwHQIOAAAwHQIOAAAwnWsScBYvXqyYmBgFBgYqISFB27Ztu2jbfv36yWKxXLDcddddjjYjR468YHtKSsq1OBUAANAINHH3AdatW6fJkydryZIlSkhI0IIFC5ScnKz9+/erdevWF7Rfv369qqurHa+//fZbde/eXQ888IBTu5SUFC1fvtzxOiAgwH0nAQAAGhW3X8F55ZVXNGbMGI0aNUqxsbFasmSJgoKCtGzZsjrbt2zZUpGRkY4lMzNTQUFBFwScgIAAp3YtWrRw96kAAIBGwq0Bp7q6Wnl5eUpKSvr+gD4+SkpKUk5Ojkt9LF26VMOGDVNwcLDT+uzsbLVu3VqdO3fWuHHj9O233160j6qqKtlsNqcFAACYl1sDzqlTp1RTU6OIiAin9RERESoqKrrs/tu2bdPu3bv16KOPOq1PSUnRqlWrlJWVpblz52rr1q0aNGiQampq6uwnPT1doaGhjsVqtV75SQEAAK/n9ntwrsbSpUsVFxenPn36OK0fNmyY4+9xcXHq1q2bOnbsqOzsbA0YMOCCfmbMmKHJkyc7XttsNkIOAAAm5tYrOOHh4fL19VVxcbHT+uLiYkVGRl5y38rKSq1du1ajR4++7HE6dOig8PBwHThwoM7tAQEBCgkJcVoAAIB5uTXg+Pv7Kz4+XllZWY51drtdWVlZSkxMvOS+77zzjqqqqvTQQw9d9jjHjh3Tt99+q6ioqKuuGQAANH5uf4pq8uTJeuutt7Ry5Url5+dr3Lhxqqys1KhRoyRJaWlpmjFjxgX7LV26VKmpqbrhhhuc1ldUVGjq1Kn6/PPPdfjwYWVlZWnIkCHq1KmTkpOT3X06AACgEXD7PThDhw7VyZMnNWvWLBUVFalHjx7KyMhw3HhcUFAgHx/nnLV//3793//9nzZv3nxBf76+vvryyy+1cuVKlZWVKTo6WgMHDtScOXOYCwcAAEiSLIZhGJ4u4lqz2WwKDQ1VeXk59+MAANBI1Ofzm9+iAgAApkPAAQAApkPAAQAApkPAAQAApkPAAQAApkPAAQAApkPAAQAApkPAAQAApkPAAQAApkPAAQAApkPAAQAApkPAAQAApkPAAQAApkPAAQAApkPAAQAApkPAAQAApkPAAQAApkPAAQAApkPAAQAApkPAAQAApkPAAQAApkPAAQAApkPAAQAApkPAAQAApkPAAQAApkPAAQAApkPAAQAApkPAAQAApkPAAQAApkPAAQAApkPAAQAApnNNAs7ixYsVExOjwMBAJSQkaNu2bRdtu2LFClksFqclMDDQqY1hGJo1a5aioqLUtGlTJSUl6ZtvvnH3aQAAgEbC7QFn3bp1mjx5smbPnq0dO3aoe/fuSk5OVklJyUX3CQkJUWFhoWM5cuSI0/Z58+Zp4cKFWrJkiXJzcxUcHKzk5GSdPXvW3acDAAAaAbcHnFdeeUVjxozRqFGjFBsbqyVLligoKEjLli276D4Wi0WRkZGOJSIiwrHNMAwtWLBATz/9tIYMGaJu3bpp1apVOnHihDZu3Oju0wEAAI2AWwNOdXW18vLylJSU9P0BfXyUlJSknJyci+5XUVGhdu3ayWq1asiQIdqzZ49j26FDh1RUVOTUZ2hoqBISEi7aZ1VVlWw2m9MCAADMy60B59SpU6qpqXG6AiNJERERKioqqnOfzp07a9myZdq0aZPefvtt2e123XrrrTp27JgkOfarT5/p6ekKDQ11LFar9WpPDQAAeDGve4oqMTFRaWlp6tGjh+644w6tX79erVq10htvvHHFfc6YMUPl5eWO5ejRow1YMQAA8DZuDTjh4eHy9fVVcXGx0/ri4mJFRka61Iefn5969uypAwcOSJJjv/r0GRAQoJCQEKcFAACYl1sDjr+/v+Lj45WVleVYZ7fblZWVpcTERJf6qKmp0VdffaWoqChJUvv27RUZGenUp81mU25urst9AgAAc2vi7gNMnjxZI0aMUK9evdSnTx8tWLBAlZWVGjVqlCQpLS1NN954o9LT0yVJzz33nG655RZ16tRJZWVlmj9/vo4cOaJHH31UUu0TVpMmTdLvfvc73XTTTWrfvr1mzpyp6Ohopaamuvt0AABAI+D2gDN06FCdPHlSs2bNUlFRkXr06KGMjAzHTcIFBQXy8fn+QtK//vUvjRkzRkVFRWrRooXi4+P12WefKTY21tHmqaeeUmVlpcaOHauysjL17dtXGRkZF0wICAAArk8WwzAMTxdxrdlsNoWGhqq8vJz7cQAAaCTq8/ntdU9RAQAAXC0CDgAAMB0CDgAAMB0CDgAAMB0CDgAAMB0CDgAAMB0CDgAAMB0CDgAAMB0CDgAAMB0CDgAAMB0CDgAAMB0CDgAAMB0CDgAAMB0CDgAAMB0CDgAAMB0CDgAAMB0CDgAAMB0CDgAAMB0CDgAAMB0CDgAAMB0CDgAAMB0CDgAAMB0CDgAAMB0CDgAAMB0CDgAAMB0CDgAAMB0CDgAAMB0CDgAAMB0CDgAAMB0CDgAAMB0CDgAAMJ1rEnAWL16smJgYBQYGKiEhQdu2bbto27feeku33367WrRooRYtWigpKemC9iNHjpTFYnFaUlJS3H0aAACgkXB7wFm3bp0mT56s2bNna8eOHerevbuSk5NVUlJSZ/vs7Gw9+OCD2rJli3JycmS1WjVw4EAdP37cqV1KSooKCwsdy5/+9Cd3nwoAAGgkLIZhGO48QEJCgnr37q1FixZJkux2u6xWqyZOnKjp06dfdv+amhq1aNFCixYtUlpamqTaKzhlZWXauHHjFdVks9kUGhqq8vJyhYSEXFEfAADg2qrP57dbr+BUV1crLy9PSUlJ3x/Qx0dJSUnKyclxqY8zZ87o3LlzatmypdP67OxstW7dWp07d9a4ceP07bffXrSPqqoq2Ww2pwUAAJiXWwPOqVOnVFNTo4iICKf1ERERKioqcqmPadOmKTo62ikkpaSkaNWqVcrKytLcuXO1detWDRo0SDU1NXX2kZ6ertDQUMditVqv/KQAAIDXa+LpAi7lxRdf1Nq1a5Wdna3AwEDH+mHDhjn+HhcXp27duqljx47Kzs7WgAEDLuhnxowZmjx5suO1zWYj5AAAYGJuvYITHh4uX19fFRcXO60vLi5WZGTkJfd96aWX9OKLL2rz5s3q1q3bJdt26NBB4eHhOnDgQJ3bAwICFBIS4rQAAADzcmvA8ff3V3x8vLKyshzr7Ha7srKylJiYeNH95s2bpzlz5igjI0O9evW67HGOHTumb7/9VlFRUQ1SNwAAaNzc/pj45MmT9dZbb2nlypXKz8/XuHHjVFlZqVGjRkmS0tLSNGPGDEf7uXPnaubMmVq2bJliYmJUVFSkoqIiVVRUSJIqKio0depUff755zp8+LCysrI0ZMgQderUScnJye4+HQAA0Ai4/R6coUOH6uTJk5o1a5aKiorUo0cPZWRkOG48LigokI/P9znr9ddfV3V1te6//36nfmbPnq1nnnlGvr6++vLLL7Vy5UqVlZUpOjpaAwcO1Jw5cxQQEODu0wEAAI2A2+fB8UbMgwMAQOPjNfPgAAAAeAIBBwAAmA4BBwAAmA4BBwAAmA4BBwAAmA4BBwAAmA4BBwAAmA4BBwAAmA4BBwAAmA4BBwAAmA4BBwAAmI7bf2wTAIDrid1uV3V1tafLaJT8/Pzk6+vbIH0RcAAAaCDV1dU6dOiQ7Ha7p0tptMLCwhQZGSmLxXJV/RBwAABoAIZhqLCwUL6+vrJarfLx4S6Q+jAMQ2fOnFFJSYkkKSoq6qr6I+AAANAAzp8/rzNnzig6OlpBQUGeLqdRatq0qSSppKRErVu3vqqvq4iXAAA0gJqaGkmSv7+/hytp3L4Lh+fOnbuqfgg4AAA0oKu9d+R611DjR8ABAACmQ8ABAACmw03GAAB4kZoa6dNPpcJCKSpKuv12qYGmhrmuEHAAAPAS69dLv/61dOzY9+vatJFefVW6917P1dUY8RUVAABeYP166f77ncONJB0/Xrt+/Xr3HLdfv356/PHH9dRTT6lly5aKjIzUM88849heUFCgIUOGqFmzZgoJCdEvf/lLFRcXO/Xx/vvvq3fv3goMDFR4eLjuuecex7aqqipNmzZNVqtVAQEB6tSpk5YuXeqek/kBAg4AAB5WU1N75cYwLtz23bpJk2rbucPKlSsVHBys3NxczZs3T88995wyMzNlt9s1ZMgQlZaWauvWrcrMzNQ///lPDR061LHvhx9+qHvuuUeDBw/Wzp07lZWVpT59+ji2p6Wl6U9/+pMWLlyo/Px8vfHGG2rWrJl7TuQHLIZR13Cam81mU2hoqMrLyxUSEuLpcgAAJnD27FkdOnRI7du3V2BgYL32zc6W7rzz8u22bJH69bui8i6qX79+qqmp0aeffupY16dPH/Xv318DBgzQoEGDdOjQIVmtVknS3r179aMf/Ujbtm1T7969deutt6pDhw56++23L+j766+/VufOnZWZmamkpCSX6rnUONbn85srOAAAeFhhYcO2q69u3bo5vY6KilJJSYny8/NltVod4UaSYmNjFRYWpvz8fEnSrl27NGDAgDr73bVrl3x9fXXHHXe4p/BLIOAAAOBhrv7s0lX+PNNF+fn5Ob22WCwu/2Dodz+vUN9t7kbAAQDAw26/vfZpqYtN4muxSFZrbbtrqWvXrjp69KiOHj3qWLd3716VlZUpNjZWUu3Vn6ysrDr3j4uLk91u19atW69JvT9EwAEAwMN8fWsfBZcuDDnfvV6w4NrPh5OUlKS4uDgNHz5cO3bs0LZt25SWlqY77rhDvXr1kiTNnj1bf/rTnzR79mzl5+frq6++0ty5cyVJMTExGjFihB555BFt3LhRhw4dUnZ2tv785z+7vXYCDgAAXuDee6V335VuvNF5fZs2tes9MQ+OxWLRpk2b1KJFC/30pz9VUlKSOnTooHXr1jna9OvXT++8847ee+899ejRQ/3799e2bdsc219//XXdf//9+tWvfqUuXbpozJgxqqysdH/tPEXFU1QAgKt3NU9R/dD1PpNxQz1FxUzGAEyj5pxd2/9SoNIjp9WyXXP1vq+tfP24UI3Gxde34R8Fvx5dk3/zFy9erJiYGAUGBiohIcHp0lVd3nnnHXXp0kWBgYGKi4vTX//6V6fthmFo1qxZioqKUtOmTZWUlKRvvvnGnacAwMtlLszX/PAXtf/BWaqYPkf7H5yl+eEvKnNhvqdLA+ABbg8469at0+TJkzV79mzt2LFD3bt3V3JyskpKSups/9lnn+nBBx/U6NGjtXPnTqWmpio1NVW7d+92tJk3b54WLlyoJUuWKDc3V8HBwUpOTtbZs2fdfToAvFDmwnwd+PVCdbTt1LcK13511rcKV0fbTh349UJCDnAdcvs9OAkJCerdu7cWLVokSbLb7bJarZo4caKmT59+QfuhQ4eqsrJSH3zwgWPdLbfcoh49emjJkiUyDEPR0dGaMmWKnnzySUlSeXm5IiIitGLFCg0bNuyyNXEPDmAeNefsmh/+ojradmqPYiX98BEUQz/SXh0M/YmmnpzG11Vwq4a6B+d61yhmMq6urlZeXp7T9Mw+Pj5KSkpSTk5Onfvk5ORcMJ1zcnKyo/2hQ4dUVFTk1CY0NFQJCQkX7bOqqko2m81pAWAO2/9SoCjbPh2VVc7hRpIsOqo2iirP1/a/FHiiPAAe4taAc+rUKdXU1CgiIsJpfUREhIqKiurcp6io6JLtv/uzPn2mp6crNDTUsfxwymkAjVvpkdNqqrOqVHCd2ysVrECdVemR09e4MgCedF1cr50xY4bKy8sdyw9nZATQuLVs11z/VqCCVfe8GsGq1FkFqmW75te4MgCe5NaAEx4eLl9fXxUXFzutLy4uVmRkZJ37REZGXrL9d3/Wp8+AgACFhIQ4LQDMofd9bVUY0kVWHZX0n7cUGrLqmApDu6r3fW09UR4AD3FrwPH391d8fLzTb1TY7XZlZWUpMTGxzn0SExMv+E2LzMxMR/v27dsrMjLSqY3NZlNubu5F+wRgXr5+Poqfc49OKVw/0l6FqFy+Oq8QletH2qtTClf8c6ncYAxcZ9w+0d/kyZM1YsQI9erVS3369NGCBQtUWVmpUaNGSZLS0tJ04403Kj09XZL061//WnfccYdefvll3XXXXVq7dq2++OILvfnmm5Jqp42eNGmSfve73+mmm25S+/btNXPmTEVHRys1NdXdpwPAC/3s8a6SHlfezA2Ksu3TjTquswrUwdCfKP651P+/HcD1xO0BZ+jQoTp58qRmzZqloqIi9ejRQxkZGY6bhAsKCuTj8/3/Wd16661as2aNnn76af3mN7/RTTfdpI0bN+rHP/6xo81TTz2lyspKjR07VmVlZerbt68yMjJ4LA+4jv3s8a7qP66zYybjdu2a6yFmMkZjZLdLBQXS6dNS8+ZS27aSD+/j+uK3qLgfBwDQABpkHpz8fGnDBmnfPunsWSkwUOrSRbrnHqlr47gSWV1dLX9//wvWnzt3Tn5+fpfdv1HMgwMAAFyUny8tXCjt3CmFh0udO9f+uXNn7fp8983IbbfbNW/ePHXq1EkBAQFq27atnn/+eUnSV199pf79+6tp06a64YYbNHbsWFVUVDj2HTlypFJTU/X8888rOjpanTt31uHDh2WxWLRu3TrdcccdCgwM1OrVq91Wf134sU0AADzNbq+9cnPqlBQbK1n+/6SVISG1r/fulTZurA09bvi6asaMGXrrrbf0+9//Xn379lVhYaH27dunyspKJScnKzExUdu3b1dJSYkeffRRTZgwQStWrHDsn5WVpZCQEGVmZjr1O336dL388svq2bPnNb+NhIADAICnFRTUfi1ltX4fbr5jsUht2tRewSkokGJiGvTQp0+f1quvvqpFixZpxIgRkqSOHTuqb9++euutt3T27FmtWrVKwcG1k2kuWrRId999t+bOneu4nzY4OFh//OMfHV9NHT58WJI0adIk3XvvvQ1ar6v4igoAAE87fbr2npvgumfkVnBw7fbTDT8jd35+vqqqqjRgwIA6t3Xv3t0RbiTptttuk91u1/79+x3r4uLi6rzvplevXg1er6sIOAAAeFrz5rU3FFfWPSO3Kitrtzdv+Bm5mzZtetV9BF8kmF1s/bVAwAEAwNPatq19WuroUek/H242DOnYsdqnqNo2/IzcN910k5o2bXrBJLuS1LVrV/3jH/9Q5Q+C19///nf5+Pioc+fODV5LQyLgAADgaT4+tY+Ch4fX3lBcXi6dP1/75969tetTU91yg3FgYKCmTZump556SqtWrdLBgwf1+eefa+nSpRo+fLgCAwM1YsQI7d69W1u2bNHEiRP18MMPX/Cj196Gm4wBAPAGXbtKjz/+/Tw4x4/Xfi31k5/Uhhs3zoMzc+ZMNWnSRLNmzdKJEycUFRWlxx57TEFBQfroo4/061//Wr1791ZQUJDuu+8+vfLKK26rpaEw0R8T/QEAGkCDTPQnXfczGTfURH9cwQEAwJv4+DT4o+DXo+snEgIAgOsGAQcAAJgOAQcAAJgOAQcAgAZ0HT6706AaavwIOAAANABfX19JUnV1tYcradzOnDkjSfLz87uqfniKCgCABtCkSRMFBQXp5MmT8vPzk8919Gh3QzAMQ2fOnFFJSYnCwsIcgfFKEXAAAGgAFotFUVFROnTokI4cOeLpchqtsLAwRUZGXnU/BBwAABqIv7+/brrpJr6mukJ+fn5XfeXmOwQcAAAakI+Pz9XNZIwGwReEAADAdAg4AADAdAg4AADAdAg4AADAdAg4AADAdAg4AADAdAg4AADAdAg4AADAdAg4AADAdAg4AADAdAg4AADAdAg4AADAdNwacEpLSzV8+HCFhIQoLCxMo0ePVkVFxSXbT5w4UZ07d1bTpk3Vtm1bPf744yovL3dqZ7FYLljWrl3rzlMBAACNiFt/TXz48OEqLCxUZmamzp07p1GjRmns2LFas2ZNne1PnDihEydO6KWXXlJsbKyOHDmixx57TCdOnNC7777r1Hb58uVKSUlxvA4LC3PnqQAAgEbEYhiG4Y6O8/PzFRsbq+3bt6tXr16SpIyMDA0ePFjHjh1TdHS0S/288847euihh1RZWakmTWrzmMVi0YYNG5SamnpFtdlsNoWGhqq8vFwhISFX1AcAALi26vP57bavqHJychQWFuYIN5KUlJQkHx8f5ebmutzPdyfxXbj5zvjx4xUeHq4+ffpo2bJlulROq6qqks1mc1oAAIB5ue0rqqKiIrVu3dr5YE2aqGXLlioqKnKpj1OnTmnOnDkaO3as0/rnnntO/fv3V1BQkDZv3qxf/epXqqio0OOPP15nP+np6Xr22Wev7EQAAECjU+8rONOnT6/zJt8fLvv27bvqwmw2m+666y7FxsbqmWeecdo2c+ZM3XbbberZs6emTZump556SvPnz79oXzNmzFB5ebljOXr06FXXBwAAvFe9r+BMmTJFI0eOvGSbDh06KDIyUiUlJU7rz58/r9LSUkVGRl5y/9OnTyslJUXNmzfXhg0b5Ofnd8n2CQkJmjNnjqqqqhQQEHDB9oCAgDrXAwAAc6p3wGnVqpVatWp12XaJiYkqKytTXl6e4uPjJUkff/yx7Ha7EhISLrqfzWZTcnKyAgIC9N577ykwMPCyx9q1a5datGhBiAEAAJLceA9O165dlZKSojFjxmjJkiU6d+6cJkyYoGHDhjmeoDp+/LgGDBigVatWqU+fPrLZbBo4cKDOnDmjt99+2+mG4FatWsnX11fvv/++iouLdcsttygwMFCZmZl64YUX9OSTT7rrVAAAQCPj1nlwVq9erQkTJmjAgAHy8fHRfffdp4ULFzq2nzt3Tvv379eZM2ckSTt27HA8YdWpUyenvg4dOqSYmBj5+flp8eLFeuKJJ2QYhjp16qRXXnlFY8aMceepAACARsRt8+B4M+bBAQCg8fGKeXAAAAA8hYADAABMh4ADAABMh4ADAABMh4ADAABMh4ADAABMh4ADAABMh4ADAABMh4ADAABMh4ADAABMh4ADAABMh4ADAABMh4ADAABMh4ADAABMh4ADAABMh4ADAABMh4ADAABMh4ADAABMh4ADAABMh4ADAABMh4ADAABMh4ADAABMh4ADAABMh4ADAABMh4ADAABMh4ADAABMh4ADAABMh4ADAABMh4ADAABMh4ADAABMh4ADAABMx60Bp7S0VMOHD1dISIjCwsI0evRoVVRUXHKffv36yWKxOC2PPfaYU5uCggLdddddCgoKUuvWrTV16lSdP3/enacCAAAakSbu7Hz48OEqLCxUZmamzp07p1GjRmns2LFas2bNJfcbM2aMnnvuOcfroKAgx99ramp01113KTIyUp999pkKCwuVlpYmPz8/vfDCC247FwAA0HhYDMMw3NFxfn6+YmNjtX37dvXq1UuSlJGRocGDB+vYsWOKjo6uc79+/fqpR48eWrBgQZ3b//a3v+nnP/+5Tpw4oYiICEnSkiVLNG3aNJ08eVL+/v6Xrc1msyk0NFTl5eUKCQm5shMEAADXVH0+v932FVVOTo7CwsIc4UaSkpKS5OPjo9zc3Evuu3r1aoWHh+vHP/6xZsyYoTNnzjj1GxcX5wg3kpScnCybzaY9e/bU2V9VVZVsNpvTAgAAzMttX1EVFRWpdevWzgdr0kQtW7ZUUVHRRff7r//6L7Vr107R0dH68ssvNW3aNO3fv1/r16939PvDcCPJ8fpi/aanp+vZZ5+9mtMBAACNSL0DzvTp0zV37txLtsnPz7/igsaOHev4e1xcnKKiojRgwAAdPHhQHTt2vKI+Z8yYocmTJzte22w2Wa3WK64RAAB4t3oHnClTpmjkyJGXbNOhQwdFRkaqpKTEaf358+dVWlqqyMhIl4+XkJAgSTpw4IA6duyoyMhIbdu2zalNcXGxJF2034CAAAUEBLh8TAAA0LjVO+C0atVKrVq1umy7xMRElZWVKS8vT/Hx8ZKkjz/+WHa73RFaXLFr1y5JUlRUlKPf559/XiUlJY6vwDIzMxUSEqLY2Nh6ng0AADAjt91k3LVrV6WkpGjMmDHatm2b/v73v2vChAkaNmyY4wmq48ePq0uXLo4rMgcPHtScOXOUl5enw4cP67333lNaWpp++tOfqlu3bpKkgQMHKjY2Vg8//LD+8Y9/6KOPPtLTTz+t8ePHc5UGAABIcvNEf6tXr1aXLl00YMAADR48WH379tWbb77p2H7u3Dnt37/f8ZSUv7+//vd//1cDBw5Uly5dNGXKFN133316//33Hfv4+vrqgw8+kK+vrxITE/XQQw8pLS3Nad4cAABwfXPbPDjejHlwAABofLxiHhwAAABPIeAAAADTIeAAAADTIeAAAADTIeAAAADTIeAAAADTIeAAAADTIeAAAADTIeAAAADTIeAAAADTIeAAAADTIeAAAADTIeAAAADTIeAAAADTIeAAAADTIeAAAADTIeAAAADTIeAAAADTIeAAAADTIeAAAADTIeAAAADTIeAAAADTIeAAAADTIeAAAADTIeAAAADTIeAAAADTIeAAAADTIeAAAADTIeAAAADTIeAAAADTIeAAAADTcWvAKS0t1fDhwxUSEqKwsDCNHj1aFRUVF21/+PBhWSyWOpd33nnH0a6u7WvXrnXnqQAAgEakiTs7Hz58uAoLC5WZmalz585p1KhRGjt2rNasWVNne6vVqsLCQqd1b775pubPn69BgwY5rV++fLlSUlIcr8PCwhq8fgAA0Di5LeDk5+crIyND27dvV69evSRJr732mgYPHqyXXnpJ0dHRF+zj6+uryMhIp3UbNmzQL3/5SzVr1sxpfVhY2AVtAQAAJDd+RZWTk6OwsDBHuJGkpKQk+fj4KDc316U+8vLytGvXLo0ePfqCbePHj1d4eLj69OmjZcuWyTCMi/ZTVVUlm83mtAAAAPNy2xWcoqIitW7d2vlgTZqoZcuWKioqcqmPpUuXqmvXrrr11lud1j/33HPq37+/goKCtHnzZv3qV79SRUWFHn/88Tr7SU9P17PPPntlJwIAABqdel/BmT59+kVvBP5u2bdv31UX9u9//1tr1qyp8+rNzJkzddttt6lnz56aNm2annrqKc2fP/+ifc2YMUPl5eWO5ejRo1ddHwAA8F71voIzZcoUjRw58pJtOnTooMjISJWUlDitP3/+vEpLS126d+bdd9/VmTNnlJaWdtm2CQkJmjNnjqqqqhQQEHDB9oCAgDrXAwAAc6p3wGnVqpVatWp12XaJiYkqKytTXl6e4uPjJUkff/yx7Ha7EhISLrv/0qVL9Ytf/MKlY+3atUstWrQgxAAAAEluvAena9euSklJ0ZgxY7RkyRKdO3dOEyZM0LBhwxxPUB0/flwDBgzQqlWr1KdPH8e+Bw4c0CeffKK//vWvF/T7/vvvq7i4WLfccosCAwOVmZmpF154QU8++aS7TgUAADQybp0HZ/Xq1ZowYYIGDBggHx8f3XfffVq4cKFj+7lz57R//36dOXPGab9ly5apTZs2Gjhw4AV9+vn5afHixXriiSdkGIY6deqkV155RWPGjHHnqQAAgEbEYlzq+WqTstlsCg0NVXl5uUJCQjxdDgAAcEF9Pr/5LSoAAGA6BBwAAGA6BBwAAGA6BBwAAGA6BBwAAGA6BBwAAGA6BBwAAGA6BBwAAGA6BBwAAGA6BBwAAGA6BBwAAGA6BBwAAGA6BBwAAGA6BBwAAGA6BBwAAGA6BBwAAGA6BBwAAGA6BBwAAGA6BBwAAGA6BBwAAGA6BBwAAGA6BBwAAGA6BBwAAGA6BBwAAGA6BBwAAGA6BBwAAGA6BBwAAGA6BBwAAGA6BBwAAGA6BBwAAGA6TTxdgKnY7VJBgXT6tNS8udS2reRDhqwTY+U6xgruwPvKdYxV/XjJeLkt4Dz//PP68MMPtWvXLvn7+6usrOyy+xiGodmzZ+utt95SWVmZbrvtNr3++uu66aabHG1KS0s1ceJEvf/++/Lx8dF9992nV199Vc2aNXPXqbgmP1/asEHat086e1YKDJS6dJHuuUfq2tWztXkbxsp1jBXcgfeV6xir+vGi8XJbpKqurtYDDzygcePGubzPvHnztHDhQi1ZskS5ubkKDg5WcnKyzp4962gzfPhw7dmzR5mZmfrggw/0ySefaOzYse44Bdfl50sLF0o7d0rh4VLnzrV/7txZuz4/37P1eRPGynWMFdyB95XrGKv68bLxclvAefbZZ/XEE08oLi7OpfaGYWjBggV6+umnNWTIEHXr1k2rVq3SiRMntHHjRklSfn6+MjIy9Mc//lEJCQnq27evXnvtNa1du1YnTpxw16lcmt1em1ZPnZJiY6WQEMnXt/bP2Nja9Rs31ra73jFWrmOs4A68r1zHWNWPF46X13yJeOjQIRUVFSkpKcmxLjQ0VAkJCcrJyZEk5eTkKCwsTL169XK0SUpKko+Pj3Jzcy/ad1VVlWw2m9PSYAoKai/FWa2SxeK8zWKR2rSpTa0FBQ13zMaKsXIdYwV34H3lOsaqfrxwvLwm4BQVFUmSIiIinNZHREQ4thUVFal169ZO25s0aaKWLVs62tQlPT1doaGhjsVqtTZc4adP137PGBxc9/bg4Nrtp0833DEbK8bKdYwV3IH3lesYq/rxwvGqV8CZPn26LBbLJZd9+/a5q9YrNmPGDJWXlzuWo0ePNlznzZvX3kRVWVn39srK2u3NmzfcMRsrxsp1jBXcgfeV6xir+vHC8apXwJkyZYry8/MvuXTo0OGKComMjJQkFRcXO60vLi52bIuMjFRJSYnT9vPnz6u0tNTRpi4BAQEKCQlxWhpM27a1d4gfPSoZhvM2w5COHau9c7xt24Y7ZmPFWLmOsYI78L5yHWNVP144XvV6TLxVq1Zq1aqVWwpp3769IiMjlZWVpR49ekiSbDabcnNzHU9iJSYmqqysTHl5eYqPj5ckffzxx7Lb7UpISHBLXZfl41P7+NvRo9LevbXfMwYH16bVY8dq7yBPTWXOBImxqg/GCu7A+8p1jFX9eOF4WQzjP6NWwygoKFBpaanee+89zZ8/X59++qkkqVOnTo45a7p06aL09HTdc889kqS5c+fqxRdf1MqVK9W+fXvNnDlTX375pfbu3avAwEBJ0qBBg1RcXKwlS5bo3LlzGjVqlHr16qU1a9a4XJvNZlNoaKjKy8sb7mpOXc/+d+1a+w+UuRKcMVauY6zgDryvXMdY1Y+bx6s+n99uCzgjR47UypUrL1i/ZcsW9evXr/bgFouWL1+ukSNHSvp+or8333xTZWVl6tu3r/7whz/o5ptvduxfWlqqCRMmOE30t3DhwnpN9OeWgCN5zeyNjQJj5TrGCu7A+8p1jFX9uHG8vCLgeDO3BRwAAOA29fn8JoICAADTIeAAAADTIeAAAADTIeAAAADTIeAAAADTIeAAAADTIeAAAADTIeAAAADTIeAAAADTqdePbZrFd5M322w2D1cCAABc9d3ntis/wnBdBpzTp09LkqxWq4crAQAA9XX69GmFhoZess11+VtUdrtdJ06cUPPmzWWxWBq0b5vNJqvVqqNHj/I7V5fBWLmOsXIdY+U6xsp1jFX9uGu8DMPQ6dOnFR0dLZ/L/IDndXkFx8fHR23atHHrMUJCQviXwEWMlesYK9cxVq5jrFzHWNWPO8brclduvsNNxgAAwHQIOAAAwHQIOA0sICBAs2fPVkBAgKdL8XqMlesYK9cxVq5jrFzHWNWPN4zXdXmTMQAAMDeu4AAAANMh4AAAANMh4AAAANMh4AAAANMh4LjRL37xC7Vt21aBgYGKiorSww8/rBMnTni6LK9z+PBhjR49Wu3bt1fTpk3VsWNHzZ49W9XV1Z4uzSs9//zzuvXWWxUUFKSwsDBPl+N1Fi9erJiYGAUGBiohIUHbtm3zdEle55NPPtHdd9+t6OhoWSwWbdy40dMlea309HT17t1bzZs3V+vWrZWamqr9+/d7uiyv9Prrr6tbt26Oyf0SExP1t7/9zWP1EHDc6M4779Sf//xn7d+/X3/5y1908OBB3X///Z4uy+vs27dPdrtdb7zxhvbs2aPf//73WrJkiX7zm994ujSvVF1drQceeEDjxo3zdCleZ926dZo8ebJmz56tHTt2qHv37kpOTlZJSYmnS/MqlZWV6t69uxYvXuzpUrze1q1bNX78eH3++efKzMzUuXPnNHDgQFVWVnq6NK/Tpk0bvfjii8rLy9MXX3yh/v37a8iQIdqzZ49nCjJwzWzatMmwWCxGdXW1p0vxevPmzTPat2/v6TK82vLly43Q0FBPl+FV+vTpY4wfP97xuqamxoiOjjbS09M9WJV3k2Rs2LDB02U0GiUlJYYkY+vWrZ4upVFo0aKF8cc//tEjx+YKzjVSWlqq1atX69Zbb5Wfn5+ny/F65eXlatmypafLQCNSXV2tvLw8JSUlOdb5+PgoKSlJOTk5HqwMZlJeXi5J/PfpMmpqarR27VpVVlYqMTHRIzUQcNxs2rRpCg4O1g033KCCggJt2rTJ0yV5vQMHDui1117Tf//3f3u6FDQip06dUk1NjSIiIpzWR0REqKioyENVwUzsdrsmTZqk2267TT/+8Y89XY5X+uqrr9SsWTMFBAToscce04YNGxQbG+uRWgg49TR9+nRZLJZLLvv27XO0nzp1qnbu3KnNmzfL19dXaWlpMq6TyaPrO1aSdPz4caWkpOiBBx7QmDFjPFT5tXclYwXg2ho/frx2796ttWvXeroUr9W5c2ft2rVLubm5GjdunEaMGKG9e/d6pBZ+qqGeTp48qW+//faSbTp06CB/f/8L1h87dkxWq1WfffaZxy7ZXUv1HasTJ06oX79+uuWWW7RixQr5+Fw/+ftK3lcrVqzQpEmTVFZW5ubqGofq6moFBQXp3XffVWpqqmP9iBEjVFZWxtXTi7BYLNqwYYPTmOFCEyZM0KZNm/TJJ5+offv2ni6n0UhKSlLHjh31xhtvXPNjN7nmR2zkWrVqpVatWl3Rvna7XZJUVVXVkCV5rfqM1fHjx3XnnXcqPj5ey5cvv67CjXR17yvU8vf3V3x8vLKyshwf1na7XVlZWZowYYJni0OjZRiGJk6cqA0bNig7O5twU092u91jn3kEHDfJzc3V9u3b1bdvX7Vo0UIHDx7UzJkz1bFjx+vi6k19HD9+XP369VO7du300ksv6eTJk45tkZGRHqzMOxUUFKi0tFQFBQWqqanRrl27JEmdOnVSs2bNPFuch02ePFkjRoxQr1691KdPHy1YsECVlZUaNWqUp0vzKhUVFTpw4IDj9aFDh7Rr1y61bNlSbdu29WBl3mf8+PFas2aNNm3apObNmzvu5woNDVXTpk09XJ13mTFjhgYNGqS2bdvq9OnTWrNmjbKzs/XRRx95piCPPLt1Hfjyyy+NO++802jZsqUREBBgxMTEGI899phx7NgxT5fmdZYvX25IqnPBhUaMGFHnWG3ZssXTpXmF1157zWjbtq3h7+9v9OnTx/j88889XZLX2bJlS53voREjRni6NK9zsf82LV++3NOleZ1HHnnEaNeuneHv72+0atXKGDBggLF582aP1cM9OAAAwHSurxsdAADAdYGAAwAATIeAAwAATIeAAwAATIeAAwAATIeAAwAATIeAAwAATIeAAwAATIeAA+CaGTlypFt+1HHFihUKCwu76n769eunSZMmXXU/ADyP36ICcM28+uqr+uHk6f369VOPHj20YMGCq+p36NChGjx48FVWB8BMCDgArpnQ0FC39Nu0aVN++BCAE76iAtCg3n33XcXFxalp06a64YYblJSUpMrKSknOX1GNHDlSW7du1auvviqLxSKLxaLDhw9Lknbv3q1BgwapWbNmioiI0MMPP6xTp05d9Jj/+RXVM888ox49euh//ud/FBMTo9DQUA0bNkynT592tKmsrFRaWpqaNWumqKgovfzyyxf0W1VVpSeffFI33nijgoODlZCQoOzsbEnS2bNn9aMf/Uhjx451tD948KCaN2+uZcuWXeHoAWgoBBwADaawsFAPPvigHnnkEeXn5ys7O1v33nuv6vpN31dffVWJiYkaM2aMCgsLVVhYKKvVqrKyMvXv3189e/bUF198oYyMDBUXF+uXv/xlvWo5ePCgNm7cqA8++EAffPCBtm7dqhdffNGxferUqdq6das2bdqkzZs3Kzs7Wzt27HDqY8KECcrJydHatWv15Zdf6oEHHlBKSoq++eYbBQYGavXq1Vq5cqU2bdqkmpoaPfTQQ/rZz36mRx555MoGEECD4SsqAA2msLBQ58+f17333qt27dpJkuLi4upsGxoaKn9/fwUFBSkyMtKxftGiRerZs6deeOEFx7ply5bJarXq66+/1s033+xSLXa7XStWrFDz5s0lSQ8//LCysrL0/PPPq6KiQkuXLtXbb7+tAQMGSJJWrlypNm3aOPYvKCjQ8uXLVVBQoOjoaEnSk08+qYyMDC1fvlwvvPCCevTood/97nd69NFHNWzYMB05ckQffPBBPUYMgLsQcAA0mO7du2vAgAGKi4tTcnKyBg4cqPvvv18tWrRwuY9//OMf2rJli5o1a3bBtoMHD7occGJiYhzhRpKioqJUUlLi6Ke6uloJCQmO7S1btlTnzp0dr7/66ivV1NRccLyqqirdcMMNjtdTpkzRxo0btWjRIv3tb39z2gbAcwg4ABqMr6+vMjMz9dlnn2nz5s167bXX9Nvf/la5ublq3769S31UVFTo7rvv1ty5cy/YFhUV5XItfn5+Tq8tFovsdrvL+1dUVMjX11d5eXny9fV12vbD8FVSUqKvv/5avr6++uabb5SSkuLyMQC4D/fgAGhQFotFt912m5599lnt3LlT/v7+2rBhQ51t/f39VVNT47TuJz/5ifbs2aOYmBh16tTJaQkODm6QGjt27Cg/Pz/l5uY61v3rX//S119/7Xjds2dP1dTUqKSk5II6fviV2iOPPKK4uDitXLlS06ZNU35+foPUCODqEHAANJjc3Fy98MIL+uKLL1RQUKD169fr5MmT6tq1a53tY2JilJubq8OHD+vUqVOy2+0aP368SktL9eCDD2r79u06ePCgPvroI40aNeqCMHSlmjVrptGjR2vq1Kn6+OOPtXv3bo0cOVI+Pt//J/Hmm2/W8OHDlZaWpvXr1+vQoUPatm2b0tPT9eGHH0qSFi9erJycHK1cuVLDhw9Xamqqhg8frurq6gapE8CVI+AAaDAhISH65JNPNHjwYN188816+umn9fLLL2vQoEF1tn/yySfl6+ur2NhYtWrVynFD79///nfV1NRo4MCBiouL06RJkxQWFuYUQK7W/Pnzdfvtt+vuu+9WUlKS+vbtq/j4eKc2y5cvV1pamqZMmaLOnTsrNTVV27dvV9u2bbVv3z5NnTpVf/jDH2S1WiVJf/jDH3Tq1CnNnDmzweoEcGUsRl3PbwIAADRiXMEBAACmQ8ABAACmQ8ABAACmQ8ABAACmQ8ABAACmQ8ABAACmQ8ABAACmQ8ABAACmQ8ABAACmQ8ABAACmQ8ABAACm8/8AVbfNFraGRcEAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "f, ax = plt.subplots()\n",
    "\n",
    "ax.plot(np.arange(Nsites) - dmax, nocc, \"bo\", label=\"nocc\")\n",
    "ax.plot(np.arange(Nsites) - dmax, corr, \"ro\", label=\"corr\", alpha = 0.5)\n",
    "ax.legend()\n",
    "ax.set_xlabel(\"site index\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see, we have created a state with alternating spins. However, in the middle we have created a hole. This is the state that we will use as the initial state for the propagation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Time evolution\n",
    "\n",
    "To directly simulate the experiment, we run a list of circuits with different evolution times on the simulator backend. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the experimental parameters\n",
    "J = 2 * np.pi * 250  # hopping parameter in units of hbar\n",
    "U = 15 * J  # interaction parameter\n",
    "\n",
    "Ntimes = 8\n",
    "tmax = 1.5\n",
    "times = np.linspace(0, tmax, Ntimes) * 1e-3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "Nshots = 500\n",
    "# create list of circuits with interaction\n",
    "circuits_1 = [qc_1.bind_parameters({hop_param: J * t, int_param: U * t}) for t in times]\n",
    "# measure the observable from simulated shots\n",
    "jobs_1 = backend.run(circuits_1, shots=Nshots, num_species=2)\n",
    "\n",
    "# create list of circuits with interaction\n",
    "circuits_2 = [qc_2.bind_parameters({hop_param_down: J * t, int_param_down: U * t}) for t in times]\n",
    "# measure the observable from simulated shots\n",
    "jobs_2 = backend.run(circuits_2, shots=Nshots, num_species=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Put it together"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Duration for of the two jobs: 54.99120879173279 seconds\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/_v/vw5sh5fd4lv4x0cnk8cxz__w0000gn/T/ipykernel_82435/3602395813.py:110: RuntimeWarning: divide by zero encountered in divide\n",
      "  var_weight = 1 / Nobs_var\n",
      "/var/folders/_v/vw5sh5fd4lv4x0cnk8cxz__w0000gn/T/ipykernel_82435/3602395813.py:113: RuntimeWarning: divide by zero encountered in divide\n",
      "  mean_weight = 1 / Nobs\n"
     ]
    }
   ],
   "source": [
    "noccs = np.zeros((Ntimes, Nsites))\n",
    "corrs = np.zeros((Ntimes, Nsites))\n",
    "\n",
    "start_time = time.time()\n",
    "count_list_1 = jobs_1.result().get_counts()\n",
    "count_list_2 = jobs_2.result().get_counts()\n",
    "end_time = time.time()\n",
    "print(f\"Duration for of the two jobs: {end_time - start_time} seconds\")\n",
    "\n",
    "for ii in range(Ntimes):\n",
    "    counts_1 = count_list_1[ii]\n",
    "    counts_2 = count_list_2[ii]\n",
    "\n",
    "    total_counts = counts_1\n",
    "    total_counts.update(counts_2)\n",
    "\n",
    "    nocc = nocc_from_counts(total_counts, Nsites)\n",
    "    corr, _,_ = squeezed_spin_corr_from_counts(total_counts, Nsites)\n",
    "\n",
    "    noccs[ii, :] = nocc\n",
    "    corrs[ii, :] = corr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can directly visualize them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Density')"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAHHCAYAAABtF1i4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAwhUlEQVR4nO3de1iUdd7H8c8wyAyeUFNAXRSsNE1DQiU6rLqhZC1lPZarJoZluyVtRllSKZpP0UnTNV23NqWTaXbloXI1F0O3xAwMt1rNI0EmIKac2sCGef5onR5WUEaBn8O8X9d1X1dzz+833+9M8PPDPfc9Y3E6nU4BAAAY4mO6AQAA4N0IIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMwIi0tTRaLRbm5uaZbAdCMWCwWzZw503QbcBNhBPX2xRdfaNSoUerevbvsdru6du2qYcOGacGCBaZbA9DITv4BcXKz2+3q0qWLYmNj9ac//UllZWWmW6zV1q1bNXPmTB0/ftx0KzgNC99Ng/rYunWrhg4dqm7dumnChAkKDg5Wfn6+tm3bpv3792vfvn1uPZ7D4dCJEydks9lksVgaqWsADSUtLU0JCQl64oknFBYWphMnTqigoEAZGRnauHGjunXrprVr1+qyyy4z2uePP/4oX19f+fr6SpKef/55TZ06VQcPHlRoaKjR3lA3X9MNwDM8+eSTCggI0GeffaZ27drVuK+oqMjtx7NarbJarQ3UHYCmMmLECA0YMMB1Ozk5WZs2bdJvf/tb3Xjjjdq1a5f8/f2N9We3243VxtnjbRrUy/79+3XppZeeEkQkKTAw0PXfFotFiYmJevPNN9WrVy/Z7XZFRkZqy5YtNebUds5IaGiofvvb3+rjjz/WoEGDZLfb1aNHD7322muN9bQANIDf/OY3mj59ur755hu98cYbrv27d+/WqFGj1KFDB9ntdg0YMEBr166tMffkWvDJJ58oKSlJnTp1UqtWrXTzzTfryJEjNcZmZWUpNjZWHTt2lL+/v8LCwjRx4sQaY/7/OSMzZ87U1KlTJUlhYWGut5hyc3M1ePBghYeH1/p8evXqpdjY2HN9WeAGwgjqpXv37srOztaXX355xrGbN2/WlClTdPvtt+uJJ57Q0aNHdd1119Vr7r59+zRq1CgNGzZMc+bMUfv27XXHHXfoq6++aoinAaCRjB8/XpL04YcfSpK++uorXXHFFdq1a5emTZumOXPmqFWrVho5cqRWrVp1yvz77rtPO3fuVEpKiu655x699957SkxMdN1fVFSk4cOHKzc3V9OmTdOCBQs0btw4bdu2rc6ebrnlFo0ZM0aS9MILL+j111/X66+/rk6dOmn8+PH65z//ecq69Nlnn2nPnj26/fbbz/k1gRucQD18+OGHTqvV6rRarc7o6Gjnww8/7NywYYOzqqqqxjhJTknOrKws175vvvnGabfbnTfffLNr39KlS52SnAcPHnTt6969u1OSc8uWLa59RUVFTpvN5nzwwQcb78kBOKOTv7OfffZZnWMCAgKcERERTqfT6bz22mud/fr1c/7444+u+6urq51XXnml8+KLLz7lcWNiYpzV1dWu/Q888IDTarU6jx8/7nQ6nc5Vq1adsb7T+fMalJKS4rr93HPPnbLWOJ1O5/Hjx512u935yCOP1Nj/xz/+0dmqVStneXn5aeugYXFkBPUybNgwZWZm6sYbb9TOnTv17LPPKjY2Vl27dj3lsGt0dLQiIyNdt7t166abbrpJGzZskMPhOG2dPn366JprrnHd7tSpk3r16qUDBw407BMC0OBat26tsrIyff/999q0aZNuu+02lZWVqbi4WMXFxTp69KhiY2O1d+9eHTp0qMbcu+++u8bJ7Ndcc40cDoe++eYbSXK9Rfz+++/rxIkT59xrQECAbrrpJr311lty/uc6DofDoRUrVmjkyJFq1arVOddA/RFGUG8DBw7Uu+++q2PHjmn79u1KTk5WWVmZRo0apX/961+ucRdffPEpc3v27KkffvjhlPeA/1u3bt1O2de+fXsdO3bs3J8AgEZVXl6uNm3aaN++fXI6nZo+fbo6depUY0tJSZF06onv//273759e0ly/e4PHjxY//M//6NZs2apY8eOuummm7R06VJVVlaedb/x8fHKy8vTP/7xD0nS3//+dxUWFrreckLT4WoauM3Pz08DBw7UwIED1bNnTyUkJGjlypWuReZc1HWFjZMr0IHz2rfffquSkhJddNFFqq6uliQ99NBDdZ4IetFFF9W4fabffYvFonfeeUfbtm3Te++9pw0bNmjixImaM2eOtm3bptatW7vdc2xsrIKCgvTGG2/o17/+td544w0FBwcrJibG7cfCuSGM4JycvMTv8OHDrn179+49ZdyePXvUsmVLderUqcl6A9B0Xn/9dUk//wPfo0cPSVKLFi0a/B/2K664QldccYWefPJJLVu2TOPGjdPy5ct111131Tr+dJ9jZLVaNXbsWKWlpemZZ57R6tWrNWnSJD52wADepkG9fPTRR7UenVi3bp2kny+FOykzM1M7duxw3c7Pz9eaNWs0fPhwfsmBZmjTpk2aPXu2wsLCNG7cOAUGBmrIkCH6y1/+UuMPlZPO9HZtbY4dO3bKGtS/f39JOu1bNSfP/ajrE1jHjx+vY8eO6fe//73Ky8u5isYQjoygXu677z798MMPuvnmm3XJJZeoqqpKW7du1YoVKxQaGqqEhATX2L59+yo2NlZ//OMfZbPZtGjRIknSrFmzTLUPoIH87W9/0+7du/XTTz+psLBQmzZt0saNG9W9e3etXbvW9aFjCxcu1NVXX61+/fpp0qRJ6tGjhwoLC5WZmalvv/1WO3fudKvuq6++qkWLFunmm2/WhRdeqLKyMr388stq27atrr/++jrnnTyZ/rHHHtPvfvc7tWjRQnFxca6QEhERob59+2rlypXq3bu3Lr/88rN8ZXAuCCOol+eff14rV67UunXr9NJLL6mqqkrdunXTvffeq8cff7zGh6ENHjxY0dHRmjVrlvLy8tSnTx+lpaUZ/5hoAOduxowZkn4+d6xDhw7q16+f5s2bp4SEBLVp08Y1rk+fPsrKytKsWbOUlpamo0ePKjAwUBEREa7HcMfgwYO1fft2LV++XIWFhQoICNCgQYP05ptvKiwsrM55AwcO1OzZs7V48WKtX79e1dXVOnjwYI2rZeLj4/Xwww9z4qpBfDcNGpTFYtHkyZP14osvmm4FAOpl/vz5euCBB5Sbm1vrFX1ofJwzAgDwWk6nU6+88ooGDx5MEDGIt2kAAF6noqJCa9eu1UcffaQvvvhCa9asMd2SVyOMAAC8zpEjRzR27Fi1a9dOjz76qG688UbTLXk1zhkBAABGuX3OyJYtWxQXF6cuXbrIYrFo9erVZ5xTWVmpxx57TN27d5fNZlNoaKiWLFlyNv0CAIBmxu23aSoqKhQeHq6JEyfqlltuqdec2267TYWFhXrllVd00UUX6fDhw66PCwYAAN7N7TAyYsQIjRgxot7j169fr82bN+vAgQPq0KGDJCk0NNStmtXV1fruu+/Upk2b0360L4DG4XQ6VVZWpi5dusjHxzMuwmPdAMyr79rR6Cewrl27VgMGDNCzzz6r119/Xa1atdKNN96o2bNny9/fv16P8d133ykkJKSROwVwJvn5+frVr35luo16Yd0Azh9nWjsaPYwcOHBAH3/8sex2u1atWqXi4mLde++9Onr0qJYuXVrrnMrKyhrfNXDyHNu8rM1qexbfzNggWviZqfsf1Z9nGK3vWLPKaH3rTTcbre8T/muj9SVJlf82Vrq0vELdrri2xidsnm/qWjd+k/qefO2t6prWqDYvXWGkrku7QLP1Jcn0NRJOw6cElB41Wn7oxNFG6//0Y4X+Pi3ujGtHo4eR6upqWSwWvfnmmwoICJAkzZ07V6NGjdKiRYtqPTqSmppa6/eYtG3dWm3bmAojNjN1/6O6Vf2OIjUWh5/Zq8Cthp+/z/nwj7Cf+S8ZPJ/f7qhr3fC1t1ILfzPrhsVq9o8Y+drN1pcII4Z/Bkz97P+3M60djf7mb+fOndW1a1dXEJGk3r17y+l06ttvv611TnJyskpKSlxbfn5+Y7cJwMOxbgCeq9H/3L3qqqu0cuVKlZeXq/V/3mLZs2ePfHx86nz/yGazyWYzeyQCgGdh3QA8l9tHRsrLy5WTk6OcnBxJ0sGDB5WTk6O8vDxJP/91Eh8f7xo/duxYXXDBBUpISNC//vUvbdmyRVOnTtXEiRPrfQIrAABovtwOI1lZWYqIiFBERIQkKSkpqcZXQh8+fNgVTCSpdevW2rhxo44fP64BAwZo3LhxiouL05/+9KcGegoAAMCTuf02zZAhQ3S6T5BPS0s7Zd8ll1yijRs3ulsKAAB4Ac/49CIAANBsEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGOV2GNmyZYvi4uLUpUsXWSwWrV69ut5zP/nkE/n6+qp///7ulgUAAM2U22GkoqJC4eHhWrhwoVvzjh8/rvj4eF177bXulgQAAM2Yr7sTRowYoREjRrhd6A9/+IPGjh0rq9Xq1tEUAADQvDXJOSNLly7VgQMHlJKS0hTlAACAB3H7yIi79u7dq2nTpukf//iHfH3rV66yslKVlZWu26WlpY3VHoBmgnUD8FyNemTE4XBo7NixmjVrlnr27FnveampqQoICHBtISEhjdglgOaAdQPwXI0aRsrKypSVlaXExET5+vrK19dXTzzxhHbu3ClfX19t2rSp1nnJyckqKSlxbfn5+Y3ZJoBmgHUD8FyN+jZN27Zt9cUXX9TYt2jRIm3atEnvvPOOwsLCap1ns9lks9lOvaOFn9Silv1NwGIx/JEsfnaj5X3atDRaX63aGi1v8W1htL4kOS2tzBU/UW2udj3VtW4UFJXLanca6EhSyKVm6p6U/5XZ+pKsF11utL7jwE6j9dWl/u8KNIbDRyqM1nf8+EO9xrkdRsrLy7Vv3z7X7YMHDyonJ0cdOnRQt27dlJycrEOHDum1116Tj4+P+vbtW2N+YGCg7Hb7KfsBAIB3cjuMZGVlaejQoa7bSUlJkqQJEyYoLS1Nhw8fVl5eXsN1CAAAmjW3w8iQIUPkdNZ9yDMtLe2082fOnKmZM2e6WxYAADRTfDcNAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMMrtMLJlyxbFxcWpS5cuslgsWr169WnHv/vuuxo2bJg6deqktm3bKjo6Whs2bDjbfgEAQDPjdhipqKhQeHi4Fi5cWK/xW7Zs0bBhw7Ru3TplZ2dr6NChiouL0+eff+52swAAoPnxdXfCiBEjNGLEiHqPnzdvXo3bTz31lNasWaP33ntPERER7pYHAADNjNth5FxVV1errKxMHTp0qHNMZWWlKisrXbdLS0ubojUAHox1A/BcTX4C6/PPP6/y8nLddtttdY5JTU1VQECAawsJCWnCDgF4ItYNwHM1aRhZtmyZZs2apbfffluBgYF1jktOTlZJSYlry8/Pb8IuAXgi1g3AczXZ2zTLly/XXXfdpZUrVyomJua0Y202m2w2WxN1BqA5YN0APFeTHBl56623lJCQoLfeeks33HBDU5QEAAAewu0jI+Xl5dq3b5/r9sGDB5WTk6MOHTqoW7duSk5O1qFDh/Taa69J+vmtmQkTJmj+/PmKiopSQUGBJMnf318BAQEN9DQAAICncjuMZGVlaejQoa7bSUlJkqQJEyYoLS1Nhw8fVl5enuv+l156ST/99JMmT56syZMnu/afHO8Oi4+vLD5NfgHQz3ysZuqedKLKaHnHUbNXJvj+u8JofeP//yVZLOY+MNli9TNW+1wdPnRMPn6VZx7YCNp0bG+k7kllx+q+arGpOKrMrl1qYfatu3ZdOhutfyj/e6P1q6t+qNc4t/9lHzJkiJxOZ533/3fAyMjIcLcEAADwInw3DQAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAot8PIli1bFBcXpy5dushisWj16tVnnJORkaHLL79cNptNF110kdLS0s6iVQAA0By5HUYqKioUHh6uhQsX1mv8wYMHdcMNN2jo0KHKycnRlClTdNddd2nDhg1uNwsAAJofX3cnjBgxQiNGjKj3+MWLFyssLExz5syRJPXu3Vsff/yxXnjhBcXGxrpbHgAANDNuhxF3ZWZmKiYmpsa+2NhYTZkypc45lZWVqqysdN0uLS1trPYANBOsG4DnavQTWAsKChQUFFRjX1BQkEpLS/Xvf/+71jmpqakKCAhwbSEhIY3dJgAPx7oBeK7z8mqa5ORklZSUuLb8/HzTLQE4z7FuAJ6r0d+mCQ4OVmFhYY19hYWFatu2rfz9/WudY7PZZLPZGrs1AM0I6wbguRr9yEh0dLTS09Nr7Nu4caOio6MbuzQAAPAAboeR8vJy5eTkKCcnR9LPl+7m5OQoLy9P0s+HSuPj413j//CHP+jAgQN6+OGHtXv3bi1atEhvv/22HnjggYZ5BgAAwKO5HUaysrIUERGhiIgISVJSUpIiIiI0Y8YMSdLhw4ddwUSSwsLC9MEHH2jjxo0KDw/XnDlz9Ne//pXLegEAgKSzOGdkyJAhcjqddd5f26erDhkyRJ9//rm7pQAAgBc4L6+mAQAA3qPRr6ZpUC1skp+hs+VPczSoSRwtMlr+ULbZyyS7x31vtL6s58GvismfwfPh+Z+l0u8OyeJrN1K7fffuRuqe1DKsl9H6kvTDnp1G67fpfbnR+lZfq9H6Rw8eNFrf+dOP9RrHkREAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGHVWYWThwoUKDQ2V3W5XVFSUtm/fftrx8+bNU69eveTv76+QkBA98MAD+vHHH8+qYQAA0Ly4HUZWrFihpKQkpaSkaMeOHQoPD1dsbKyKiopqHb9s2TJNmzZNKSkp2rVrl1555RWtWLFCjz766Dk3DwAAPJ/bYWTu3LmaNGmSEhIS1KdPHy1evFgtW7bUkiVLah2/detWXXXVVRo7dqxCQ0M1fPhwjRkz5oxHUwAAgHdwK4xUVVUpOztbMTExvzyAj49iYmKUmZlZ65wrr7xS2dnZrvBx4MABrVu3Ttdff32ddSorK1VaWlpjA4DTYd0APJdbYaS4uFgOh0NBQUE19gcFBamgoKDWOWPHjtUTTzyhq6++Wi1atNCFF16oIUOGnPZtmtTUVAUEBLi2kJAQd9oE4IVYNwDP1ehX02RkZOipp57SokWLtGPHDr377rv64IMPNHv27DrnJCcnq6SkxLXl5+c3dpsAPBzrBuC5fN0Z3LFjR1mtVhUWFtbYX1hYqODg4FrnTJ8+XePHj9ddd90lSerXr58qKip0991367HHHpOPz6l5yGazyWazudMaAC/HugF4LreOjPj5+SkyMlLp6emufdXV1UpPT1d0dHStc3744YdTAofVapUkOZ1Od/sFAADNjFtHRiQpKSlJEyZM0IABAzRo0CDNmzdPFRUVSkhIkCTFx8era9euSk1NlSTFxcVp7ty5ioiIUFRUlPbt26fp06crLi7OFUoAAID3cjuMjB49WkeOHNGMGTNUUFCg/v37a/369a6TWvPy8mocCXn88cdlsVj0+OOP69ChQ+rUqZPi4uL05JNPNtyzAAAAHsvtMCJJiYmJSkxMrPW+jIyMmgV8fZWSkqKUlJSzKQUAAJo5vpsGAAAYRRgBAABGEUYAAIBRhBEAAGDUWZ3AaoyPj+TjpZcDH//eaPnPD5v9no/uhd8ZrS/LeZDbfSzmanvyZfglxZKvmQ9DO1bYxkjdk3pH9DBaX5J2HbQbrX9x785G63++ba/R+io2/EnEjsp6DTsPVlgAAODNCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADDqrMLIwoULFRoaKrvdrqioKG3fvv20448fP67Jkyerc+fOstls6tmzp9atW3dWDQMAgObF190JK1asUFJSkhYvXqyoqCjNmzdPsbGx+vrrrxUYGHjK+KqqKg0bNkyBgYF655131LVrV33zzTdq165dQ/QPAAA8nNthZO7cuZo0aZISEhIkSYsXL9YHH3ygJUuWaNq0aaeMX7Jkib7//ntt3bpVLVq0kCSFhoaeW9cAAKDZcOttmqqqKmVnZysmJuaXB/DxUUxMjDIzM2uds3btWkVHR2vy5MkKCgpS37599dRTT8nhcNRZp7KyUqWlpTU2ADgd1g3Ac7kVRoqLi+VwOBQUFFRjf1BQkAoKCmqdc+DAAb3zzjtyOBxat26dpk+frjlz5uh///d/66yTmpqqgIAA1xYSEuJOmwC8EOsG4Lka/Wqa6upqBQYG6qWXXlJkZKRGjx6txx57TIsXL65zTnJyskpKSlxbfn5+Y7cJwMOxbgCey61zRjp27Cir1arCwsIa+wsLCxUcHFzrnM6dO6tFixayWq2ufb1791ZBQYGqqqrk5+d3yhybzSabzeZOawC8HOsG4LncOjLi5+enyMhIpaenu/ZVV1crPT1d0dHRtc656qqrtG/fPlVXV7v27dmzR507d641iAAAAO/i9ts0SUlJevnll/Xqq69q165duueee1RRUeG6uiY+Pl7Jycmu8ffcc4++//573X///dqzZ48++OADPfXUU5o8eXLDPQsAAOCx3L60d/To0Tpy5IhmzJihgoIC9e/fX+vXr3ed1JqXlycfn18yTkhIiDZs2KAHHnhAl112mbp27ar7779fjzzySMM9CwAA4LHcDiOSlJiYqMTExFrvy8jIOGVfdHS0tm3bdjalAABAM8d30wAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMOqsPPTPH8p/NC/UON1o+p7zSaP2Rl15utL58vDy3Wzz4966yQvrphJnahQfM1P2PqqruRutL0oCbYozWLyn50Wh95+F9Ruvrpyqz9R31+93z8hUWAACYRhgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUWcVRhYuXKjQ0FDZ7XZFRUVp+/bt9Zq3fPlyWSwWjRw58mzKAgCAZsjtMLJixQolJSUpJSVFO3bsUHh4uGJjY1VUVHTaebm5uXrooYd0zTXXnHWzAACg+XE7jMydO1eTJk1SQkKC+vTpo8WLF6tly5ZasmRJnXMcDofGjRunWbNmqUePHufUMAAAaF7cCiNVVVXKzs5WTEzMLw/g46OYmBhlZmbWOe+JJ55QYGCg7rzzznrVqaysVGlpaY0NAE6HdQPwXG6FkeLiYjkcDgUFBdXYHxQUpIKCglrnfPzxx3rllVf08ssv17tOamqqAgICXFtISIg7bQLwQqwbgOdq1KtpysrKNH78eL388svq2LFjveclJyerpKTEteXn5zdilwCaA9YNwHP5ujO4Y8eOslqtKiwsrLG/sLBQwcHBp4zfv3+/cnNzFRcX59pXXV39c2FfX3399de68MILT5lns9lks9ncaQ2Al2PdADyXW0dG/Pz8FBkZqfT0dNe+6upqpaenKzo6+pTxl1xyib744gvl5OS4thtvvFFDhw5VTk4Oh1EBAIB7R0YkKSkpSRMmTNCAAQM0aNAgzZs3TxUVFUpISJAkxcfHq2vXrkpNTZXdblffvn1rzG/Xrp0knbIfAAB4J7fDyOjRo3XkyBHNmDFDBQUF6t+/v9avX+86qTUvL08+PnywKwAAqB+3w4gkJSYmKjExsdb7MjIyTjs3LS3tbEoCAIBmikMYAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjDqr76YxxmL5efNCPn2vNlr/yAmH0fqmn7/knT93v/Dg5++s/nkzoerfZur+x/6de43Wl6SPF9xutP7V971htL4qK8zW9xAcGQEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYdVZhZOHChQoNDZXdbldUVJS2b99e59iXX35Z11xzjdq3b6/27dsrJibmtOMBAIB3cTuMrFixQklJSUpJSdGOHTsUHh6u2NhYFRUV1To+IyNDY8aM0UcffaTMzEyFhIRo+PDhOnTo0Dk3DwAAPJ/bYWTu3LmaNGmSEhIS1KdPHy1evFgtW7bUkiVLah3/5ptv6t5771X//v11ySWX6K9//auqq6uVnp5+zs0DAADP51YYqaqqUnZ2tmJiYn55AB8fxcTEKDMzs16P8cMPP+jEiRPq0KGDe50CAIBmydedwcXFxXI4HAoKCqqxPygoSLt3767XYzzyyCPq0qVLjUDz3yorK1VZWem6XVpa6k6bALwQ6wbguZr0apqnn35ay5cv16pVq2S32+scl5qaqoCAANcWEhLShF0C8ESsG4DnciuMdOzYUVarVYWFhTX2FxYWKjg4+LRzn3/+eT399NP68MMPddlll512bHJyskpKSlxbfn6+O20C8EKsG4DnciuM+Pn5KTIyssbJpydPRo2Ojq5z3rPPPqvZs2dr/fr1GjBgwBnr2Gw2tW3btsYGAKfDugF4LrfOGZGkpKQkTZgwQQMGDNCgQYM0b948VVRUKCEhQZIUHx+vrl27KjU1VZL0zDPPaMaMGVq2bJlCQ0NVUFAgSWrdurVat27dgE8FAAB4IrfDyOjRo3XkyBHNmDFDBQUF6t+/v9avX+86qTUvL08+Pr8ccPnzn/+sqqoqjRo1qsbjpKSkaObMmefWPQAA8HhuhxFJSkxMVGJiYq33ZWRk1Lidm5t7NiUAAICX4LtpAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABglK/pBlA/llYBRusvmHSl0fqmn7/Xs1hMd+CZHCfM1v/2X2brS7r0V23NNnAevAY4M46MAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjziqMLFy4UKGhobLb7YqKitL27dtPO37lypW65JJLZLfb1a9fP61bt+6smgUAAM2P22FkxYoVSkpKUkpKinbs2KHw8HDFxsaqqKio1vFbt27VmDFjdOedd+rzzz/XyJEjNXLkSH355Zfn3DwAAPB8boeRuXPnatKkSUpISFCfPn20ePFitWzZUkuWLKl1/Pz583Xddddp6tSp6t27t2bPnq3LL79cL7744jk3DwAAPJ9bYaSqqkrZ2dmKiYn55QF8fBQTE6PMzMxa52RmZtYYL0mxsbF1jgcAAN7F153BxcXFcjgcCgoKqrE/KChIu3fvrnVOQUFBreMLCgrqrFNZWanKykrX7dLSUnfaBOCFWDcAz3VeXk2TmpqqgIAA1xYSEmK6JQDnOdYNwHO5FUY6duwoq9WqwsLCGvsLCwsVHBxc65zg4GC3xktScnKySkpKXFt+fr47bQLwQqwbgOdyK4z4+fkpMjJS6enprn3V1dVKT09XdHR0rXOio6NrjJekjRs31jlekmw2m9q2bVtjA4DTYd0APJdb54xIUlJSkiZMmKABAwZo0KBBmjdvnioqKpSQkCBJio+PV9euXZWamipJuv/++zV48GDNmTNHN9xwg5YvX66srCy99NJLDftMAACAR3I7jIwePVpHjhzRjBkzVFBQoP79+2v9+vWuk1Tz8vLk4/PLAZcrr7xSy5Yt0+OPP65HH31UF198sVavXq2+ffs23LMAAAAey+0wIkmJiYlKTEys9b6MjIxT9t1666269dZbz6YUAABo5s7Lq2kAAID3IIwAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADDK13QD9eF0OiVJpWVlhjvxXieqfjJav0VpqdH63u7k797J30VPcLJXp6PKcCferdTw7y7//806+fqfae2wOD1gdfn2228VEhJiug3A6+Xn5+tXv/qV6TbqhXUDOH+cae3wiDBSXV2t7777Tm3atJHFYjHdjttKS0sVEhKi/Px8tW3b1nQ7Rnj7a+Dpz9/pdKqsrExdunSRj49nvLvr6euG5Pk/N+fK25+/5PmvQX3XDo94m8bHx8dj/ho7nbZt23rkD1ND8vbXwJOff0BAgOkW3NJc1g3Js39uGoK3P3/Js1+D+qwdnvEnDgAAaLYIIwAAwCjCSBOw2WxKSUmRzWYz3Yox3v4aePvzx9nx9p8bb3/+kve8Bh5xAisAAGi+ODICAACMIowAAACjCCMAAMAowggAADCKMNKEcnNzdeeddyosLEz+/v668MILlZKSoqqq5v3dCQsXLlRoaKjsdruioqK0fft20y01mdTUVA0cOFBt2rRRYGCgRo4cqa+//tp0W/AwrB3etXZ447pBGGlCu3fvVnV1tf7yl7/oq6++0gsvvKDFixfr0UcfNd1ao1mxYoWSkpKUkpKiHTt2KDw8XLGxsSoqKjLdWpPYvHmzJk+erG3btmnjxo06ceKEhg8froqKCtOtwYOwdnjX2uGN6waX9hr23HPP6c9//rMOHDhgupVGERUVpYEDB+rFF1+U9PP3hYSEhOi+++7TtGnTDHfX9I4cOaLAwEBt3rxZv/71r023Aw/G2uE9vGHd4MiIYSUlJerQoYPpNhpFVVWVsrOzFRMT49rn4+OjmJgYZWZmGuzMnJKSEklqtv/P0XRYO7yHN6wbhBGD9u3bpwULFuj3v/+96VYaRXFxsRwOh4KCgmrsDwoKUkFBgaGuzKmurtaUKVN01VVXqW/fvqbbgQdj7fAe3rJuEEYawLRp02SxWE677d69u8acQ4cO6brrrtOtt96qSZMmGeocTWny5Mn68ssvtXz5ctOt4DzB2oEz8ZZ1w9d0A83Bgw8+qDvuuOO0Y3r06OH67++++05Dhw7VlVdeqZdeeqmRuzOnY8eOslqtKiwsrLG/sLBQwcHBhroyIzExUe+//762bNnSbL7WHueOtaN2rB0/86Z1gzDSADp16qROnTrVa+yhQ4c0dOhQRUZGaunSpfLxab4Hp/z8/BQZGan09HSNHDlS0s+HHNPT05WYmGi2uSbidDp13333adWqVcrIyFBYWJjplnAeYe2onbevHd64bhBGmtChQ4c0ZMgQde/eXc8//7yOHDniuq+5pv2kpCRNmDBBAwYM0KBBgzRv3jxVVFQoISHBdGtNYvLkyVq2bJnWrFmjNm3auN7vDggIkL+/v+Hu4ClYO7xr7fDKdcOJJrN06VKnpFq35mzBggXObt26Of38/JyDBg1ybtu2zXRLTaau/99Lly413Ro8CGuHd60d3rhu8DkjAADAqOb7piMAAPAIhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYge644w7XRy43pLS0NLVr1+6cH2fIkCGaMmXKOT8OgIbDuoGGxMfBQ/Pnz9f//+y7IUOGqH///po3b945Pe7o0aN1/fXXn2N3AM5HrBtoSIQRKCAgoFEe19/fv/l+jwLg5Vg30JB4m8ZLvPPOO+rXr5/8/f11wQUXKCYmRhUVFZJqHm694447tHnzZs2fP18Wi0UWi0W5ubmSpC+//FIjRoxQ69atFRQUpPHjx6u4uLjOmv99uHXmzJnq37+/Xn/9dYWGhiogIEC/+93vVFZW5hpTUVGh+Ph4tW7dWp07d9acOXNOedzKyko99NBD6tq1q1q1aqWoqChlZGRIkn788Uddeumluvvuu13j9+/frzZt2mjJkiVn+eoB3ol1g3WjqRBGvMDhw4c1ZswYTZw4Ubt27VJGRoZuueUW1fa1RPPnz1d0dLQmTZqkw4cP6/DhwwoJCdHx48f1m9/8RhEREcrKytL69etVWFio2267za1e9u/fr9WrV+v999/X+++/r82bN+vpp5923T916lRt3rxZa9as0YcffqiMjAzt2LGjxmMkJiYqMzNTy5cv1z//+U/deuutuu6667R3717Z7Xa9+eabevXVV7VmzRo5HA7dfvvtGjZsmCZOnHh2LyDghVg3WDealMlv6UPTyM7Odkpy5ubm1nr/hAkTnDfddJPr9uDBg533339/jTGzZ892Dh8+vMa+/Px8pyTn119/XevjLl261BkQEOC6nZKS4mzZsqWztLTUtW/q1KnOqKgop9PpdJaVlTn9/Pycb7/9tuv+o0ePOv39/V39fPPNN06r1eo8dOhQjVrXXnutMzk52XX72WefdXbs2NGZmJjo7Ny5s7O4uLjWHgHUjnWDdaMpcc6IFwgPD9e1116rfv36KTY2VsOHD9eoUaPUvn37ej/Gzp079dFHH6l169an3Ld//3717NmzXo8TGhqqNm3auG537txZRUVFrsepqqpSVFSU6/4OHTqoV69erttffPGFHA7HKfUqKyt1wQUXuG4/+OCDWr16tV588UX97W9/q3EfgDNj3WDdaEqEES9gtVq1ceNGbd26VR9++KEWLFigxx57TJ9++qnCwsLq9Rjl5eWKi4vTM888c8p9nTt3rncvLVq0qHHbYrGourq63vPLy8tltVqVnZ0tq9Va477/v+AVFRVpz549slqt2rt3r6677rp61wDAusG60bQ4Z8RLWCwWXXXVVZo1a5Y+//xz+fn5adWqVbWO9fPzk8PhqLHv8ssv11dffaXQ0FBddNFFNbZWrVo1SI8XXnihWrRooU8//dS179ixY9qzZ4/rdkREhBwOh4qKik7pIzg42DVu4sSJ6tevn1599VU98sgj2rVrV4P0CHgT1g3WjaZCGPECn376qZ566illZWUpLy9P7777ro4cOaLevXvXOj40NFSffvqpcnNzVVxcrOrqak2ePFnff/+9xowZo88++0z79+/Xhg0blJCQcMoCdLZat26tO++8U1OnTtWmTZv05Zdf6o477pCPzy8/pj179tS4ceMUHx+vd999VwcPHtT27duVmpqqDz74QJK0cOFCZWZm6tVXX9W4ceM0cuRIjRs3TlVVVQ3SJ+ANWDdYN5oSYcQLtG3bVlu2bNH111+vnj176vHHH9ecOXM0YsSIWsc/9NBDslqt6tOnjzp16qS8vDx16dJFn3zyiRwOh4YPH65+/fppypQpateuXY1f+nP13HPP6ZprrlFcXJxiYmJ09dVXKzIyssaYpUuXKj4+Xg8++KB69eqlkSNH6rPPPlO3bt20e/duTZ06VYsWLVJISIgkadGiRSouLtb06dMbrE+guWPdYN1oShans5brtAAAAJoIR0YAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABG/R+2UCv64tU2OAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "f, [ax1, ax2] = plt.subplots(1, 2, sharey=True)\n",
    "# color plot of the occupation in a heatmap\n",
    "velo = J\n",
    "ax1.pcolor(np.arange(Nsites) - dmax, times * 1e3, corrs, cmap=\"Reds\")\n",
    "\n",
    "ax1.set_xlabel(\"time [ms]\")\n",
    "ax1.set_xlabel(\"site index\")\n",
    "ax1.set_title(\"Spin\")\n",
    "\n",
    "ax2.pcolor(np.arange(Nsites) - dmax, times * 1e3, noccs, cmap=\"Blues\")\n",
    "ax2.set_xlabel(\"site index\")\n",
    "ax2.set_title(\"Density\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can nicely observe a much slower spread in the spin than in the charge density. This is the spin-charge separation. Now it is time to quantify the difference in the speed of the spread"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "stds_spin = np.zeros(Ntimes)\n",
    "for ii, corr in enumerate(corrs):\n",
    "    mass_func = corr + 1\n",
    "    mean_position = np.average(np.arange(len(mass_func)), weights=mass_func)\n",
    "    # Calculate the standard deviation of the position\n",
    "    std_position = np.sqrt(\n",
    "        np.average((np.arange(len(mass_func)) - mean_position) ** 2, weights=mass_func)\n",
    "    )\n",
    "    stds_spin[ii] = std_position\n",
    "\n",
    "stds_density = np.zeros(Ntimes)\n",
    "for ii, nocc in enumerate(noccs):\n",
    "    # print(nocc)\n",
    "    mass_func = 1 - nocc\n",
    "    mean_position = np.average(np.arange(len(mass_func)), weights=mass_func)\n",
    "    # Calculate the standard deviation of the position\n",
    "    std_position = np.sqrt(\n",
    "        np.average((np.arange(len(mass_func)) - mean_position) ** 2, weights=mass_func)\n",
    "    )\n",
    "    stds_density[ii] = std_position"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x16a0d2140>"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAGwCAYAAABVdURTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA/p0lEQVR4nO3dfVxUZf7/8fcwCogC3hSCgkKa5g3el6KxYZlWriuxbq72yzStLbU0sza3NrV2s++WrdZabrre1KZWLtZummkm6npTKdpqa1aKgghiN4JCeTOc3x+zTI6iceDMDBxez8djHuNcc805nyPFvL3Oda7jMAzDEAAAgE0EBboAAAAAKxFuAACArRBuAACArRBuAACArRBuAACArRBuAACArRBuAACArdQJdAH+VlpaqiNHjig8PFwOhyPQ5QAAgAowDEMnTpxQs2bNFBR06bGZWhdujhw5ori4uECXAQAAKiEnJ0exsbGX7FPrwk14eLgk919OREREgKsBAAAVUVRUpLi4OM/3+KXUunBTdioqIiKCcAMAQA1TkSklTCgGAAC2QrgBAAC2QrgBAAC2Uuvm3FSUy+XSmTNnAl1GrVa3bl05nc5AlwEAqGEIN+cxDEP5+fk6fvx4oEuBpIYNGyo6Opo1iQAAFUa4OU9ZsImKilJYWBhfqgFiGIZKSkpUUFAgSYqJiQlwRQCAmoJwcw6Xy+UJNk2aNAl0ObVevXr1JEkFBQWKioriFBUAoEKYUHyOsjk2YWFhAa4EZcp+Fsx/AgBUFOGmHJyKqj74WQAAzOK0FAAAsITLJW3aJOXlSTExUnKyFIgZBYQbAABQZenp0oQJ0uHDP7bFxkqzZ0tpaf6thdNSvuJySRkZ0tKl7meXK9AVXVRKSoomTpwY6DIAADVUero0ZIh3sJGk3Fx3e3q6f+th5MYXqlN8rYD09HTVrVs30GUAAGogl8v9lWcYF75nGJLDIU2cKA0e7L9TVIzcWK26xdcKaNy4cYVuIQ8AwPk2bbrwK+9chiHl5Lj7+Qvhxko/FV8ld3z1wSmq5cuXKzExUfXq1VOTJk3Ur18/FRcXa+TIkUpNTdX06dN1+eWXKyIiQvfee69Onz7t+ez5p6Xi4+P19NNP66677lJ4eLhatGihV155xfKaAQA1X16etf2sQLixUoDia15enoYNG6a77rpLe/fuVUZGhtLS0mT8L1CtW7fO07506VKlp6dr+vTpl9zmzJkz1aNHD+3cuVNjx47Vfffdp3379llaNwCg5qvoAvL+XGiecGOlAMXXvLw8nT17VmlpaYqPj1diYqLGjh2rBg0aSJKCg4O1YMECdejQQQMHDtSTTz6pF154QaWlpRfd5i233KKxY8eqdevW+u1vf6vLLrtM69evt7RuAEDNl5zsnlZ6sWXJHA4pLs7dz18IN1YKUHzt3LmzbrjhBiUmJupXv/qV5s2bp++++87r/XNXXU5KStLJkyeVk5Nz0W126tTJ82eHw6Ho6GjPfZ4AACjjdLqvl5EuDDhlr2fN8u96N4QbKwUovjqdTq1du1bvvfee2rdvrxdffFFt27ZVVlZWpbd5/tVTDofjkiM9AIDaKy1NWr5cat7cuz021t3OOjc1WQDjq8PhUJ8+fTR9+nTt3LlTwcHBWrFihSTp008/1ffff+/pu23bNjVo0EBxcXGW1wEAqJ3S0qSDB6X166UlS9zPWVmBWQGFcGO1AMTXjz76SE8//bS2b9+u7Oxspaen69ixY2rXrp0k6fTp0xo9erT++9//atWqVZo6darGjx+voCB+/AAA6zidUkqKNGyY+zkQt16QWMTPN9LS3KsV+ekGGxEREdq4caNmzZqloqIitWzZUjNnztTNN9+sN954QzfccIOuvPJK/exnP9OpU6c0bNgwTZs2zSe1AAAQaA7DKG9RFvsqKipSZGSkCgsLFRER4fXeDz/8oKysLCUkJCg0NDRAFVpr5MiROn78uN5+++1Al1IpdvyZAADMu9T39/k4LwEAAGyFcAMAAGyFOTc2t2jRokCXAACAXzFyAwAAbIVwAwAAbIVwAwAAbIU5NwAAWMDl8tvyZvgJhBsAAKooPV2aMEE6fPjHtthY9x15AnH7gdqO01I2lpKSookTJ/plX9OmTVOXLl38si8AqE7S06UhQ7yDjSTl5rrb09MDU1dtRrjxEZdLysiQli51P7tcga7ItyZPnqx169Z5Xo8cOVKpqamBKwgA/MDlco/YlLfWf1nbxIn2/w6obgg3PpCeLsXHS337SsOHu5/j4+2d3hs0aKAmTZoEugwA8KtNmy4csTmXYUg5Oe5+8B/CjcUCNTxZXFysESNGqEGDBoqJidHMmTO93j916pQmT56s5s2bq379+urZs6cyMjI87y9atEgNGzbU+++/r3bt2qlBgwa66aablJeX5+mTkZGha665RvXr11fDhg3Vp08fHTp0SJL3aalp06Zp8eLFeuedd+RwOORwOJSRkaHrr79e48eP96rr2LFjCg4O9hr1AYCa4pxfkZb0gzUINxYK5PDkww8/rA0bNuidd97RmjVrlJGRoczMTM/748eP19atW7Vs2TL95z//0a9+9SvddNNN+vLLLz19SkpK9Nxzz+m1117Txo0blZ2drcmTJ0uSzp49q9TUVF133XX6z3/+o61bt+qee+6Rw+G4oJbJkyfrtttu84SjvLw89e7dW2PGjNGSJUt06tQpT9+///3vat68ua6//nrr/1IAwMdiYqztB2sQbiwUqOHJkydP6m9/+5uee+453XDDDUpMTNTixYt19uxZSVJ2drYWLlyot956S8nJyWrVqpUmT56sa6+9VgsXLvRs58yZM5o7d6569Oihbt26afz48Z4RlaKiIhUWFurnP/+5WrVqpXbt2unOO+9UixYtLqinQYMGqlevnkJCQhQdHa3o6GgFBwcr7X+XDLzzzjuevosWLdLIkSPLDUkAUN0lJ7uvirrYrzCHQ4qLc/eD/xBuLBSo4cn9+/fr9OnT6tmzp6etcePGatu2rSRp9+7dcrlcatOmjRo0aOB5bNiwQfv37/d8JiwsTK1atfK8jomJUUFBgWd7I0eO1IABAzRo0CDNnj3b65RVRYSGhuqOO+7QggULJEmZmZnas2ePRo4cWdlDB4CAcjrdl3tLFwacstezZrHejb8FNNzMmDFDV199tcLDwxUVFaXU1FTt27fvkp9ZtGiRZx5H2SM0NNRPFV9adR2ePHnypJxOp3bs2KFdu3Z5Hnv37tXssv8rJdWtW9frcw6HQ8Y559gWLlyorVu3qnfv3nrjjTfUpk0bbdu2zVQtY8aM0dq1a3X48GEtXLhQ119/vVq2bFm1AwSAAEpLk5Yvl5o3926PjXW3s86N/wU03GzYsEHjxo3Ttm3btHbtWp05c0b9+/dXcXHxJT8XERHhmcuRl5fnmdQaaIEanmzVqpXq1q2rjz76yNP23Xff6YsvvpAkde3aVS6XSwUFBWrdurXXIzo62tS+unbtqilTpmjLli3q2LGjlixZUm6/4OBgucqZXJSYmKgePXpo3rx5WrJkie666y5T+weA6igtTTp4UFq/XlqyxP2clUWwCZSArlC8evVqr9eLFi1SVFSUduzYoZ/97GcX/ZzD4ajwl/KpU6e8JrAWFRVVrtgKKBueHDLEHWTOnVjsy+HJBg0aaPTo0Xr44YfVpEkTRUVF6bHHHlNQkDu7tmnTRrfffrtGjBihmTNnqmvXrjp27JjWrVunTp06aeDAgT+5j6ysLL3yyiv6xS9+oWbNmmnfvn368ssvNWLEiHL7x8fH6/3339e+ffvUpEkTRUZGekaGxowZo/Hjx6t+/fq69dZbrfuLAIAAcjqllJRAVwGpms25KSwslOSe33EpJ0+eVMuWLRUXF6fBgwfrs88+u2jfGTNmKDIy0vOIi4uztObzBWp48tlnn1VycrIGDRqkfv366dprr1X37t097y9cuFAjRozQQw89pLZt2yo1NVWffPJJuROCyxMWFqbPP/9cv/zlL9WmTRvdc889GjdunH7zm9+U2//uu+9W27Zt1aNHD11++eXavHmz571hw4apTp06GjZsWLU5pQgAsA+HYZR34bL/lZaW6he/+IWOHz+uf//73xftt3XrVn355Zfq1KmTCgsL9dxzz2njxo367LPPFBsbe0H/8kZu4uLiVFhYqIiICK++P/zwg7KyspSQkFDlL11uoHZxBw8eVKtWrfTJJ5+oW7dul+xr5c8EAFBzFRUVKTIystzv7/NVmxtnjhs3Tnv27LlksJGkpKQkJSUleV737t1b7dq101//+lc99dRTF/QPCQlRSEiI5fX+FIYnL3TmzBl98803evzxx9WrV6+fDDYAAFRGtTgtNX78eL377rtav359uaMvl1K3bl117dpVX331lY+qg1U2b96smJgYffLJJ5o7d26gywEA2FRAR24Mw9D999+vFStWKCMjQwkJCaa34XK5tHv3bt1yyy0+qBBWSklJUTU5CwoAsLGAhptx48ZpyZIleueddxQeHq78/HxJUmRkpOrVqydJGjFihJo3b64ZM2ZIkp588kn16tVLrVu31vHjx/Xss8/q0KFDGjNmTMCOAwAAVB8BDTcvv/yyJPe/6M+1cOFCz6q12dnZnkuaJff6LXfffbfy8/PVqFEjde/eXVu2bFH79u0tq4vRheqDnwUAwKxqc7WUv1xqtrXL5dIXX3yhqKgoNWnSJEAV4lzffPONCgoK1KZNGzm53AwAaq0aebVUdeB0OtWwYUPP/ZTCwsK4oWOAGIahkpISFRQUqGHDhgQbAECFEW7OU7bycVnAQWA1bNjQ9C0iAAC1G+HmPA6HQzExMYqKitKZM2cCXU6tVrduXUZsAACmEW4uwul08sUKAEANVC0W8QMAALAK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK69wAACzhckmbNkl5eVJMjJScLLFcGAKBcAMAqLL0dGnCBOnw4R/bYmOl2bOltLTA1YXaidNSAIAqSU+XhgzxDjaSlJvrbk9PD0xdqL0INwCASnO53CM2hnHhe2VtEye6+wH+QrgBAFTapk0XjticyzCknBx3P8BfCDcAgErLy7O2H2AFwg0AoNJiYqztB1iBcAMAqLTkZPdVUQ5H+e87HFJcnLsf4C+EGwBApTmd7su9pQsDTtnrWbNY7wb+RbgBAFRJWpq0fLnUvLl3e2ysu511buBvLOIHAKiytDRp8GBWKEb1QLgBAFjC6ZRSUgJdBcBpKQAAYDOEGwAAYCuEGwAAYCuEGwAAYCuEGwAAYCuEGwAAYCuEGwAAYCuEGwAAYCuEGwAAYCuEGwAAYCuEGwAAYCuEGwAAYCuEGwAAYCuEGwAAYCuEGwAAYCuEGwAAYCuEGwAAYCuEGwAAYCuEGwAAYCuEGwAAYCuEGwAAYCuEGwAAYCuEGwAAYCuEGwAAYCuEGwAAYCuEGwAAYCuEGwAAYCsBDTczZszQ1VdfrfDwcEVFRSk1NVX79u37yc+99dZbuuqqqxQaGqrExEStWrXKD9UCAICaIKDhZsOGDRo3bpy2bdumtWvX6syZM+rfv7+Ki4sv+pktW7Zo2LBhGj16tHbu3KnU1FSlpqZqz549fqwcAABUVw7DMIxAF1Hm2LFjioqK0oYNG/Szn/2s3D5Dhw5VcXGx3n33XU9br1691KVLF82dO/cn91FUVKTIyEgVFhYqIiLCstoBAIDvmPn+rlZzbgoLCyVJjRs3vmifrVu3ql+/fl5tAwYM0NatW8vtf+rUKRUVFXk9AACAfVWbcFNaWqqJEyeqT58+6tix40X75efnq2nTpl5tTZs2VX5+frn9Z8yYocjISM8jLi7O0roBAED1Um3Czbhx47Rnzx4tW7bM0u1OmTJFhYWFnkdOTo6l2wcAANVLnUAXIEnjx4/Xu+++q40bNyo2NvaSfaOjo3X06FGvtqNHjyo6Orrc/iEhIQoJCbGsVgAAUL0FdOTGMAyNHz9eK1as0IcffqiEhISf/ExSUpLWrVvn1bZ27VolJSX5qkwAAFCDBHTkZty4cVqyZIneeecdhYeHe+bNREZGql69epKkESNGqHnz5poxY4YkacKECbruuus0c+ZMDRw4UMuWLdP27dv1yiuvBOw4AABA9RHQkZuXX35ZhYWFSklJUUxMjOfxxhtvePpkZ2crLy/P87p3795asmSJXnnlFXXu3FnLly/X22+/fclJyADgay6XlJEhLV3qfna5Al0RUHtVq3Vu/IF1bgBYLT1dmjBBOnz4x7bYWGn2bCktLXB1AXZSY9e5AYCaJj1dGjLEO9hIUm6uuz09PTB1AbUZ4QYAKsnlco/YlDf+XdY2cSKnqAB/I9wAQCVt2nThiM25DEPKyXH3A+A/hBsAqKRzrnWwpB8AaxBuAKCSYmKs7QfAGoQbAKik5GT3VVEOR/nvOxxSXJy7HwD/IdwAQCU5ne7LvaULA07Z61mz3P0A+A/hBgCqIC1NWr5cat7cuz021t3OOjeA/1WLG2cCQE2WliYNHuy+Kiovzz3HJjmZERsgUEyHm+LiYj3zzDNat26dCgoKVFpa6vX+gQMHLCsOAGoKp1NKSQl0FQCkSoSbMWPGaMOGDbrjjjsUExMjx8Vm0gEAAASA6XDz3nvvaeXKlerTp48v6gEAAKgS0xOKGzVqpMaNG/uiFgAAgCozHW6eeuopPfHEEyopKfFFPQAAAFVi+rTUzJkztX//fjVt2lTx8fGqW7eu1/uZmZmWFQcAAGCW6XCTmprqgzIAAACs4TAMwwh0Ef5UVFSkyMhIFRYWKiIiItDlAACACjDz/V3pRfx27NihvXv3SpI6dOigrl27VnZTAAAAljEdbgoKCvTrX/9aGRkZatiwoSTp+PHj6tu3r5YtW6bLL7/c6hoBAAAqzPTVUvfff79OnDihzz77TN9++62+/fZb7dmzR0VFRXrggQd8USMAAECFmZ5zExkZqQ8++EBXX321V/vHH3+s/v376/jx41bWZznm3AAAUPOY+f42PXJTWlp6weXfklS3bt0L7jMFAADgb6bDzfXXX68JEyboyJEjnrbc3Fw9+OCDuuGGGywtDgAAwCzT4eYvf/mLioqKFB8fr1atWqlVq1ZKSEhQUVGRXnzxRV/UCAAAUGGmr5aKi4tTZmamPvjgA33++eeSpHbt2qlfv36WFwcAAGAWi/gBAIBqz/JF/F544QXdc889Cg0N1QsvvHDJvlwODgAAAqlCIzcJCQnavn27mjRpooSEhItvzOHQgQMHLC3QaozcAABQ81g+cpOVlVXunwEAAKob01dLPfnkkyopKbmg/fvvv9eTTz5pSVEAAACVZXpCsdPpVF5enqKiorzav/nmG0VFRcnlcllaoNU4LQUAQM3j0xWKDcOQw+G4oP3TTz9V48aNzW4OAADAUhVe56ZRo0ZyOBxyOBxq06aNV8BxuVw6efKk7r33Xp8UCQAAUFEVDjezZs2SYRi66667NH36dEVGRnreCw4OVnx8vJKSknxSJAAAQEVVONzceeedktyXhffu3bvcm2cCAAAEmunbL1x33XWeP//www86ffq01/tM0gUAAIFkekJxSUmJxo8fr6ioKNWvX1+NGjXyegAAAASS6XDz8MMP68MPP9TLL7+skJAQzZ8/X9OnT1ezZs306quv+qJGAACACjN9Wupf//qXXn31VaWkpGjUqFFKTk5W69at1bJlS73++uu6/fbbfVEnAABAhZgeufn22291xRVXSHLPr/n2228lSddee602btxobXUAAAAmmQ43V1xxhef+UldddZXefPNNSe4RnYYNG1paHAAAgFmmw82oUaP06aefSpIeffRRzZkzR6GhoXrwwQf18MMPW14gAACAGabvLXW+Q4cOaceOHWrdurU6depkVV0+w72lAACoecx8f5ueUHy+li1bqmXLllXdDAAAgCUqFG5eeOEF3XPPPQoNDdULL7xwyb4PPPCAJYUBAABURoVOSyUkJGj79u1q0qSJEhISLr4xh0MHDhywtECrcVoKAICax/LTUmVXR53/ZwAAgOrG9NVS//73v31RBwAAgCVMh5vrr79eCQkJ+t3vfqfPPvvMFzUBAABUmulwc+TIET300EPasGGDEhMT1aVLFz377LM6fPiwL+oDAAAwpUrr3GRlZWnJkiVaunSpPv/8c/3sZz/Thx9+aGV9lmNCMQAANY+Z72/TIzfnSkhI0KOPPqpnnnlGiYmJ2rBhg6nPb9y4UYMGDVKzZs3kcDj09ttvX7J/RkaGHA7HBY/8/PwqHAUAALCTSoebzZs3a+zYsYqJidHw4cPVsWNHrVy50tQ2iouL1blzZ82ZM8fU5/bt26e8vDzPIyoqytTnAQCAfZleoXjKlClatmyZjhw5ohtvvFGzZ8/W4MGDFRYWZnrnN998s26++WbTn4uKiuImnQAAoFymw83GjRv18MMP67bbbtNll13mi5p+UpcuXXTq1Cl17NhR06ZNU58+fS7a99SpUzp16pTndVFRkT9KBACg9nG5pE2bpLw8KSZGSk6WnE6/l2E63GzevNkXdVRITEyM5s6dqx49eujUqVOaP3++UlJS9NFHH6lbt27lfmbGjBmaPn26nysFAKCWSU+XJkyQzr16OjZWmj1bSkvzaymVulrqtdde09y5c5WVlaWtW7eqZcuWmjVrlhISEjR48ODKFeJwaMWKFUpNTTX1ueuuu04tWrTQa6+9Vu775Y3cxMXFcbUUAABWSU+XhgyRzo8UDof7efnyKgccn14t9fLLL2vSpEm65ZZbdPz4cblcLklSw4YNNWvWrEoVXBXXXHONvvrqq4u+HxISooiICK8HAOu5XFJGhrR0qfv5f78aANidy+UesSlvrKSsbeJEv/5SMB1uXnzxRc2bN0+PPfaYnOecR+vRo4d2795taXEVsWvXLsXExPh9vwB+lJ4uxcdLfftKw4e7n+Pj3e0AbG7TJu9TUeczDCknx93PT0zPucnKylLXrl0vaA8JCVFxcbGpbZ08edJr1CUrK0u7du1S48aN1aJFC02ZMkW5ubl69dVXJclz6qtDhw764YcfNH/+fH344Ydas2aN2cMAYJGLjUbn5rrbLRiNBlCd5eVZ288CpsNNQkKCdu3apZYtW3q1r169Wu3atTO1re3bt6tv376e15MmTZIk3XnnnVq0aJHy8vKUnZ3tef/06dN66KGHlJubq7CwMHXq1EkffPCB1zYA+M9PjUY7HO7R6MGDA3LBBAB/qOjZEz+eZTE9oXj+/PmaNm2aZs6cqdGjR2v+/Pnav3+/ZsyYofnz5+vXv/61r2q1BLdfAKyTkeE+BfVT1q+XUlJ8XQ2AgHC53Oehc3PL/5eOw+G+aiorq0r/yjHz/W165GbMmDGqV6+eHn/8cZWUlGj48OFq1qyZZs+eXe2DDQBrVcPRaAD+5nS6L/ceMsQdZM4NOGVXS82a5dfh20rdfuH222/Xl19+qZMnTyo/P1+HDx/W6NGjra4NQDVXDUejAQRCWpp7gl3z5t7tsbEBmXhXpbuC10SclgKs46fRaAA1hQ9XKLb8tFTXrl3lKBta+gmZmZkV6geg5quGo9EAAsnprBYT7CoUbs5dNfiHH37QSy+9pPbt2yspKUmStG3bNn322WcaO3asT4oEUH2VjUaXt+r6rFlcBg7A/0yflhozZoxiYmL01FNPebVPnTpVOTk5WrBggaUFWo3TUoBvVJP75QGwKTPf36bDTWRkpLZv364rr7zSq/3LL79Ujx49VFhYaL5iPyLcAAB8goTvUz69t1S9evXKvTP45s2bFRoaanZzAADUfNyDpFoxvc7NxIkTdd999ykzM1PXXHONJOmjjz7SggUL9Pvf/97yAgEAqNa4B0m1U6lLwd98803Nnj1be/fulSS1a9dOEyZM0G233WZ5gVbjtBQAwDJl6yFc7MaRrIdgGZ/OuanpCDcAAMtwDxK/8emcGwAA8D/cg6RaItwAAFBZ3IOkWiLcAABQWcnJ7jk1F1vF3+GQ4uLc/eA3hBsAACqr7B4k0oUBh3uQBAzhBgCAqqhmd8RGBde5mTRpUoU3+Pzzz1e6GAAAaqS0NGnwYFYoriYqFG527tzp9TozM1Nnz55V27ZtJUlffPGFnE6nunfvbn2FAADUBNXkjtioYLhZv36958/PP/+8wsPDtXjxYjVq1EiS9N1332nUqFFKZsIUAAAIMNOL+DVv3lxr1qxRhw4dvNr37Nmj/v3768iRI5YWaDUW8QMAoObx6SJ+RUVFOnbs2AXtx44d04kTJ8xuDgAAwFKmw82tt96qUaNGKT09XYcPH9bhw4f1j3/8Q6NHj1YaM8IBAECAmb4r+Ny5czV58mQNHz5cZ86ccW+kTh2NHj1azz77rOUFAgAAmGFqzo3L5dLmzZuVmJio4OBg7d+/X5LUqlUr1a9f32dFWok5NwAA1Dxmvr9Njdw4nU71799fe/fuVUJCgjp16lSlQgEAAKxmes5Nx44ddeDAAV/UAgAAUGWmw80f/vAHTZ48We+++67y8vJUVFTk9QAA1FIul5SRIS1d6n52uQJdEWop0+vcBAX9mIcc59wkzDAMORwOuar5f8zMuQEAH0hPlyZMkA4f/rEtNtZ9U0mupIUFfDbnRvJerRgAAKWnS0OGSOf/Wzk3193OzSPhZ6ZHbmo6Rm4AwEIulxQf7z1icy6Hwz2Ck5XFTSRRJT4duSlTUlKi7OxsnT592qudK6gAoBbZtOniwUZyj+bk5Lj7cVNJ+InpcHPs2DGNGjVK7733XrnvV/c5NwAAC+XlWdsPsIDpq6UmTpyo48eP66OPPlK9evW0evVqLV68WFdeeaX++c9/+qJGAEB1FRNjbT/AAqZHbj788EO988476tGjh4KCgtSyZUvdeOONioiI0IwZMzRw4EBf1AkAqI6Sk91zanJzL5xQLP045yY52f+1odYyPXJTXFysqKgoSVKjRo08dwhPTExUZmamtdUBAKo3p9N9ubfkDjLnKns9axaTieFXpsNN27ZttW/fPklS586d9de//lW5ubmaO3euYhh2BIDaJy3Nfbl38+be7bGxXAaOgDB9WmrChAnK+9/EsKlTp+qmm27S66+/ruDgYC1atMjq+gAANUFamjR4sPuqqLw89xyb5GRGbBAQVV7npqSkRJ9//rlatGihyy67zKq6fIZ1bgAAqHn8ss5NmbCwMHXr1q2qmwEAALBEhcLNpEmTKrzB559/vtLFAAAAVFWFws3OnTu9XmdmZurs2bNq27atJOmLL76Q0+lU9+7dra8QAADAhAqFm3Nvlvn8888rPDxcixcvVqNGjSRJ3333nUaNGqVk1jEAAAABZnpCcfPmzbVmzRp16NDBq33Pnj3q37+/jhw5YmmBVmNCMQAANY+Z72/T69wUFRV5Fu4717Fjx3TixAmzmwMAALCU6XBz6623atSoUUpPT9fhw4d1+PBh/eMf/9Do0aOVxkJNAAAgwExfCj537lxNnjxZw4cP15kzZ9wbqVNHo0eP1rPPPmt5gQAAAGZUehG/4uJi7d+/X5LUqlUr1a9f39LCfIU5NwAA1Dx+WcSvfv366tSpU2U/DgAA4BOmw01xcbGeeeYZrVu3TgUFBSotLfV6/8CBA5YVBwAAYJbpcDNmzBht2LBBd9xxh2JiYuQ4/xb3AAAAAWQ63Lz33ntauXKl+vTp44t6AAAAqsT0peCNGjVS48aNLdn5xo0bNWjQIDVr1kwOh0Nvv/32T34mIyND3bp1U0hIiFq3bq1FixZZUgsAALAH0+Hmqaee0hNPPKGSkpIq77y4uFidO3fWnDlzKtQ/KytLAwcOVN++fbVr1y5NnDhRY8aM0fvvv1/lWgAAgD2YvhS8a9eu2r9/vwzDUHx8vOrWrev1fmZmZuUKcTi0YsUKpaamXrTPb3/7W61cuVJ79uzxtP3617/W8ePHtXr16grth0vBAQCoeXx6Kfilwoevbd26Vf369fNqGzBggCZOnHjRz5w6dUqnTp3yvC4qKvJVeQAAoBowHW6mTp3qizoqJD8/X02bNvVqa9q0qYqKivT999+rXr16F3xmxowZmj59ur9KBFBbuVzSpk1SXp4UEyMlJ0tOZ6CrAmol03NuapopU6aosLDQ88jJyQl0SQDsJj1dio+X+vaVhg93P8fHu9sB+J3pkRuXy6U///nPevPNN5Wdna3Tp097vf/tt99aVtz5oqOjdfToUa+2o0ePKiIiotxRG0kKCQlRSEiIz2oCUMulp0tDhkjnT1/MzXW3L18ucVNhwK9Mj9xMnz5dzz//vIYOHarCwkJNmjRJaWlpCgoK0rRp03xQ4o+SkpK0bt06r7a1a9cqKSnJp/sFKsLlkjIypKVL3c8uV6Args+5XNKECRcGG+nHtokT+Y8B8DPT4eb111/XvHnz9NBDD6lOnToaNmyY5s+fryeeeELbtm0zta2TJ09q165d2rVrlyT3pd67du1Sdna2JPcppREjRnj633vvvTpw4IAeeeQRff7553rppZf05ptv6sEHHzR7GIClOCtRS23aJB0+fPH3DUPKyXH3A+A3psNNfn6+EhMTJUkNGjRQYWGhJOnnP/+5Vq5caWpb27dvV9euXdW1a1dJ0qRJk9S1a1c98cQTkqS8vDxP0JGkhIQErVy5UmvXrlXnzp01c+ZMzZ8/XwMGDDB7GIBlys5KnP8dV3ZWgoBjY3l51vYDYAnTc25iY2OVl5enFi1aqFWrVlqzZo26deumTz75xPTclpSUFF1qmZ3yVh9OSUnRzp07zZYN+MRPnZVwONxnJQYP5sIZW4qJsbYfAEuYHrm59dZbPfNe7r//fv3+97/XlVdeqREjRuiuu+6yvECgOuOsRC2XnCzFxrpTbHkcDikuzt0PgN+YHrl55plnPH8eOnSoWrZsqS1btujKK6/UoEGDLC0OqO44K1HLOZ3S7Nnu848Oh/cQXlngmTWLYTvAz0yP3GzcuFFnz571vO7Vq5cmTZqkm2++WRs3brS0OKC646wElJbmvty7eXPv9thYLgMHAsT0vaWcTqfy8vIUFRXl1f7NN98oKipKrmp+ySP3loKVXC73VVG5ueXPu3E43N9xWVn84932WKEY8Cmf3lvKMAw5yjm//M0336h+/fpmNwfUaJyVgIfTKaWkBLoKADIRbtL+N7TqcDg0cuRIryujXC6X/vOf/6h3797WVwhUc2VnJSZM8J5cHBvrDjaclQAA/6pwuImMjJTkHrkJDw/3ut1BcHCwevXqpbvvvtv6CoEaIC3Nfbk3ZyUAIPAqHG4WLlwoSYqPj9fkyZM5BQWch7MSAFA9mL5a6pFHHvGac3Po0CHNmjVLa9assbQwAACAyjAdbgYPHqxXX31VknT8+HFdc801mjlzpgYPHqyXX37Z8gIBAADMMB1uMjMzlfy/1TaXL1+u6OhoHTp0SK+++qpeeOEFywsEAAAww3S4KSkpUXh4uCRpzZo1SktLU1BQkHr16qVDhw5ZXiAAAIAZpsNN69at9fbbbysnJ0fvv/+++vfvL0kqKChgUTygNnO5pIwMaelS93M1X9ATgH2ZDjdPPPGEJk+erPj4ePXs2VNJSUmS3KM4Xbt2tbxAADVAerp7qea+faXhw93P8fHudgDwM9O3X5Ck/Px85eXlqXPnzgoKcuejjz/+WBEREbrqqqssL9JK3H4BsFh6unuJ5vN/lZRdVcn9lQBYwMz3d6XCTU1GuAEsVHZzrXOXZj4XN9cCYBEz39+mT0sBgMemTRcPNpJ7NCcnx90PAPyEcAOg8vLyrO0HABYg3ACovJgYa/sBgAUINwAqLznZPafmnFuyeHE4pLg4dz8A8BPCDYDKczql2bPdfz4/4JS9njWLycQA/IpwA6Bq0tLcl3s3b+7dHhvLZeAAAqJOoAsAYANpadLgwe6rovLy3HNskpMZsQEQEIQbANZwOqWUlEBXAQCclgIAAPZCuAEAALbCaSnAKi4Xc04AoBog3ABWSE+XJkzwvhVBbKz7MmmuFgIAv+K0FFBVZXfFPv8eS7m57vb09MDUBQC1FOEGqAqXyz1iYxgXvlfWNnGiux8AwC8IN0BVcFdsAKh2CDdAVXBXbACodgg3QFVwV2wAqHYIN0BVcFdsAKh2CDdAVXBXbACodgg3QFVxV2wAqFZYxA+wAnfFBoBqg3ADWIW7YgNAtcBpKQAAYCuEGwAAYCuEGwAAYCuEGwAAYCuEGwAAYCuEGwAAYCuEGwAAYCuEGwAAYCuEGwAAYCuEGwAAYCuEGwAAYCuEGwAAYCuEGwAAYCvVItzMmTNH8fHxCg0NVc+ePfXxxx9ftO+iRYvkcDi8HqGhoX6sFuVxuaSMDGnpUvezyxXoigAAtVXAw80bb7yhSZMmaerUqcrMzFTnzp01YMAAFRQUXPQzERERysvL8zwOHTrkx4pxvvR0KT5e6ttXGj7c/Rwf724HAMDfAh5unn/+ed19990aNWqU2rdvr7lz5yosLEwLFiy46GccDoeio6M9j6ZNm16076lTp1RUVOT1gHXS06UhQ6TDh73bc3Pd7QQcAIC/BTTcnD59Wjt27FC/fv08bUFBQerXr5+2bt160c+dPHlSLVu2VFxcnAYPHqzPPvvson1nzJihyMhIzyMuLs7SY6jNXC5pwgTJMC58r6xt4kROUQEA/Cug4ebrr7+Wy+W6YOSladOmys/PL/czbdu21YIFC/TOO+/o73//u0pLS9W7d28dPn/o4H+mTJmiwsJCzyMnJ8fy46itNm26cMTmXIYh5eS4+wEA4C91Al2AWUlJSUpKSvK87t27t9q1a6e//vWveuqppy7oHxISopCQEH+WWGvk5VnbDwAAKwR05Oayyy6T0+nU0aNHvdqPHj2q6OjoCm2jbt266tq1q7766itflIhLiImq2PmmivYDAMAKAQ03wcHB6t69u9atW+dpKy0t1bp167xGZy7F5XJp9+7diomJ8VWZuIhkbVKscuRQabnvO1SqOGUrWZyXAgD4T8Cvlpo0aZLmzZunxYsXa+/evbrvvvtUXFysUaNGSZJGjBihKVOmePo/+eSTWrNmjQ4cOKDMzEz9v//3/3To0CGNGTMmUIdQazkL8jRbEyTpgoBT9nqWJspZwHkpAID/BHzOzdChQ3Xs2DE98cQTys/PV5cuXbR69WrPJOPs7GwFBf2Ywb777jvdfffdys/PV6NGjdS9e3dt2bJF7du3D9Qh1F4xMUrTCi3XEE3QbB3Wj1eixeqwZmmi0rRCinkggEUCAGobh2GUdyGvfRUVFSkyMlKFhYWKiIgIdDk1m8vlXq0vN1cuw6FNSlaeYhSjPCVrk5wOQ4qNlbKyJKcz0NUCAGowM9/fAR+5QQ3mdEqzZ0tDhsjpMJRibPjxPYfD/TxrFsEGAOBXAZ9zgxouLU1avlxq3ty7PTbW3Z6WFpi6AAC1FiM3qLq0NGnwYPdqfXl5UkyMlJzMiA0AICAIN7CG0ymlpAS6CgAAOC0FAADshXADAABshXADAABshXADAABshXADAABshXADAABshXADAABshXADAABshXADAABshXADAABshXADAABshXADAABshXADAABshXADAABshXADAABshXADAABshXADAABshXADAABshXADAABshXADAABshXADAABshXADAABshXADAABshXADAABshXADAABshXADAABshXADAABshXADAABshXADAABshXADAABshXADAABshXADAABshXADAABshXADAABshXADAABshXADAABshXADAABshXADAABshXADAABshXADAABshXADAABshXADAABshXADAABshXADAABshXADAABspU6gC7AL12mXNr20W3n7SxTTKkzJYxPlDHYGuiwAAGqdajFyM2fOHMXHxys0NFQ9e/bUxx9/fMn+b731lq666iqFhoYqMTFRq1at8lOl5Ut/ZJviw46q74NdNPwvvdX3wS6KDzuq9Ee2BbQuAABqo4CHmzfeeEOTJk3S1KlTlZmZqc6dO2vAgAEqKCgot/+WLVs0bNgwjR49Wjt37lRqaqpSU1O1Z88eP1fulv7INg159hoddkV7tee6ojXk2WsIOAAA+JnDMAwjkAX07NlTV199tf7yl79IkkpLSxUXF6f7779fjz766AX9hw4dquLiYr377ruetl69eqlLly6aO3fuT+6vqKhIkZGRKiwsVERERJVqd512KT7s6P+CzYU50aFSxTrzlFUSzSkqAACqwMz3d0BHbk6fPq0dO3aoX79+nragoCD169dPW7duLfczW7du9eovSQMGDLho/1OnTqmoqMjrYZVNL+3WYVczXeyv0VCQclzNteml3ZbtEwAAXFpAw83XX38tl8ulpk2berU3bdpU+fn55X4mPz/fVP8ZM2YoMjLS84iLi7OmeEl5+0ss7QcAAKou4HNufG3KlCkqLCz0PHJycizbdkyrMEv7AQCAqgtouLnsssvkdDp19OhRr/ajR48qOjq63M9ER0eb6h8SEqKIiAivh1WSxyYq1nlEDpWW+75DpYpz5ip5bKJl+wQAAJcW0HATHBys7t27a926dZ620tJSrVu3TklJSeV+Jikpyau/JK1du/ai/X3JGezU7EnZknRBwCl7PWtSDpOJAQDwo4Cflpo0aZLmzZunxYsXa+/evbrvvvtUXFysUaNGSZJGjBihKVOmePpPmDBBq1ev1syZM/X5559r2rRp2r59u8aPHx+Q+tP+1EvLH/5YzZ3ec35inXla/vDHSvtTr4DUBQBAbRXwFYqHDh2qY8eO6YknnlB+fr66dOmi1atXeyYNZ2dnKyjoxwzWu3dvLVmyRI8//rh+97vf6corr9Tbb7+tjh07BuoQlPanXhr8B5c2vbTrvBWKmwesJgAAaquAr3Pjb1aucwMAAPyjxqxzAwAAYDXCDQAAsBXCDQAAsBXCDQAAsBXCDQAAsBXCDQAAsBXCDQAAsBXCDQAAsBXCDQAAsJWA337B38oWZC4qKgpwJQAAoKLKvrcrcmOFWhduTpw4IUmKi4sLcCUAAMCsEydOKDIy8pJ9at29pUpLS3XkyBGFh4fL4XBYuu2ioiLFxcUpJyenVty3iuO1N47X3mrb8Uq175jtdryGYejEiRNq1qyZ1w21y1PrRm6CgoIUGxvr031ERETY4j+kiuJ47Y3jtbfadrxS7TtmOx3vT43YlGFCMQAAsBXCDQAAsBXCjYVCQkI0depUhYSEBLoUv+B47Y3jtbfadrxS7Tvm2na856p1E4oBAIC9MXIDAABshXADAABshXADAABshXADAABshXBj0pw5cxQfH6/Q0FD17NlTH3/88SX7v/XWW7rqqqsUGhqqxMRErVq1yk+VWsPM8c6bN0/Jyclq1KiRGjVqpH79+v3k3091Y/bnW2bZsmVyOBxKTU31bYEWM3u8x48f17hx4xQTE6OQkBC1adOmRv03bfZ4Z82apbZt26pevXqKi4vTgw8+qB9++MFP1VbNxo0bNWjQIDVr1kwOh0Nvv/32T34mIyND3bp1U0hIiFq3bq1Fixb5vE6rmD3e9PR03Xjjjbr88ssVERGhpKQkvf/++/4p1gKV+fmW2bx5s+rUqaMuXbr4rL5AI9yY8MYbb2jSpEmaOnWqMjMz1blzZw0YMEAFBQXl9t+yZYuGDRum0aNHa+fOnUpNTVVqaqr27Nnj58orx+zxZmRkaNiwYVq/fr22bt2quLg49e/fX7m5uX6uvHLMHm+ZgwcPavLkyUpOTvZTpdYwe7ynT5/WjTfeqIMHD2r58uXat2+f5s2bp+bNm/u58soxe7xLlizRo48+qqlTp2rv3r3629/+pjfeeEO/+93v/Fx55RQXF6tz586aM2dOhfpnZWVp4MCB6tu3r3bt2qWJEydqzJgxNeYL3+zxbty4UTfeeKNWrVqlHTt2qG/fvho0aJB27tzp40qtYfZ4yxw/flwjRozQDTfc4KPKqgkDFXbNNdcY48aN87x2uVxGs2bNjBkzZpTb/7bbbjMGDhzo1dazZ0/jN7/5jU/rtIrZ4z3f2bNnjfDwcGPx4sW+KtFSlTnes2fPGr179zbmz59v3HnnncbgwYP9UKk1zB7vyy+/bFxxxRXG6dOn/VWipcwe77hx44zrr7/eq23SpElGnz59fFqnL0gyVqxYcck+jzzyiNGhQwevtqFDhxoDBgzwYWW+UZHjLU/79u2N6dOnW1+Qj5k53qFDhxqPP/64MXXqVKNz584+rSuQGLmpoNOnT2vHjh3q16+fpy0oKEj9+vXT1q1by/3M1q1bvfpL0oABAy7avzqpzPGer6SkRGfOnFHjxo19VaZlKnu8Tz75pKKiojR69Gh/lGmZyhzvP//5TyUlJWncuHFq2rSpOnbsqKeffloul8tfZVdaZY63d+/e2rFjh+fU1YEDB7Rq1SrdcsstfqnZ32ry7ysrlJaW6sSJEzXi91VlLVy4UAcOHNDUqVMDXYrP1bobZ1bW119/LZfLpaZNm3q1N23aVJ9//nm5n8nPzy+3f35+vs/qtEpljvd8v/3tb9WsWbMLfmFWR5U53n//+9/629/+pl27dvmhQmtV5ngPHDigDz/8ULfffrtWrVqlr776SmPHjtWZM2eq/S/Lyhzv8OHD9fXXX+vaa6+VYRg6e/as7r333hpzWsqsi/2+Kioq0vfff6969eoFqDL/eO6553Ty5EnddtttgS7FJ7788ks9+uij2rRpk+rUsf9XPyM38IlnnnlGy5Yt04oVKxQaGhrocix34sQJ3XHHHZo3b54uu+yyQJfjF6WlpYqKitIrr7yi7t27a+jQoXrsscc0d+7cQJfmExkZGXr66af10ksvKTMzU+np6Vq5cqWeeuqpQJcGiy1ZskTTp0/Xm2++qaioqECXYzmXy6Xhw4dr+vTpatOmTaDL8Qv7xzeLXHbZZXI6nTp69KhX+9GjRxUdHV3uZ6Kjo031r04qc7xlnnvuOT3zzDP64IMP1KlTJ1+WaRmzx7t//34dPHhQgwYN8rSVlpZKkurUqaN9+/apVatWvi26Cirz842JiVHdunXldDo9be3atVN+fr5Onz6t4OBgn9ZcFZU53t///ve64447NGbMGElSYmKiiouLdc899+ixxx5TUJC9/m14sd9XERERth61WbZsmcaMGaO33nqrRowyV8aJEye0fft27dy5U+PHj5fk/n1lGIbq1KmjNWvW6Prrrw9wlday1/+dPhQcHKzu3btr3bp1nrbS0lKtW7dOSUlJ5X4mKSnJq78krV279qL9q5PKHK8k/elPf9JTTz2l1atXq0ePHv4o1RJmj/eqq67S7t27tWvXLs/jF7/4hedKk7i4OH+Wb1plfr59+vTRV1995QlxkvTFF18oJiamWgcbqXLHW1JSckGAKQt2hg1vyVeTf19V1tKlSzVq1CgtXbpUAwcODHQ5PhMREXHB76t7771Xbdu21a5du9SzZ89Al2i9AE9orlGWLVtmhISEGIsWLTL++9//Gvfcc4/RsGFDIz8/3zAMw7jjjjuMRx991NN/8+bNRp06dYznnnvO2Lt3rzF16lSjbt26xu7duwN1CKaYPd5nnnnGCA4ONpYvX27k5eV5HidOnAjUIZhi9njPV9OuljJ7vNnZ2UZ4eLgxfvx4Y9++fca7775rREVFGX/4wx8CdQimmD3eqVOnGuHh4cbSpUuNAwcOGGvWrDFatWpl3HbbbYE6BFNOnDhh7Ny509i5c6chyXj++eeNnTt3GocOHTIMwzAeffRR44477vD0P3DggBEWFmY8/PDDxt69e405c+YYTqfTWL16daAOwRSzx/v6668bderUMebMmeP1++r48eOBOgRTzB7v+ex+tRThxqQXX3zRaNGihREcHGxcc801xrZt2zzvXXfddcadd97p1f/NN9802rRpYwQHBxsdOnQwVq5c6eeKq8bM8bZs2dKQdMFj6tSp/i+8ksz+fM9V08KNYZg/3i1bthg9e/Y0QkJCjCuuuML44x//aJw9e9bPVVeemeM9c+aMMW3aNKNVq1ZGaGioERcXZ4wdO9b47rvv/F94Jaxfv77c/x/LjvHOO+80rrvuugs+06VLFyM4ONi44oorjIULF/q97soye7zXXXfdJftXd5X5+Z7L7uHGYRg2HF8FAAC1FnNuAACArRBuAACArRBuAACArRBuAACArRBuAACArRBuAACArRBuAACArRBuAACArRBuAPhcRkaGHA6Hjh8/7vd9OxwOORwONWzY0Kf7SUlJ8exr165dPt0XgEsj3ACwVEpKiiZOnOjV1rt3b+Xl5SkyMjIgNS1cuFBffPGFT/eRnp6ujz/+2Kf7AFAxdQJdAAD7Cw4OVnR0dMD237BhQ0VFRfl0H40bN1ZRUZFP9wGgYhi5AWCZkSNHasOGDZo9e7bnFM3BgwcvOC21aNEiNWzYUO+++67atm2rsLAwDRkyRCUlJVq8eLHi4+PVqFEjPfDAA3K5XJ7tnzp1SpMnT1bz5s1Vv3599ezZUxkZGabrnDZtmrp06aIFCxaoRYsWatCggcaOHSuXy6U//elPio6OVlRUlP74xz96PmMYhqZNm6YWLVooJCREzZo10wMPPFDVvzIAPsDIDQDLzJ49W1988YU6duyoJ598UpJ0+eWX6+DBgxf0LSkp0QsvvKBly5bpxIkTSktL06233qqGDRtq1apVOnDggH75y1+qT58+Gjp0qCRp/Pjx+u9//6tly5apWbNmWrFihW666Sbt3r1bV155pala9+/fr/fee0+rV6/W/v37NWTIEB04cEBt2rTRhg0btGXLFt11113q16+fevbsqX/84x/685//rGXLlqlDhw7Kz8/Xp59+WuW/MwDWI9wAsExkZKSCg4MVFhb2k6ehzpw5o5dfflmtWrWSJA0ZMkSvvfaajh49qgYNGqh9+/bq27ev1q9fr6FDhyo7O1sLFy5Udna2mjVrJkmaPHmyVq9erYULF+rpp582VWtpaakWLFig8PBwz7727dunVatWKSgoSG3bttX//d//af369erZs6eys7MVHR2tfv36qW7dumrRooWuueaayv1FAfApwg2AgAgLC/MEG0lq2rSp4uPj1aBBA6+2goICSdLu3bvlcrnUpk0br+2cOnVKTZo0Mb3/+Ph4hYeHe+3L6XQqKCjIq61s/7/61a80a9YsXXHFFbrpppt0yy23aNCgQapTh1+jQHXD/5UAAqJu3bperx0OR7ltpaWlkqSTJ0/K6XRqx44dcjqdXv3ODUS+2n9cXJz27dunDz74QGvXrtXYsWP17LPPasOGDRd8DkBgEW4AWCo4ONhrErBVunbtKpfLpYKCAiUnJ1u+/YqoV6+eBg0apEGDBmncuHG66qqrtHv3bnXr1i0g9QAoH+EGgKXi4+P10Ucf6eDBg2rQoIEaN25syXbbtGmj22+/XSNGjNDMmTPVtWtXHTt2TOvWrVOnTp00cOBAS/ZzMYsWLZLL5VLPnj0VFhamv//976pXr55atmzp0/0CMI9LwQFYavLkyXI6nWrfvr0uv/xyZWdnW7bthQsXasSIEXrooYfUtm1bpaam6pNPPlGLFi0s28fFNGzYUPPmzVOfPn3UqVMnffDBB/rXv/5Vqfk+AHzLYRiGEegiAMBXHA6HVqxYodTUVJ/v6+DBg0pISNDOnTvVpUsXn+8PQPkYuQFge8OGDVNsbKxP93HzzTerQ4cOPt0HgIph5AaArX311VeSJKfTqYSEBJ/tJzc3V99//70kqUWLFgoODvbZvgBcGuEGAADYCqelAACArRBuAACArRBuAACArRBuAACArRBuAACArRBuAACArRBuAACArRBuAACArfx/LtKCYJgb4zAAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "f, ax = plt.subplots()\n",
    "ax.plot(times * 1e3, stds_spin, \"ro\", label=\"spin\")\n",
    "ax.plot(times * 1e3, stds_density, \"bo\", label=\"density\")\n",
    "ax.set_xlabel(\"time [ms]\")\n",
    "ax.set_ylabel(\"standard deviation\")\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is exactly the expected behavior for the decoupling of spin and density. The density spreads much faster than the spin. This is the spin-charge separation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The ffsim backend\n",
    "\n",
    "As we can see the time of the simulation scales exponentially with the number of sites. To increase the possible we can also run the simulations with the `ffsim` backend, which is exactly optimized for such tasks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Duration for of the two jobs: 62.67665886878967 seconds\n"
     ]
    }
   ],
   "source": [
    "# prepare the jobs\n",
    "jobs_1_ffsim = backend_ffsim.run(circuits_1, shots=Nshots, num_species=2)\n",
    "jobs_2_ffsim = backend_ffsim.run(circuits_2, shots=Nshots, num_species=2)\n",
    "\n",
    "# run them and time it\n",
    "start_time = time.time()\n",
    "count_list_1 = jobs_1_ffsim.result().get_counts()\n",
    "count_list_2 = jobs_2_ffsim.result().get_counts()\n",
    "end_time = time.time()\n",
    "print(f\"Duration for of the two jobs: {end_time - start_time} seconds\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We observe here that the code might not run quite as fast as we might like it. This is because the `ffsim` backend is optimized for trotterized approaches. Let us therefore try it out. Most important is to find time steps that are small compared to the large energy scale in the system, which is $U$ in this case. So let us run one long time evolution in two different ways:\n",
    "\n",
    "1. With the full Fermi Hubbard Hamiltonian applied as a single gate.\n",
    "2. With the trotterized version of the Fermi Hubbard Hamiltonian."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "hop_long = Parameter(\"hop_long\")\n",
    "int_long = Parameter(\"int_long\")\n",
    "# initialize one spin-up and one spin-down atom in the left tweezer\n",
    "qc_long = backend_ffsim.initialize_circuit(init_state_1)\n",
    "\n",
    "# apply a global Fermi-Hubbard gate with interaction\n",
    "qc_long.fhubbard(\n",
    "    j=hop_long * np.ones(Nsites - 1), u=int_long, mu=np.zeros(Nsites), modes=range(2 * Nsites)\n",
    ")\n",
    "qc_long.measure_all()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Duration for of the two jobs: 1.97 seconds\n"
     ]
    }
   ],
   "source": [
    "tmax = 1 / U\n",
    "\n",
    "circuit_long = qc_long.bind_parameters({hop_long: J * tmax, int_long: U * tmax})\n",
    "job_long = backend.run(circuit_long, shots=5, num_species=2)\n",
    "start_time = time.time()\n",
    "count_long = job_long.result().get_counts()\n",
    "end_time = time.time()\n",
    "print(f\"Duration for of the two jobs: {end_time - start_time:.2f} seconds\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And now we will trotterize the circuit to make it more efficient."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "time step = 0.01 ms\n",
      "Number of trotter steps: 4\n"
     ]
    }
   ],
   "source": [
    "time_step = 1/(4*U)\n",
    "print(f\"time step = {time_step*1e3:.2f} ms\")\n",
    "\n",
    "Ntrott = int(tmax / time_step)\n",
    "print(f\"Number of trotter steps: {Ntrott}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "trott_hop = Parameter(\"trott_hop\")\n",
    "trott_int = Parameter(\"trott_int\")\n",
    "# initialize one spin-up and one spin-down atom in the left tweezer\n",
    "qc_trott = backend_ffsim.initialize_circuit(init_state_1)\n",
    "\n",
    "for ii in range(Ntrott):\n",
    "    # apply a global Fermi-Hubbard gate with interaction\n",
    "    qc_trott.fhop(\n",
    "        j=trott_hop * np.ones(Nsites - 1),\n",
    "        modes=range(2 * Nsites),\n",
    "    )\n",
    "    qc_trott.fint(u=trott_int, modes=range(2 * Nsites))\n",
    "# apply a global Fermi-Hubbard gate with interaction\n",
    "qc_trott.measure_all()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Duration for of the two jobs: 0.52 ms\n"
     ]
    }
   ],
   "source": [
    "circuit_trott = qc_trott.bind_parameters({trott_hop: J * time_step, trott_int: U * time_step})\n",
    "job_trott = backend_ffsim.run(circuit_trott, shots=5, num_species=2)\n",
    "start_time = time.time()\n",
    "count_long = job_long.result().get_counts()\n",
    "end_time = time.time()\n",
    "print(f\"Duration for of the two jobs: {(end_time - start_time)*1e3:.2f} ms\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So as you can see nicely here the trotterized version is much faster and will enable us to simulate much larger systems in the future."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Conclusion\n",
    "\n",
    "In this notebook we have seen how to simulate the dynamics of a fermionic system in a 1D lattice. We have seen how to initialize the system in an anti-ferromagnetic configuration and how to measure the spin and density dynamics. We have seen how the spin and density dynamics decouple and how this can be quantified.\n",
    "\n",
    "Finally, we have seen that the `ffsim` backend is much faster for trotterized approaches and how it might potentially be used for larger systems in the future."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "[1] [Jayadev Vijayan et al.](https://arxiv.org/abs/1905.13638) *Time-Resolved Observation of Spin-Charge Deconfinement in Fermionic Hubbard Chains*. Science 367, 186 (2020)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<h3>Version Information</h3><table><tr><th>Qiskit Software</th><th>Version</th></tr><tr><td><code>qiskit-terra</code></td><td>0.24.0</td></tr><tr><td><code>qiskit-aer</code></td><td>0.12.2</td></tr><tr><td><code>qiskit-nature</code></td><td>0.6.2</td></tr><tr><th>System information</th></tr><tr><td>Python version</td><td>3.10.0</td></tr><tr><td>Python compiler</td><td>Clang 12.0.0 </td></tr><tr><td>Python build</td><td>default, Mar  3 2022 03:54:28</td></tr><tr><td>OS</td><td>Darwin</td></tr><tr><td>CPUs</td><td>8</td></tr><tr><td>Memory (Gb)</td><td>8.0</td></tr><tr><td colspan='2'>Fri Jan 12 11:17:23 2024 CET</td></tr></table>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import qiskit.tools.jupyter\n",
    "\n",
    "%qiskit_version_table"
   ]
  }
 ],
 "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.10.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}