{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "toc": true
   },
   "source": [
    "<h1>Table of Contents<span class=\"tocSkip\"></span></h1>\n",
    "<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#Path-Integrals\" data-toc-modified-id=\"Path-Integrals-1\"><span class=\"toc-item-num\">1&nbsp;&nbsp;</span>Path Integrals</a></span><ul class=\"toc-item\"><li><span><a href=\"#Free-Particle\" data-toc-modified-id=\"Free-Particle-1.1\"><span class=\"toc-item-num\">1.1&nbsp;&nbsp;</span>Free Particle</a></span></li></ul></li><li><span><a href=\"#Discretization\" data-toc-modified-id=\"Discretization-2\"><span class=\"toc-item-num\">2&nbsp;&nbsp;</span>Discretization</a></span><ul class=\"toc-item\"><li><span><a href=\"#Discrete-Unitarity\" data-toc-modified-id=\"Discrete-Unitarity-2.1\"><span class=\"toc-item-num\">2.1&nbsp;&nbsp;</span>Discrete Unitarity</a></span></li><li><span><a href=\"#Free-Particle-on-a-Periodic-Lattice\" data-toc-modified-id=\"Free-Particle-on-a-Periodic-Lattice-2.2\"><span class=\"toc-item-num\">2.2&nbsp;&nbsp;</span>Free Particle on a Periodic Lattice</a></span></li></ul></li><li><span><a href=\"#Harmonic-Oscillator\" data-toc-modified-id=\"Harmonic-Oscillator-3\"><span class=\"toc-item-num\">3&nbsp;&nbsp;</span>Harmonic Oscillator</a></span></li><li><span><a href=\"#Harmonic-Oscillator\" data-toc-modified-id=\"Harmonic-Oscillator-4\"><span class=\"toc-item-num\">4&nbsp;&nbsp;</span>Harmonic Oscillator</a></span></li></ul></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Path Integrals\n",
    "\n",
    "Based on a discussion with Fred Gittes, we compute the propagator for small systems using a path integral.\n",
    "\n",
    "<!-- END_TEASER -->"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>.grade {\n",
       "   background-color: #66FFCC;\n",
       "}\n",
       "</style>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<script id=\"MathJax-Element-48\" type=\"math/tex\">\\newcommand{\\vect}[1]{\\mathbf{#1}}\n",
       "\\newcommand{\\uvect}[1]{\\hat{#1}}\n",
       "\\newcommand{\\abs}[1]{\\lvert#1\\rvert}\n",
       "\\newcommand{\\norm}[1]{\\lVert#1\\rVert}\n",
       "\\newcommand{\\I}{\\mathrm{i}}\n",
       "\\newcommand{\\ket}[1]{\\left|#1\\right\\rangle}\n",
       "\\newcommand{\\bra}[1]{\\left\\langle#1\\right|}\n",
       "\\newcommand{\\braket}[1]{\\langle#1\\rangle}\n",
       "\\newcommand{\\Braket}[1]{\\left\\langle#1\\right\\rangle}\n",
       "\\newcommand{\\op}[1]{\\mathbf{#1}}\n",
       "\\newcommand{\\mat}[1]{\\mathbf{#1}}\n",
       "\\newcommand{\\d}{\\mathrm{d}}\n",
       "\\newcommand{\\pdiff}[3][]{\\frac{\\partial^{#1} #2}{\\partial {#3}^{#1}}}\n",
       "\\newcommand{\\diff}[3][]{\\frac{\\d^{#1} #2}{\\d {#3}^{#1}}}\n",
       "\\newcommand{\\ddiff}[3][]{\\frac{\\delta^{#1} #2}{\\delta {#3}^{#1}}}\n",
       "\\newcommand{\\floor}[1]{\\left\\lfloor#1\\right\\rfloor}\n",
       "\\newcommand{\\ceil}[1]{\\left\\lceil#1\\right\\rceil}\n",
       "\\DeclareMathOperator{\\Tr}{Tr}\n",
       "\\DeclareMathOperator{\\erf}{erf}\n",
       "\\DeclareMathOperator{\\erfi}{erfi}\n",
       "\\DeclareMathOperator{\\sech}{sech}\n",
       "\\DeclareMathOperator{\\sgn}{sgn}\n",
       "\\DeclareMathOperator{\\order}{O}\n",
       "\\DeclareMathOperator{\\diag}{diag}\n",
       "\n",
       "\\newcommand{\\mylabel}[1]{\\label{#1}\\tag{#1}}\n",
       "\\newcommand{\\degree}{\\circ}</script>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import mmf_setup;mmf_setup.nbinit(quiet=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One of Feynman great insights was that the propagator for quantum systems can be constructed through a path integral.  He expressed this by saying that the probability for a particle starting at point and time $(x_i, t_i)$ to be later at a point and time $(x_f, t_f)$ could be obtained by averaging the phases $e^{\\I S[x]}$ over all possible trajectories $x(t)$ that the particle might take from $x(t_i) = x_i$ to $x(t_f) = x_f$ where $S[x]$ is the classical action:\n",
    "\n",
    "$$\n",
    "  S[x] = \\int_{t_i}^{t_f} L(x, \\dot{x}, t) \\d{t}, \\qquad\n",
    "  L(x, \\dot{x}(t), t) = \\frac{m \\dot{x}^2}{2} - V(x, t),\\qquad\n",
    "  x(t_i) = x_i, \\qquad x(t_f) = x_f.\n",
    "$$\n",
    "\n",
    "Formally, this can be expressed in terms of a new type of integral:\n",
    "\n",
    "$$\n",
    "  \\newcommand{\\D}[1]{\\mathcal{D}[#1]\\;}\n",
    "  U(x_f, t_f;x_i, t_i) \n",
    "  = \\braket{x_f|\\op{U}(t_f,t_i|x_i}\n",
    "  = \\int \\D{x}e^{\\I S[x]/\\hbar},\n",
    "$$\n",
    "\n",
    "where $\\op{U}(t_f, t_i)$ is the quantum mechanical propagator:\n",
    "\n",
    "$$\n",
    "  \\ket{\\psi(t_f)} = \\op{U}(t_f, t_i)\\ket{\\psi(t_i)}, \\qquad\n",
    "  \\op{U}(t_f, t_i) = \\mathcal{T}\\exp\\left(\n",
    "    \\int_{t_i}^{t_f}\\frac{\\op{H}(t)}{\\I\\hbar}\\d{t}\n",
    "  \\right), \\qquad\n",
    "  \\op{H} = \\frac{\\op{P}^2}{2m} + V(\\op{X}, t).\n",
    "$$\n",
    "\n",
    "In the expression for $\\op{U}(t_f, t_i)$, the integral must be time-ordered as signified by the operator $\\mathcal{T}$ which means that in every term, operators must appear in descending time order.  Formally this can be understood by Taylor expanding the exponential and then manually rearranging all terms so that they are in the correct order.\n",
    "\n",
    "For example, expanding only to quadratic order, we would have:\n",
    "\n",
    "$$\n",
    "  \\op{U}(t_f, t_i) = \\mat{1} \n",
    "  +\n",
    "  \\int_{t_i}^{t_f}\\frac{\\op{H}(t)}{\\I\\hbar}\\d{t}\n",
    "  +\n",
    "  \\frac{1}{2!}\n",
    "  \\mathcal{T}\\left(\n",
    "    \\int_{t_i}^{t_f}\\frac{\\op{H}(t)}{\\I\\hbar}\\d{t}\n",
    "  \\right)^2\n",
    "  +\n",
    "  \\cdots =\\\\\n",
    "  \\mathcal{T}\\left(\n",
    "    \\int_{t_i}^{t_f}\\frac{\\op{H}(t)}{\\I\\hbar}\\d{t}\n",
    "  \\right)^2\n",
    "  =\n",
    "  \\frac{1}{(\\I\\hbar)^{2}}\n",
    "  \\int_{t_i}^{t_f}\\d{t_1}\n",
    "  \\int_{t_i}^{t_f}\\d{t_2}\\;\n",
    "  \\mathcal{T}\n",
    "  \\op{H}(t_1)\n",
    "  \\op{H}(t_2)\n",
    "  =\\\\\n",
    "  =\n",
    "  \\frac{1}{(\\I\\hbar)^{2}}\n",
    "  \\int_{t_i}^{t_f}\\d{t_1}\n",
    "  \\int_{t_i}^{t_f}\\d{t_2}\\;\n",
    "  \\begin{cases}\n",
    "      \\op{H}(t_1)\\op{H}(t_2) & t_2\\leq t_1\\\\\n",
    "      \\op{H}(t_2)\\op{H}(t_1) & t_1<t_2\n",
    "  \\end{cases}=\\\\\n",
    "  =\n",
    "  \\frac{1}{(\\I\\hbar)^{2}}\n",
    "  \\int_{t_i}^{t_f}\\d{t_1}\n",
    "  \\left[\n",
    "      \\int_{t_i}^{t_1}\\d{t_2}\\;\n",
    "      \\op{H}(t_1)\\op{H}(t_2)\n",
    "      +\n",
    "      \\int_{t_1}^{t_f}\\d{t_2}\\;\n",
    "      \\op{H}(t_2)\\op{H}(t_1)\n",
    "  \\right].\n",
    "$$\n",
    "\n",
    "*Of course, if the Hamiltonian commutes at different times $[\\op{H}(t_1), \\op{H}(t_2)] = 0$ – for example, if $V(\\op{X}, t) = (\\op{X})$ is time independent – then there is no need to worry about the time-ordering and one can simply ignore this complication.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "Consider the path integral representation of the propagator corresponding to a time-independent Hamiltonian:\n",
    "\n",
    "$$\n",
    "  \\newcommand{\\D}[1]{\\mathcal{D}[#1]\\;}\n",
    "  U(x_N, T;x_0, 0) = \\int \\D{x}e^{\\I S/\\hbar}, \\qquad\n",
    "  S = \\int_{0}^{T}\\d{t} L[x(t), \\dot{x}(t)], \\qquad\n",
    "  L[x, \\dot{x}] = \\frac{m}{2}\\dot{x}^2 - V(x), \\qquad\n",
    "  H[x, p] = \\frac{p^2}{2m} + V(x).\n",
    "$$\n",
    "\n",
    "Discretizing in time we have\n",
    "\n",
    "$$\n",
    "  S_N = \\sum_{i=0}^{N-1} \\d{t}\\;\\left[\n",
    "    \\frac{m}{2}\\left(\\frac{x_{i+1}-x_{i}}{\\d{t}}\\right)^2\n",
    "    - \\frac{V(x_{i}) + V(x_{i+1})}{2}\n",
    "  \\right], \\qquad\n",
    "  \\d{t} = \\frac{T}{N}.\\\\\n",
    "  U(x_N, T;x_0, 0) = \\lim_{N\\rightarrow \\infty} A^N \\int\\d{x_1}\\int\\d{x_2}\\cdots\\int\\d{x_{N-1}}\n",
    "  \\exp\\left\\{\n",
    "    \\frac{\\I}{\\hbar}\\sum_{i=0}^{N-1}\\d{t}\n",
    "    \\left[\n",
    "      \\frac{m}{2}\\frac{(x_{i+1}-x_{i})^2}{\\d{t}^2}\n",
    "      - \n",
    "      \\frac{V(x_{i}) + V(x_{i+1})}{2}\n",
    "    \\right]\n",
    "  \\right\\}\n",
    "    = \\\\\n",
    "    = \\lim_{N\\rightarrow \\infty} A^N \\int\\d{x_1}\\int\\d{x_2}\\cdots\\int\\d{x_{N-1}}\n",
    "    \\prod_{i=0}^{N-1}\n",
    "  \\exp\\left\\{\n",
    "    \\frac{\\I}{\\hbar}\n",
    "    \\left[\n",
    "      \\frac{m}{2}\\frac{(x_{i+1}-x_{i})^2}{\\d{t}}\n",
    "      - \n",
    "      \\frac{V(x_{i}) + V(x_{i+1})}{2}\\d{t}\n",
    "    \\right]\n",
    "  \\right\\} \n",
    "    = \\\\\n",
    "    = \\lim_{N\\rightarrow \\infty} A^N\n",
    "    \\int\\d{x_1}\\int\\d{x_2}\\cdots\\int\\d{x_{N-1}}\n",
    "    M_{x_Nx_{N-1}}M_{x_{N-1}x_{N-1}}\\cdots M_{x_2x_1}M_{x_1x_0}\n",
    "    = \\lim_{N\\rightarrow \\infty} (A \\mat{M})^{N}, \\\\\n",
    "    M_{xy} = \\exp\\left\\{\n",
    "      \\frac{\\I}{\\hbar}\n",
    "      \\left[\n",
    "        \\frac{m}{2}\\frac{(x-y)^2}{\\d{t}}\n",
    "        - \n",
    "        \\frac{V(x) + V(y)}{2}\\d{t}\n",
    "      \\right]\n",
    "    \\right\\}.  \n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The matrices $A\\mat{M}$ are unitary $\\abs{A}^2\\mat{M}^\\dagger\\cdot\\mat{M} = \\mat{1}$ if properly normalized:\n",
    "\n",
    "$$\n",
    "  [\\mat{M}^\\dagger\\cdot\\mat{M}]_{xy} = \\int\\d{z}\\; M^*_{zx}M_{zy}\n",
    "  = \\int \\d{z}\n",
    "    \\exp\\left\\{\n",
    "      \\frac{\\I}{\\hbar}\n",
    "      \\left[\n",
    "        \\frac{m}{2}\\frac{(z-y)^2 - (z-x)^2}{\\d{t}}\n",
    "        + \n",
    "        \\frac{V(x) - V(y)}{2}\\d{t}\n",
    "      \\right]\n",
    "    \\right\\}=\\\\\n",
    "  = \\exp\\left\\{\n",
    "      \\frac{\\I}{\\hbar}\n",
    "      \\left[\n",
    "        \\frac{m}{2}\\frac{y^2 - x^2}{\\d{t}}\n",
    "        + \n",
    "        \\frac{V(x) - V(y)}{2}\\d{t}\n",
    "      \\right]\n",
    "    \\right\\}\n",
    "  \\int \\d{z}\n",
    "    \\exp\\left\\{\n",
    "      \\frac{\\I}{\\hbar}\n",
    "      \\left[\n",
    "        \\frac{m}{2}\\frac{2z(x-y)}{\\d{t}}\n",
    "      \\right]\n",
    "    \\right\\}=\\\\\n",
    "  =\n",
    "  \\frac{2\\pi\\hbar\\d{t}}{m}\n",
    "  \\delta(x-y).\n",
    "$$\n",
    "\n",
    "Hence (the factor of $\\I$ included in $A$ will help later):\n",
    "\n",
    "$$\n",
    "  A = \\sqrt{\\frac{m}{2\\pi \\I\\hbar\\d{t}}}, \\qquad\n",
    "  AM_{xy} = \\sqrt{\\frac{m}{2\\pi \\I\\hbar\\d{t}}}\\exp\\left\\{\n",
    "      \\frac{\\I}{\\hbar}\n",
    "      \\left[\n",
    "        \\frac{m}{2}\\frac{(x-y)^2}{\\d{t}}\n",
    "        - \n",
    "        \\frac{V(x) + V(y)}{2}\\d{t}\n",
    "      \\right]\n",
    "    \\right\\}.  \n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Free Particle\n",
    "\n",
    "If we neglect the potential $V=0$, then we can explicitly perform the matrix multiplication, completing the square in $z$ and absorbing the integrating factor into $A$:\n",
    "\n",
    "$$\n",
    "  [(A\\mat{M})^2]_{xy} = \\frac{m}{2\\pi \\I\\hbar\\d{t}}\\int \\d{z}\\; \n",
    "  \\exp\\left\\{\\frac{\\I}{\\hbar}\\frac{m}{2}\\frac{(x-z)^2 + (z-y)^2}{\\d{t}}\\right\\}\n",
    "   =\\\\\n",
    "   =\\exp\\left\\{\\frac{\\I}{2\\hbar}\\frac{m}{2}\\frac{(x-y)^2}{\\d{t}}\\right\\}\n",
    "  \\frac{m}{2\\pi \\I\\hbar\\d{t}}\\int \\d{z}\\;\n",
    "  \\exp\\left\\{\\frac{m\\I}{\\hbar}\\frac{z^2}{\\d{t}}\\right\\}\n",
    "   =\n",
    "  \\sqrt{\\frac{m}{2\\pi \\I\\hbar 2\\d{t}}}\n",
    "  \\exp\\left\\{\\frac{\\I}{2\\hbar}\\frac{m}{2}\\frac{(x-y)^2}{\\d{t}}\\right\\}\n",
    "$$\n",
    "\n",
    "Note: the integral over $z$ is a little dubious since the factor of $\\I$ makes it non-convergent.\n",
    "\n",
    "Likewise:\n",
    "\n",
    "$$\n",
    "  [(A\\mat{M})^N]_{xy} \n",
    "  =\n",
    "  \\sqrt{\\frac{m}{2\\pi \\I\\hbar N\\d{t}}}\n",
    "  \\exp\\left\\{\\frac{\\I}{N\\hbar}\\frac{m}{2}\\frac{(x-y)^2}{\\d{t}}\\right\\}\n",
    "  = \n",
    "  \\sqrt{\\frac{m}{2\\pi \\I\\hbar T}}\n",
    "  \\exp\\left\\{\\frac{\\I}{\\hbar}\\frac{m}{2}\\frac{(x-y)^2}{T}\\right\\}.\n",
    "$$\n",
    "\n",
    "Hence (including the normalization factor),\n",
    "\n",
    "$$\n",
    "  U(y,T;x, 0) = \\sqrt{\\frac{m}{2\\pi \\hbar \\I T}}\\exp\\left\\{\\frac{\\I m (y-x)^2}{2\\hbar T}\\right\\}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Discretization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To make the path integral concrete (and well defined) we can express it in terms of a discrete lattice in space and time.  A useful approach is to use periodic functions on a lattice of $N_x$ points in a box of length $L$.  With this discretization both the positions $x=x_n$ and wave-numbers $k = p/\\hbar = k_n$ become discrete with $n$ running over a finite range $n \\in \\{0, 1, 2\\cdots N-1\\}$.\n",
    "\n",
    "$$\n",
    "  x_n = n\\d{x} = n\\frac{L}{N_x}, \\qquad\n",
    "  k_n = \\frac{2\\pi n}{L},\n",
    "  \\qquad n\\in\\{-N/2, -N/2+1, \\dots N/2-1\\}\n",
    "$$\n",
    "\n",
    "Two sets of normalizations are useful and follow from the following lattice relationship\n",
    "\n",
    "$$\n",
    "  \\sum_{l}e^{\\I k_l(x_m - x_n)} = \\sum_{l}e^{\\I 2\\pi l(m - n)/N} = N\\delta_{nm}, \\\\\n",
    "  \\frac{1}{L}\\sum_{l}e^{\\I k_l(x_m - x_n)} = \\frac{\\delta_{nm}}{\\d{x}}, \\qquad\n",
    "  \\d{x}\\sum_{l}e^{\\I k_l(x_m - x_n)} = L\\delta_{nm}.\n",
    "$$\n",
    "\n",
    "**Physical**:\n",
    "\n",
    "$$\n",
    "  \\braket{x_m|k_n} = e^{\\I k_n x_m}, \\qquad\n",
    "  \\braket{x_m|x_n} = \\frac{\\delta_{nm}}{\\d{x}}, \\qquad\n",
    "  \\braket{k_m|k_n} = L\\delta_{nm},\\\\\n",
    "  \\mat{1} = \\d{x}\\sum_{n}\\ket{x_n}\\bra{x_n} = \\frac{1}{L}\\sum_{m}\\ket{k_n}\\bra{k_n}.\n",
    "$$\n",
    "\n",
    "In this normalization, continuum equations can be directly translated by replacing integrals as follows:\n",
    "\n",
    "$$\n",
    "  \\int \\d{x} \\rightarrow \\d{x}\\sum_{n}, \\qquad\n",
    "  \\int \\frac{\\d{k}}{(2\\pi)} \\rightarrow \\frac{1}{L}\\sum_{n}.\n",
    "$$\n",
    "\n",
    "**Numerical**:\n",
    "\n",
    "The **Numerical** normalization is obtained from the **Physical** normalization by rescaling:\n",
    "\n",
    "$$\n",
    "  \\ket{x_n} \\rightarrow \\frac{\\ket{x_n}}{\\sqrt{\\d{x}}}, \\qquad\n",
    "  \\ket{k_n} \\rightarrow \\sqrt{L}\\ket{k_n}\\\\\n",
    "  \\braket{x_m|k_n} = \\sqrt{\\frac{\\d{x}}{L}} e^{\\I k_n x_m} \n",
    "                   = \\frac{1}{\\sqrt{N_x}}e^{\\I k_n x_m}, \\qquad\n",
    "  \\braket{x_m|x_n} = \\braket{k_m|k_n} = \\delta_{nm},\\\\\n",
    "  \\mat{1} = \\sum_{n}\\ket{x_n}\\bra{x_n} = \\sum_{m}\\ket{k_n}\\bra{k_n}.\n",
    "$$\n",
    "\n",
    "In this normalization $\\braket{x_m|k_n}$ is unitary.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Discrete Unitarity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fred Gittes points out that in discrete form, there is a special relationship between $\\d{t}$ and $\\d{x}$ such that $\\mat{M}^\\dagger\\cdot\\mat{M}$ is unitary:\n",
    "\n",
    "$$\n",
    "  [(A\\mat{M})^\\dagger \\cdot (A\\mat{M})]_{xy} = \n",
    "  \\sum_{q} (AM)_{qx}^*(AM)_{qy}\n",
    "  =\n",
    "  \\abs{A}^2\n",
    "  \\sum_{l=0}^{N_x-1}\n",
    "  \\exp\\left\\{\n",
    "    \\frac{m\\I}{2\\hbar\\d{t}}\n",
    "    [(x_l-y)^2 - (x-x_l)^2]\n",
    "    -\\frac{\\I}{\\hbar}\\frac{V(y)-V(x)}{2}\\d{t}\n",
    "  \\right\\}\\\\\n",
    "  =\n",
    "  \\frac{m}{2\\pi\\hbar\\d{t}}\n",
    "  \\exp\\left\\{\n",
    "    \\frac{m\\I}{2\\hbar\\d{t}}(y^2-x^2)\n",
    "    -\\frac{\\I}{\\hbar}\\frac{V(y)-V(x)}{2}\\d{t}\n",
    "  \\right\\}\n",
    "  \\sum_{l=0}^{N_x-1}\n",
    "  \\exp\\left\\{\n",
    "    \\frac{m\\I}{2\\hbar\\d{t}}[2x_l(x-y)]\n",
    "  \\right\\}\n",
    "  =\n",
    "  N_x\\delta_{xy}\n",
    "$$\n",
    "\n",
    "where the last relationship is true for the lattice iff:\n",
    "\n",
    "$$\n",
    "  \\frac{m}{\\hbar\\d{t}}\\d{x}^2 = \\frac{2\\pi n}{N_x}, \\qquad\n",
    "  \\d{t} = \\frac{m L\\d{x}}{2\\pi n\\hbar}.\n",
    "$$\n",
    "\n",
    "*To see this, note that with the special $\\d{t}$, the last sum is $N_x\\delta_{xy}$.  This leaves the first factor, whose exponent vanishes when $x=y$ leaving a factor of unity.*\n",
    "\n",
    "The discrete matrix $\\mat{M}$ is thus unitary with this special value of $\\d{t}$ if the normalization factor of $A = 1/N_x$ at this special value.  As we shall see below,\n",
    "\n",
    "$$\n",
    "  A = \\sqrt{\\frac{m}{2\\pi \\I\\hbar\\d{t}}}\n",
    "  \\\\\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the numerical basis, we can diagonalize $M_{xy}$.\n",
    "$$\n",
    "  \\tilde{M}_{mn} = \\sum_{mn}\\braket{k_m|x}M_{xy}\\braket{y|k_n} \n",
    "  = \n",
    "  \\frac{1}{N_x}\n",
    "  \\sum_{xy}\n",
    "  \\exp\\left\\{\n",
    "      -\\I (k_mx - k_ny)\n",
    "      +\n",
    "      \\frac{\\I}{\\hbar}\n",
    "      \\Bigl[\n",
    "        \\frac{m}{2}\\frac{(\\overbrace{x-y}^{q})^2}{\\d{t}}\n",
    "      \\Bigr]\n",
    "    \\right\\} = \\\\\n",
    "  = \n",
    "  \\frac{1}{N_x}\n",
    "  \\sum_{qy}\n",
    "  e^{\\I(k_n-k_m)y}\n",
    "  \\exp\\left\\{\n",
    "      -\\I k_m q\n",
    "      +\n",
    "      \\frac{\\I}{\\hbar}\n",
    "      \\left[\n",
    "        \\frac{m}{2}\\frac{q^2}{\\d{t}}\n",
    "      \\right]\n",
    "    \\right\\}\n",
    "  = \n",
    "  \\delta_{mn}\n",
    "  \\sum_{q}\n",
    "  \\exp\\left\\{\n",
    "      \\frac{\\I}{\\hbar}\n",
    "      \\left[\n",
    "        \\frac{m}{2}\\frac{q^2}{\\d{t}}\n",
    "      \\right]\n",
    "      - \\I q k_n\n",
    "    \\right\\}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Free Particle on a Periodic Lattice"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The last sum can be simplified if we can shift $q \\rightarrow q + l\\d{x}$ to complete the square which requires\n",
    "\n",
    "$$\n",
    "  l = \\frac{\\hbar k_n}{m}\\frac{\\d{t}}{\\d{x}}\n",
    "$$\n",
    "\n",
    "to be an integer.  This is the case if $\\d{t}$ is chosen in terms of an integer $j$:\n",
    "\n",
    "$$\n",
    "  \\d{t} = j\\frac{m L\\d{x}}{2\\pi \\hbar}, \\qquad\n",
    "  l = jn.\n",
    "$$\n",
    "\n",
    "If this is true, then we can do this shift and we have:\n",
    "\n",
    "$$\n",
    "  \\tilde{M}_{mn} = \\delta_{mn}\n",
    "  \\exp\\left\\{\n",
    "    \\frac{\\I m l^2 \\d{x}^2}{2\\hbar \\d{t}}\n",
    "  \\right\\}\n",
    "  \\sum_{q}\n",
    "  \\exp\\left\\{\n",
    "    \\frac{\\I m q^2}{2\\hbar \\d{t}}\n",
    "  \\right\\}\n",
    "  =\n",
    "  \\delta_{mn}\n",
    "  \\exp\\left\\{\n",
    "    \\frac{\\I}{\\hbar} \\frac{\\hbar^2 k_n^2}{2m} \\d{t}\n",
    "  \\right\\}\n",
    "  \\sum_{q}\n",
    "  \\exp\\left\\{\n",
    "    \\frac{\\I m q^2}{2\\hbar \\d{t}}\n",
    "  \\right\\}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Harmonic Oscillator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline --no-import-all\n",
    "\n",
    "Nx = 62\n",
    "L = 10.0\n",
    "w = m = h = 1.0\n",
    "dx = L/Nx\n",
    "x = np.arange(Nx)*dx - L/2\n",
    "\n",
    "n_t = 0\n",
    "Nt = 2*n_t+1\n",
    "dt = m*L*dx/2/np.pi/Nt/h\n",
    "T = Nt*dt\n",
    "\n",
    "def V(x):\n",
    "    return m*(w*x)**2/2.0\n",
    "\n",
    "x_ = x[:,None]\n",
    "y_ = x[None,:]\n",
    "\n",
    "A = np.sqrt(m/2/np.pi/h/T)\n",
    "AM = dx*np.exp(1j/h*(m*(x_-y_)**2/2.0/dt - dt*(V(x_) + V(y_))/2))*A\n",
    "I = AM.conj().T.dot(AM)\n",
    "assert np.allclose(np.eye(Nx), I)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0xa1a2b6fd0>]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAD8CAYAAABzTgP2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XlclOe99/HPD4ZVEESQHRFxRREjblFj4pKYxMQkTxK1idm30zRpn67JyTnP0/bpdtLTmrZpk5rVNotmt3Vp4pLFJWrAIKISUUBAEAFlUXbmev5gzEGDojJwzzC/9+s1L2bubb6TgL+5r/u+rkuMMSillFKneVkdQCmllGvRwqCUUuoMWhiUUkqdQQuDUkqpM2hhUEopdQYtDEoppc6ghUEppdQZul0YRGSliGQ5HoUikuVYPqnD8t0icvM59h8iIjtEJM9xLN/uZlJKKXXpxJkd3ETkd0CNMebnIhIINBtjWkUkGtgNxBhjWs/a5y3gPWPMChF5HthtjHnOaaGUUkpdFKcVBhERoAiYZYzJO2vdEGA7ENuxMDj2qQCiHAVkKvBTY8w153uv8PBwk5iY6JTcSinlKTIzMyuNMRFdbWdz4nvOAMo7FgURmQy8DAwGlpx9tgAMBKo7LC8BYrt6o8TERDIyMpyTWimlPISIHL6Q7S6oMIjIBiCqk1VPGWNWOZ4vBt7suNIYswNIEZFRwHIRWWeMaex46E6O2ekpjIg8BDwEkJCQcCGxlVJKXYILKgzGmDnnWy8iNuAWYMI59t8vIqeAMUDHr/qVQKiI2BxnDXFA6TmOsQxYBpCenq4j/ymlVA9x1u2qc4BcY0zJ6QWOu41sjueDgRFAYcedTPsFjo+BWx2L7gZWoZRSyjLOKgyLOKsZCZgO7Hbcvvo+8G1jTCWAiKwVkRjHdj8Bvi8iB2m/5vCSkzIppZS6BE69XbW3pKenG734rJRSF0dEMo0x6V1tpz2flVJKnUELg1JKqTM4sx+DUh6ttc1O0fF6DlWcIr/iJIG+3gyNCGLooCAGBfvR3p9TKdenhUGpS1BR10ROaQ37SmvZW1rDgfKTHK46RUtb59fsgv1sJA0KYmRkMCmx/UmJCWFUdDCBvvonqFyP/lYq1QVjDPmVp9iSV8mWg5XsLq7mWF3T1+sTwgIZERXMnFGRJA8KYmhEP5IigmhobuNQxUkOHjv59c+P9h1lZUYxAF4CSRFBpA8ewPRh4Vw+NJywfjqGpLKeFgalzmF7fhXvZpaw9WAlpTXtHfbjBgQwLTmclJj2b/2jY/oTEuDT6f4hAT5EhfgzLTn862XGGEprGtl7pIac0lpyjtSwJruMFV8UIwIpMf2ZlhzO3VMTiQkN6JXPqdTZ9HZVpTqxr7SWBX/eQoCPN9OSw5k+LJzpyeEMHtjP6e/V2mYn+0jN12ckuw6fIHlQEKu+Mw0/m7fT3095rgu9XVULg1JnaWxp46Y/b6XqVDMffu+KXm/e2ZRbzn2vZvDwzCSevHZUr7636tu0H4NSl+j36w+Qe7SOp29NtaTNf9bISBZPSmDZZ/nsyK/q9fdXSguDUh1sz6/ihc353DE5gatGDLIsx39cP4qEsEB+8PZu6hpbLMuhPJMWBqUcahtb+MFbuxkcFshT11vbhNPPz8bvb0+jtLqBn/9zn6VZlOfRwqCUw8/+sY+ymgZ+vzDNJfoXTBg8gG9fmczbmSV8tPeo1XGUB9HCoBTwr5wy3t1VwqNXJXNZwgCr43zt8dnDSInpz5Pv7aGiQ98JpXqSFgbl8b4sOsGP3s5mTGx/Hp89zOo4Z/C1efHMwjTqmlp58G8Z1Or1BtULtDAoj5Z5+ARLXtpJWJAvy5ak4+Pten8SwyKDeXbxePaW1rDkpZ3UNGhxUD3L9f4KlOolmYePc/fLOwkP8mXFQ1Ncuqfx1SlR/OWOCewrrWHJSzuoqdfioHqOFgblkTIKj3PXSzuJCPZjxUNTiQ5x3aJw2tzRkTx/5wRyy+q4U4uD6kFaGJRHMcbw2YEK7n55J5H9/Vnx0BSiQvytjnXBZo+K5Pkll/HV0TrueGk7ZTUNVkdSfZAOiaH6vPrmVrbkVfLxV8f4OLeCo7WNJEX0Y8WDUxjU332KQkcff3WMh/+eSXOrnTGx/Zk1YhBXjRzEuLhQvLx03gfVOR0rSSngF6v38bfth2lutRPkZ2PGsHCuGjmIeWOi6O/f+aio7qKg8hTrcsr4OPcYmYdPYDcwsJ8vv70tlVkjI62Op1yQFgbl8b46Wsc1z3zGvJQo7po6mPTEMHxtfbP1tLq+mU8PVPCHDXkAbPj+TD1zUN/QK4PoichKEclyPApFJMuxfFKH5btF5OZz7P+qiBR02DatO3mU6uiFzfkE+Hjzm/81lsuTw/tsUQAIDfRlQVos35s7nPzKU2zMPWZ1JOXGutXv3xiz8PRzEfkdUON4mQOkG2NaRSQa2C0i/zTGtHZymB8ZY97pTg6lznastpFVWUdYPCmB0EDPmRXtujFR/FdoAC9szmfuaG1OUpfGKV+hpH2W89uBNwGMMfUdioA/4H7tVcqtLf+8kFa74b5pQ6yO0qts3l7cOy2RnQXH2V1cbXUc5aacdW49Ayg3xuSdXiAik0VkL7AHeOQcZwsAvxSRbBFZKiJ+TsqjPFh9cyuvbS/imtFRJIY7f8Y1V7dwYjzBfjZe2JxvdRTlprosDCKyQURyOnks6LDZYhxnC6cZY3YYY1KAicCTItLZfYFPAiMd24QBPzlPjodEJENEMioqKi7goylP9XZGCTUNLTx4hWedLZwW7O/D4skJrMs5SvHxeqvjKDfUZWEwxswxxozp5LEKQERswC3AynPsvx84BYzpZF2ZadcEvAJMOk+OZcaYdGNMekRExIV9OuVx2uyGl7YUMD4hlAmDw6yOY5l7Lk9EgFe2FlodRbkhZzQlzQFyjTElpxeIyBBHwUBEBgMjgMKzd3RcmD59jeIm2i9aK3XJPtp7lKLj9Tw0I8nqKJaKCQ1gfmo0K78o0kH31EVzRmFYxFnNSMB02u9EygLeB75tjKkEEJG1IhLj2O51EdlD+3WIcOAXTsijPNgLm/NJCAvk6pQoq6NY7oEZSZxqbmPFziKroyg30+1pqowx93Sy7O/A38+x/XUdns/q7vsrdVrm4ePsKqrmpzeMxls7dzEmNoSpSQN5ZWsh904b0qf7cSjn0t8U1We88FkBIQE+3JYeb3UUl/HQFUkcrW1kdXap1VGUG9HCoPqEQxUn+XDfUe6ckkA/P+vna3YVM4dHMGxQEMs+y8cdh79R1tDCoPqEFz7Lx8fbi3su98xbVM/Fy0t4eOZQco/W8clXepu3ujBaGJTbK69t5L1dR7g9PY6IYO0jebYbx8UQE+LPc58csjqKchNaGJTbe3lLAa12Ow/NGGp1FJfka/Pi/hlJ7Cw8Tubh41bHUW5AC4NyazUNLby+o4jrU2NIGBhodRyXtWhiPKGBPjz3iQ6TobqmhUG5tde2H+ZkUyuPzPTsDm1d6edn4+6piWzYX05eeZ3VcZSL08Kg3FZjSxuvbC3giuERpMSEWB3H5d19eSL+Pl48/6meNajz08Kg3NY7mSVUnmzm32bqtYULEdbPl0UTE1iVdYTS6gar4ygXpoVBuaXWNjvLPstnXHwoU5I8d7C8i/XAjCEY4MXNBVZHUS5MC4NyS+ty2gfL+7eZQ2kfg1FdiLgBgSwYF8OKL4o4carZ6jjKRWlhUG7HGMNznxwiKaIfV+v0lRft4ZlDqW9uY/nnhVZHUS5KC4NyOx9/dYx9ZbU8MnMoXjpY3kUbERXMnFGRvLK1kJNN55pYUXkyLQzKrRhjeHbTQWJDA7h5fKzVcdzWd2Ylt/cB2X7Y6ijKBWlhUG7l8/wqdhVV88jMJHy89df3UqXFhzJjWDgvbC6gsaXN6jjKxehflnIrz246SESwnw6t7QSPXpVM5ckmVn5RbHUU5WK0MCi3savoBNsOVfHQjCT8fbytjuP2Jg8JI33wAP766SGaW+1Wx1EuRAuDcht/3nSQ0EAfvjU5weoofYKI8J1ZyZTWNPL+lyVd76A8hhYG5Rb2ltawMfcY908bohPxONHM4RGMjQ3huU8O0dqmZw2qnRYG5Rb+8vEhgv1s3HV5otVR+hQR4dGrkimsqmfNnjKr4ygXoYVBubyDx06yNqeMJVMHExLgY3WcPufq0ZEMjwzizx8fxG7X6T+VEwqDiKwUkSzHo1BEss5anyAiJ0Xkh+fYf4iI7BCRPMexfLubSfUtz31yCD+bF/dP12k7e4KXV/tZw4Hyk3y0r9zqOMoFdLswGGMWGmPSjDFpwLvAe2dtshRYd55D/Bew1BgzDDgB3N/dTKrvKKqq54OsI3xr0mAGBum0nT3l+rHRJA4M5NmP8zBGzxo8ndOakqR9JLPbgTc7LLsJyAf2nmefWcA7jkXLgZuclUm5v798chBvL+FhnYinR9m8vfj2VcnkHKnl46+OWR1HWcyZ1xhmAOXGmDwAEekH/AT42Xn2GQhUG2NOD9hSAug4BwqAI9UNvLurhEUT44ns7291nD7v5vGxxA0I4I8bD+pZg4e7oMIgIhtEJKeTx4IOmy2mw9kC7QVhqTHm5PkO3cmyTn8jReQhEckQkYyKiooLia3c3POfHALgEZ2Ip1f4eHvxb1cOJau4mi0HK62Ooyx0QTeEG2PmnG+9iNiAW4AJHRZPBm4VkaeBUMAuIo3GmGc7bFMJhIqIzXHWEAeUniPDMmAZQHp6un6d6eOO1jSy8otibp0QT0xogNVxPMatE+J4dtNB/rgxj+nJ4TrXhYdyVlPSHCDXGPN190ljzAxjTKIxJhF4BvjVWUUB036++jFwq2PR3cAqJ2VSbuyvnx2izRi+faWeLfQmP5s3j8wcyheFJ9ief9zqOMoizioMizizGem8RGStiMQ4Xv4E+L6IHKT9msNLTsqk3FRFXRNv7Cji5vGxxIcFWh3H4yycGE9EsB9/2pRndRRlEaeMLWCMuaeL9T896/V1HZ7nA5OckUP1DS9uzqelzc6jVyVbHcUj+ft48/AVSfxizX4yCo+Tnqhzansa7fmsXMrxU838ffthbhwXw5DwflbH8VjfmpxAWD9f/rjpoNVRlAW0MCiX8tKWfBpa2vjOLD1bsFKgr40HZyTx2YEKsoqrrY6jepkWBuUyaupbWL7tMNeNjSZ5ULDVcTzekqmDCQ304Vm91uBxtDAol/HS1gJONrXymJ4tuIQgPxv3TxvChv3HyDlSY3Uc1Yu0MCiXUNPQwitbC5iXEsXIqP5Wx1EOd09LJNjfpncoeRgtDMolLN9WSF1jK4/N1rMFV9Lf34f7pg3hw73l7C+rtTqO6iVaGJTl6hpbeGlLAXNGRZISE2J1HHWW+6YNIcjPxrN6h5LH0MKgLPe3zw9T09DC43q24JJCAn245/JE1uaUcaC8zuo4qhdoYVCWOtXUyoub87lqRASpcaFWx1HncP/0IQT4eOtZg4fQwqAs9dr2w5yob+Gx2cOsjqLOY0A/X+6amsg/s0s5eOx8AyarvkALg7JMQ3Mbyz7LZ8awcC5LGGB1HNWFB2YMwd/mzV8+1rOGvk4Lg7LM6zsOU3Wqmcf1bMEthAf5ceeUBD7IOkJh5Smr46gepIVBWaLNbnh5SwGTh4QxUQdpcxsPXpGEzduLV7YWWB1F9SAtDMoSWw9WUlrTyJKpg62Ooi7CoGB/5qVE8UFWKY0tbVbHUT1EC4OyxFsZxYQG+jB3dKTVUdRFuj09npqGFtbvK7c6iuohWhhUr6uub+ajveXclBaLn83b6jjqIl0+dCCxoQG8lVFsdRTVQ7QwqF63KquU5jY7t6XHWR1FXQIvL+HWCXFsOVjJkeoGq+OoHqCFQfW6tzKKSYnpr8NfuLFbJ8RhDLybWdL1xsrtaGFQvSrnSA17S2u5PT3e6iiqG+LDApmWPJC3M4ux243VcZSTaWFQveqdzBJ8vb1YkBZjdRTVTbenx1N8vIHtBVVWR1FOpoVB9ZrGljbe//IIV6dEEhroa3Uc1U3XpEQR7G/j7QxtTuprtDCoXrNhfzk1DS3ajNRH+Pt4syAthrV7yqhtbLE6jnKibhUGEVkpIlmOR6GIZJ21PkFETorID8+x/6siUtDhGGndyaNc21sZJcSE+DMtOdzqKMpJbk+Pp6nVzj93l1odRTlRtwqDMWahMSbNGJMGvAu8d9YmS4F1XRzmR6ePYYzJ6mJb5aZKqxvYnFfBrenxeHuJ1XGUk4yNDWFkVDBvaXNSn+KUpiQREeB24M0Oy24C8oG9zngP5d7ezSzBGLhtgvZd6EtEhNvS49ldXM1XR3USn77CWdcYZgDlxpg8ABHpB/wE+NkF7PtLEckWkaUi4neujUTkIRHJEJGMiooK56RWvcJuN7yVWczUpIHEhwVaHUc52U1pMfh4Cyu/0J7QfUWXhUFENohITiePBR02W0yHswXaC8JSY0xXM3o8CYwEJgJhtBeTThljlhlj0o0x6REREV3FVi5k26Eqio83sGiSXnTuiwYG+XH16Cje/7KEplYdWK8vsHW1gTFmzvnWi4gNuAWY0GHxZOBWEXkaCAXsItJojHn2rGOXOZ42icgrQKcXqZV7e/OLIkIDfbgmJcrqKKqHLJwYz5o9ZXy4t5wbx2kfFXfnjKakOUCuMebrq0/GmBnGmERjTCLwDPCrs4sCgIhEO34KcBOQ44Q8yoUcP9XMR3uPcvP4WPx9dMC8vmp6cjhxAwJYsbPI6ijKCZxRGBZxZjPSeYnIWhE5/ZXidRHZA+wBwoFfOCGPciHv7Sqhpc2waGKC1VFUD/LyEhamx7PtUBWHq3R2N3fXZVNSV4wx93Sx/qdnvb6uw/NZ3X1/5bqMMby5s4jxCaGMiAq2Oo7qYbemx7F0wwFWflHMj+eNtDqO6gbt+ax6TObhExyqOMViPVvwCNEhAVw1YhBvZ5bQ0ma3Oo7qBi0Mqse8ubOYID8b88dFWx1F9ZJFkxKoqGtiU+4xq6OobtDCoHpETUMLa/aUcmNaDIG+3W6xVG7iqhERDAr20z4Nbk4Lg+oR/9hdSmOLnUUTte+CJ7F5e3FbehyffHWMshqd3c1daWFQPWLFziJGR/dnbKzO0uZpFqYnYDfw1hc6fpK70sKgnG5PSfssbYsnxdPeRUV5koSBgUxPDuetjGLadHY3t6SFQTndii+K8Pfx4sa0WKujKIssmhTPkeoGthystDqKugRaGJRTNbfaWZ1dxryUKEICfKyOoywyd3QkIQE+fPDlEaujqEughUE51daDldQ0tHCDjpfj0fxs3lyTEsn6feU0tujAeu5GC4NyqtXZZfT3tzFjmI6A6+nmp8ZwsqmVTw/oMPnuRguDcpqm1jY+2neUa1Ki8LXpr5anu3zoQAYE+rA6u6zrjZVL0b9e5TSfHaikrrGV61O1p7Nq79Mwb0w0G/eX09CszUnuRAuDcprV2aUMCPRhWnK41VGUi7ghNZr65jYdIsPNaGFQTtHY0saGfeXMGxOFj7f+Wql2k5MGEh7kx5o9pVZHURdB/4KVU3yce4xTzW3MT9W7kdT/8PYSrhsbxabcY5xqarU6jrpAWhiUU6zeU8bAfr5MHhJmdRTlYq4fG01ji50N+8utjqIukBYG1W31za1s2n+Ma8dGYdNmJHWWiYlhRPb3Y43eneQ29K9YddvG/cdoaNFmJNU5Ly/hurHRfHKggrrGFqvjqAughUF12+rsUgYF+zExUZuRVOfmp8bQ3Gpn/T5tTnIHWhhUt9Q1tvDxVxVcNzYaby8dSVV1bnx8KDEh/trZzU1oYVDdsnH/MZpb7czXTm3qPLy8hOtTo9mcV0FNvTYnubpuFQYRWSkiWY5HoYhkOZYnikhDh3XPn2P/MBFZLyJ5jp8DupNH9b73vjxCdIg/lyXo/zp1fvNTY2hpM/wzW/s0uLpuFQZjzEJjTJoxJg14F3ivw+pDp9cZYx45xyGeADYaY4YBGx2vlZv4/FAVnx2o4M4pg/HSZiTVhdS4EMbFh/KnTXnUN2ufBlfmlKYkaZ+m63bgzYvcdQGw3PF8OXCTM/Konme3G365dh8xIf7cP32I1XGUGxAR/vP6UZTXNvHCZwVWx1Hn4axrDDOAcmNMXodlQ0TkSxH5VERmnGO/SGNMGYDj56BzvYGIPCQiGSKSUVGhw/ha7YOsI+QcqeVH80bg7+NtdRzlJtITw7h2TBR//ewQx2obrY6jzqHLwiAiG0Qkp5PHgg6bLebMs4UyIMEYMx74PvCGiPTvTlBjzDJjTLoxJj0iQsf6t1JDcxu//fArUuNCWDBOp+9UF+cn80bS0mbn9+sPWB1FnYOtqw2MMXPOt15EbMAtwIQO+zQBTY7nmSJyCBgOZJy1e7mIRBtjykQkGtAhGN3Ay1sLKKtpZOnCNL22oC5aYng/lkxJ5NVtBdwzLZGRUd36zqh6gDOakuYAucaYktMLRCRCRLwdz5OAYUB+J/v+A7jb8fxuYJUT8qgeVFHXxF8+Psjc0ZFMSRpodRzlph6fnUywvw+/WptrdRTVCWcUhkV886LzFUC2iOwG3gEeMcYcBxCRF0Uk3bHdb4C5IpIHzHW8Vi7smQ0HaGq18+S1I62OotxYaKAvj81K5rMDFTr1pwvqsimpK8aYezpZ9i7tt692tv0DHZ5XAbO7m0H1jrzyOt7cWcRdUxNJigiyOo5yc0umDuZvnx/mV2v2Mz05XHvOuxDt+awu2NMffkU/PxuPzx5mdRTVB/jZvHni2pF8VV7He7tKut5B9RotDOqCFFSeYv2+cu69PJGwfr5Wx1F9xLVjohgZFcxLWwowxlgdRzloYVAXZPm2Qny8hTunDrY6iupDRIT7pg0h92gdn+dXWR1HOWhhUF2qbWzh7YxibkiNYVCwv9VxVB9zY1oMYf18eXlLodVRlIMWBtWlt74o5lRzG/dO06EvlPP5+3hzx+QENuaWc7jqlNVxFFoYVBfa7IblnxcyMXEAY+NCrI6j+qg7pwzG5iW8uq3Q6igKLQyqCxv2l1N8vIH79GxB9aDI/v5cPzaatzNKdPpPF6CFQZ3XK1sLiA0NYO7oSKujqD7uvulDONnUytsZeuuq1bQwqHPaW1rD9vzj3H35YGze+quielZqXCgTBg9g+eeFtNn11lUr6V+7OqdXtxYS4OPNwvQEq6MoD3HvtEQOV9WzKVfH07SSFgbVqcqTTazKKuXWCXGEBPpYHUd5iHkpUcSE+PPKVp3Ix0paGFSn3thRRHObnXumJVodRXkQm7cXS6Ymsu1QFfvLaq2O47G0MKhvMMbwVkYxM4aFM1QHy1O9bPGkeHy9vXgro9jqKB5LC4P6hqziakpONLAgTWdnU70vNNCXmSMiWLunDLtehLaEFgb1Dauzy/D19uLqFL1FVVljfmo05bVNfFF43OooHkkLgzqD3W5Yk13GFcMj6O+vF52VNeaMisTfx4vV2WVWR/FIWhjUGTKLTnC0tpEbxkVbHUV5sH5+NmaNHMS6nDJa2+xWx/E4WhjUGVbvLsXP5sXsUdqMpKw1PzWGypPN7CjQ5qTepoVBfa3Nblibc5RZIwcR5NftWV+V6parRgwi0Neb1dmlVkfxOFoY1Nd2FFRRUdfE/NQYq6MoRYCvN3NGRbIu5ygt2pzUq7QwqK+tzi4j0NebWSMHWR1FKaD97qTq+ha2Hqy0OopH6VZhEJGVIpLleBSKSJZjeaKINHRY9/w59v+piBzpsN113cmjLl1rm51/5Rxl9qhIAny9rY6jFAAzR0QQ7GfTu5N6Wbcako0xC08/F5HfATUdVh8yxqRdwGGWGmP+uzs5VPdtO1TF8VPNzE/Vu5GU6/CzeTM3JZIP9x7llzePwc+mX1p6g1OakkREgNuBN51xPNX7VmeXEuxnY+bwCKujKHWGG1JjqGtsZfMBbU7qLc66xjADKDfG5HVYNkREvhSRT0Vkxnn2/Y6IZIvIyyIywEl51EVobm1vRpo7OhJ/H/1GplzLtORwQgJ89O6kXtRlYRCRDSKS08ljQYfNFnPm2UIZkGCMGQ98H3hDRPp3cvjngKFAmmOf350nx0MikiEiGRUVFRfw0dSF2nKwgtrGVuZrpzblgnxtXsxLiWL9vnIaW9qsjuMRuiwMxpg5xpgxnTxWAYiIDbgFWNlhnyZjTJXjeSZwCBjeybHLjTFtxhg78AIw6Tw5lhlj0o0x6RER2tzhTP/IKiUkwIfpyfrfVbmm+eOiOdXcphP49BJnNCXNAXKNMV9P1CoiESLi7XieBAwD8s/eUUQ6fkW9GchxQh51EY7VNrJmTxkL0mLwtendy8o1TU0aSGxoAK9uLbQ6ikdwxr8Ei/jmRecrgGwR2Q28AzxijDkOICIviki6Y7unRWSPiGQDVwH/2wl51EVY/nkhrXbDfdOGWB1FqXOyeXtx77REdhYeZ3dxtdVx+jwxxv3GO09PTzcZGRlWx3B79c2tTP31JqYkhfHXJeld76CUheoaW7j815uYOSKCZ791mdVx3JKIZBpjuvxj17YDD/Z2Rgk1DS08OCPJ6ihKdSnY34fFkxNYl3OUkhP1Vsfp07QweKg2u+GlLQWMTwhlwmC9S1i5h3suT0SAV/RaQ4/SwuCh1u87StHxeh6ckUR7/0SlXF9MaADXp0azYmcRNQ0tVsfps7QweKgXNhcQHxbANSlRVkdR6qI8OCOJU81trNhZZHWUPksLgwfKPHyCzMMnuH/aELy99GxBuZcxsSFMTRrIq9sKdTjuHqKFwQO9uDmf/v42bkuPtzqKUpfkwSuGUFbTyBoddbVHaGHwMEVV9Xy49yh3TBlMP52lTbmpK4cPYmhEP17YnI873nLv6rQweJgXNufj7SXcc3mi1VGUumReXsIDM5LYW1rL1oNVVsfpc7QweJD9ZbW8sbOI29Ljiezvb3Ucpbrl5vGxxIYG8PPVe/Vag5NpYfAQdrvhqff3EBLgw4+uHmF1HKW6zd/Hm5/dmMKB8pO8tKXA6jh9ihYGD/FWRjG7iqp58tqRDOjna3UcpZxizuhI5o6O5A8b8rQ3tBNpYfAAVSeb+PW6XCYNCePWCXFWx1HKqX56Y0r7z3/2sKgVAAAP5ElEQVTsszhJ36GFwQP8el0up5pa+cVNY7SXs+pzYkMD+N6cYWzYX876feVWx+kTtDD0cTvyq3gns4QHZiQxPDLY6jhK9Yj7pg9hRGQwP/3HXuqbW62O4/a0MPRhza12/uODHGJDA3h8drLVcZTqMT7eXvzi5jEcqW7gDxvzut5BnZcWhj7s5a0F5B07yc9uTCHQVzuzqb5tYmIYt6fH8dLmAg6U11kdx61pYeijKuqa+NPGPOaMGsSc0ZFWx1GqVzxx7Sj6+dn45Zr9Vkdxa1oY+qg/bDxAY6udJ68bZXUUpXpNWD9fHpuVzKcHKticV2F1HLelhaEPOnjsJG/uLOaOyQkMjQiyOo5SvWrJ1MHEhwXwyzX7abPrOEqXQgtDH/SbdbkE+Hjz3dnDrI6iVK/zs3nz42tGknu0jvd2lVgdxy1pYehjtudXsWF/Of925VAGBvlZHUcpS8xPjWZcfCj//dFXNDS3WR3H7XSrMIjIShHJcjwKRSSrw7pUEflcRPaKyB4R+caobSISJiLrRSTP8VMnH+4Gu93wq7X7iQ7x5/7pQ6yOo5RlRISnrhtFeW0TL23JtzqO2+lWYTDGLDTGpBlj0oB3gfcARMQGvAY8YoxJAa4EOpug9QlgozFmGLDR8brHlNc29umekf/MLiW7pIYfXj0Cfx9vq+MoZalJQ8K4enQkz31yiIq6JqvjuBWnNCVJ+zgLtwNvOhZdDWQbY3YDGGOqjDGdnc8tAJY7ni8HbnJGnnP5zbpcHn1jF/vLanvybSzR2NLG0//6itHR/bl5fKzVcZRyCU9cO5KmVjt/2HjA6ihuxVnXGGYA5caY010OhwNGRD4UkV0i8uNz7BdpjCkDcPwc5KQ8nXrq+lGEBPjw2Jtf9rl2x79/fpgj1Q08df0ovHQeZ6UASIoI4luTE3hzZzEHj520Oo7b6LIwiMgGEcnp5LGgw2aL+Z+zBQAbMB24w/HzZhGZ3Z2gIvKQiGSISEZFxaXdnxwe5MfS29M4VHGSn6/uOyMxtrTZeXFLPtOTw5mWHG51HKVcyndnD8PmJTpnw0XosjAYY+YYY8Z08lgFX19PuAVY2WG3EuBTY0ylMaYeWAtc1snhy0Uk2nGcaODYeXIsM8akG2PSIyIiLvwTnmX6sHAemTmUN3cW9ZmJxDfuL6e8tom7dbpOpb5hYJAfN4yLYVXWEeoaO7vUqc7mjKakOUCuMabjDcMfAqkiEugoHDOBzr6i/wO42/H8bmCVE/J06ftzhzMuPpQn3svuE5N7vL6jiJgQf2aN7NGWOKXc1p1TBlPf3MYHXx6xOopbcEZhWMSZzUgYY04Avwe+ALKAXcaYNQAi8qKIpDs2/Q0wV0TygLmO1z3Ox9uLPy0ajzHw3RVZtLrxfLEFlafYnFfJokkJeOu1BaU6NS4uhDGx/XltexHGaG/ornS7MBhj7jHGPN/J8teMMSmOZqcfd1j+gDEmw/G8yhgz2xgzzPHzeHfzXKiEgYH88uYxZB4+4dbD9L6x4zA2L2HRxHiroyjlskSEOyYP5qvyOjIPn7A6ziVpam3j6X/l9kpzmEf3fF6QFsttE+L406aD/CvH/a43NLa08XZmCVenRDKo/zf6DyqlOliQFkOwn43Xth+2OspFM8bw5Ht7+Msnh9iR3/Pfnz26MAD8v5vGkBYfyv9euZucIzVWx7koa/eUUV3fwp2TB1sdRSmXF+hr45bLYlm75yhVJ92rw9vzn+bz3q4jfHf2sF4ZRt/jC4O/jzfL7prAgEAfHlieQXlto9WRLthr2w+TFN6PqUMHWh1FKbdwx5TBNLfZeSfTfQbX+3DvUZ7+MJf5qdF8b07vDIzp8YUBYFCwPy/ePZHaxhYe/FuGW3R+21day66iar41OYH2judKqa4Mjwxm0pAw3thZhN0NhuTeW1rD91ZkkRobwn/fNq7X/ta1MDiMjunPHxaNZ8+RGn74zm6Xv3PhtR2H8bN5ceuEOKujKOVW7picwOGqejYfrLQ6ynkdq2vkgeUZhAb68MJd6b06/pkWhg7mjo7kJ/NGsia7jKUbXPdOpZNNraz68gjzU2MIDfS1Oo5SbmXemCgG9vPldRe+CN3Q3MZDf8ukur6FF+5K7/WbS7QwnOXhK5K4bUIcf9yYx4ubXXO43ve/PMKp5jbunJJgdRSl3I6fzZvbJ8azYX85ZTUNVsf5hsaWNh78Wwa7S6p5ZlEaY2JDej2DFoaziAi/vmUs14+N5hdr9vPKVtcaX6XNbnhlSwFjYvuTFh9qdRyl3NK3JrV/qXpla6G1Qc7S1NrGw3/PZOuhSn576ziuSYmyJIcWhk7YvL14ZlEa16RE8rN/7uPvLnTKuTq7lPzKUzx6ZbJedFbqEsWHBbIgLZa/f37YZW5dbW618+3XdvHpgQp+ffNYS68famE4Bx9vL/60+DJmjxzEf36Qw4qdRVZHwm43/GnTQYZHBln2TUKpvuLRq5JpbG1ziVFXW9rsPPbmLjbmHuP/3TSGRZOsbSbWwnAevjYv/nLnZcwcHsGT7+/hb58XWjqu0rqcoxw8dpLHZg3TOReU6qbkQUFcPzaa5dsKqa5vtixHTX0L31uRxYd7y/m/N4xmyRTrO6xqYeiCn82bvy6ZwPTkcP7Pqr3M/O0n/PXTQ73+i9R+tpBHUkQ/rhsb3avvrVRf9disYZxqbuNlC84a8srr+Pf39zDl1xtZs6eMp64bxb3TXGOudpvVAdyBv483r947ifX7ynl1WwG/XpfL0g0HuOWyOOanRjMg0JdgfxvB/j4E+dl6ZJTT9fvLyT1ax9KF43QUVaWcZERUMNeOieKVrYXcPyOJkAAfp79Hc6udusYW6hpbqWtspeREPa/vKGLLwUp8bV7clBbDPZcPYXRMf6e/96XSwnCBvL2EeWOimDcmin2ltSzfVsi7mSW8seOb1x6GDQriB1cP55qUKKdcIDbG8MeNeSQODOSG1JhuH08p9T++MyuZdTlHWb6tkMdnO2fIiVNNrby4uYBXthVQXf/N0VCj+vvzo2tGsHhSAmH9XK8vkhaGSzA6pj//dWsqT1w7kpzSGuoaWznZ2EptYwu1ja2s3VPGI6/tYnxCKE9eO4pJQ8K69X6bco+xt7SW396ais1bW/+UcqaUmBDmjIrkpS0F3DstkWD/Sz9raGmzs/KLYp7ZkEflySbmjo5kXFwIQX7tLQrB/jZCA30ZnxCKjwv/LYurD/3QmfT0dJORkWF1jHNqbbPz7q4Sfr/+AOW1TcwZNYgnrh1J8qDgiz6WMYab/rKN46ea2PSDK136l0kpd5VdUs2Nz27lR9eM4NGrki96f2MMH+4t5+l/5ZJfeYqJiQN44tpRTBg8oAfSXjoRyTTGpHe1nf4r0wNs3l4snJjAJz+8ih/PG8GO/ONc+4fNPPfJIdoucuCuj/aVs7u4mm9fmaxFQakekhoXypUjInhxc/5F92uoOtnEI69l8shrmXh7CS/elc5bD091uaJwMfSMoRdUnWziPz7IYV3OUdIHD+B3t49j8MB+nW5b09DC54cq2ZzX/ig6Xk/cgAA2/eBKfG1aGJTqKbuKTvC/ntuGAOPiQ5kxLIIZw8JJiz93s8/6feU8+V42tQ2t/ODq4dw/fYhLN/de6BmDFoZeYoxhVVYp/7kqhza74anrR7F4YgL5lafIKq4mq/gEWcXV7CutxW6gn683U4cOZMawCK4dE6UztCnVC/aU1LB+31E+y6sku6Qau4EgPxtp8aH/80gIxc/mxc//uY+3M0sYHd2fpQvTGBF18U3FvU0Lg4sqq2ngR29ns+VgJf4+XjS2tHeYC/KzkRoXQvrgAUwfFuHyF6eU6utq6lvYdqiSrYcq+bKomtyjdV83Bfv7eLUPYXFlMo/PHuY2Z/NaGFyY3W54K6OYvaW1jI0LIS0+lKERQdo/QSkX1tDcRk5pDVlF1eRXnuK29DguS3Cv6wgXWhi6dbuqiKwERjhehgLVxpg0x7pU4K9Af8AOTDTGNJ61/0+BB4EKx6J/N8as7U4md+DlJZaPhaKUujgBvt5MTAxjYmL3bj93B90qDMaYhaefi8jvgBrHcxvwGrDEGLNbRAYC3+zl0W6pMea/u5NDKaWU8zilg5u0d++9HZjlWHQ1kG2M2Q1gjKlyxvsopZTqec66YjIDKDfGnJ4PczhgRORDEdklIj8+z77fEZFsEXlZRNyrwU4ppfqgLguDiGwQkZxOHgs6bLYYeLPDaxswHbjD8fNmEZndyeGfA4YCaUAZ8Lvz5HhIRDJEJKOiouJcmymllOqmLpuSjDFzzrfecT3hFmBCh8UlwKfGmErHNmuBy4CNZx27vMNxXgBWnyfHMmAZtN+V1FVupZRSl8YZTUlzgFxjTEmHZR8CqSIS6CgcM4F9Z+8oIh0nFrgZyHFCHqWUUt3gjIvPizizGQljzAkR+T3wBWCAtcaYNQAi8iLwvDEmA3haRNIc2xQCDzshj1JKqW7QDm5KKeUh+nTPZxGpAA5bneMShAOVVoewgH5uz6Kf23UNNsZEdLWRWxYGdyUiGRdSrfsa/dyeRT+3+3OPkZ+UUkr1Gi0MSimlzqCFoXctszqARfRzexb93G5OrzEopZQ6g54xKKWUOoMWBouIyA9FxIhIuNVZeoOI/FZEch0DJr4vIqFWZ+opIjJPRL4SkYMi8oTVeXqDiMSLyMcisl9E9orId63O1JtExFtEvhSRcw7r4060MFhAROKBuUCR1Vl60XpgjDEmFTgAPGlxnh4hIt7An4FrgdHAYhEZbW2qXtEK/MAYMwqYAjzqIZ/7tO8C+60O4SxaGKyxFPgx7UOBeARjzEfGmFbHy+1AnJV5etAk4KAxJt8Y0wysABZ0sY/bM8aUGWN2OZ7X0f6PZKy1qXqHiMQB1wMvWp3FWbQw9DIRuRE4cnoSIw91H7DO6hA9JBYo7vC6BA/5B/I0EUkExgM7rE3Sa56h/Yue3eogzuKUGdzUmURkAxDVyaqngH+nfYa7Pud8n9sYs8qxzVO0Nzu83pvZepF0ssxjzgxFJAh4F/ieMabW6jw9TUTmA8eMMZkicqXVeZxFC0MPONccFiIyFhgC7G6fDZU4YJeITDLGHO3FiD3iAubuuBuYD8w2ffc+6RIgvsPrOKDUoiy9SkR8aC8Krxtj3rM6Ty+ZBtwoItcB/kB/EXnNGHOnxbm6RfsxWEhECoH00xMa9WUiMg/4PTDTGNNnp+BzzD9yAJgNHKF96PlvGWP2WhqshznmfV8OHDfGfM/qPFZwnDH80Bgz3+os3aXXGFRveRYIBtaLSJaIPG91oJ7guMD+Hdonq9oPvNXXi4LDNGAJMMvx/zfL8S1auSE9Y1BKKXUGPWNQSil1Bi0MSimlzqCFQSml1Bm0MCillDqDFgallFJn0MKglFLqDFoYlFJKnUELg1JKqTP8f8x3gR5uWWuZAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import scipy as sp\n",
    "import scipy.linalg\n",
    "H = 1j*h*sp.linalg.logm(AM)/dt\n",
    "assert np.allclose(H, H.T.conj())\n",
    "assert np.allclose(H.imag, 0)\n",
    "H = H.real\n",
    "plt.plot(x, np.diag(H))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0xa1a38f390>]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAD8CAYAAABzTgP2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8VNeVJ/DfqU2lfd93ITYhCSEEAgOG2DjGgMF2Yqed2HFip53JTGeZdKbbiZNMJ+lOd/Z8Oj2dtNt4yeYlsbGDNwx2jMCYRUgICYQkkIR2qbSUdpVqufOHqjDCElpque9Vne/nw8eo9Kx3ZKt03r3n3nNJCAHGGGPMRSM7AMYYY8rCiYExxtg0nBgYY4xNw4mBMcbYNJwYGGOMTcOJgTHG2DScGBhjjE3DiYExxtg0nBgYY4xNo5MdwGLExcWJrKws2WEwxpiqnDlzplcIET/XdapMDFlZWSgvL5cdBmOMqQoRXZnPdTyVxBhjbBpODIwxxqbhxMAYY2waTgyMMcam4cTAGGNsGk4MjDHGpuHEwBhjbBpODIx5wXt1PWg0jcgOg7FF4cTAmIfZ7A78j9+fwc8O1csOhbFF4cTAmIc19IxgwupAddug7FAYWxRODIx5mCshtPSPwTw2KTkaxhaOEwNjHlbd/uFIoaZ9SGIkjC0OJwbGPOxc+yDykiOcfzdLjoaxhePEwJgHTdocqO0cwpalcciKDeE6A1OlgEoMR+pNeKLssuwwmB+r7x7GpM2B/NRIFKRF4RwnBqZCAZUYyupN+Nnb9bDZHbJDYX7KVV8oTItEYWok2s3j6BuxSI6KsYUJqMRQkBoJi82Bhh7eeMS841zbICKMOmTEhKAgLRLA9GI0Y2oQUIkhP3XqjVrDb1TmJdXtZhSmRYGIsColAkTgOgNTnYBKDDlxoQg1aDkxMK+YsNpR1zV8daQQbtQjOy4U5/jnjamMW4mBiO4lovNE5CCikhtc9xQR9RBRzXWvxxDRISJqcP4z2p145qLREPJSInhoz7yirmsYVrtAoXNkCgCFqZE8YmCq4+6IoQbAPQDK5rjuGQA7Znj9MQDvCCGWAnjH+bFX5adG4kLnEBegmce5RgauEcPU36PQNTSBnqEJWWExtmBuJQYhRK0Qom4e15UB6J/hU3sBPOv8+7MA7nInnvkoSI3EhNWBxt5Rb9+KBZjqNjNiQg1IjQq++lohF6CZCsmuMSQKIToBwPnPhNkuJKJHiaiciMpNJtOib1jgHObz8J552rm2QRSkRoKIrr6WlxwBDYH3MzBVmTMxENFhIqqZ4c9eXwToIoR4QghRIoQoiY+PX/TXyYkPQ7Bey09wzKPGJ+1o6Bm5OkJwCQ3SITchjH/emKro5rpACLHdi/fvJqJkIUQnESUD6PHivQAAWmcBmlcmMU+60DkEu0NcHZFeqyA1CkfqTRBCTBtNMKZUsqeS/gLgIeffHwLwqi9uWuAsQNsdwhe3YwGgum2qWV5hWtRHPleYFoneEQu6uADNVMLd5ap3E1EbgI0AXieig87XU4jojWuuew7ABwCWE1EbET3i/NS/AbiNiBoA3Ob82OvyUyMxNmlHUy/vgGaeca59EPHhQUiMCPrI51yrlLjOwNRizqmkGxFC7Aewf4bXOwDsvObj+2f59/sA3OpODItxtQDdPojchHBf3575oeoZCs8ueckR0GoI1W2DuH1VkoToGFsY2VNJUiyJD4VRr0F1Gx+iwtw3arHhkmlkxvoCABj1WixLDOcd0Ew1AjIx6LQarEzmAjTzjPMdQxACH1mRdK2pHdBmCMF1LaZ8AZkYgKnppPMdg3BwAZq56Zyz8DzbiAGYqjMMjFnRNjDuq7AYW7SATQz5qZEYnbSjqY93QDP3VLcPIinCiIQI46zX8A5opiYBmxgKuAU385DqtsGrLd1nszwpHHotoaqNz4BmyhewiSE3IQwGnYZbYzC39I9OorF3FGsyPrp/4VpBOi3yUiJR2cKJgSlfwCYGvasA3cGJgS1exZUBAEBJ5twd40syo1HVasakjTv7MmUL2MQAAAWpETjfPsQFaLZoZ1oGoNPQjDuer7c2MxoWmwMXOnmZNFO2AE8MkRi22HClf0x2KEylzjQPYFVqJIIN2jmvXescVZQ3z9SBnjHlCOjEkJ/KK0XY4k3aHKhqM2NtxvwOHkyMMCItOhgVLQNejowx9wR0YliaEA6DVoPznBjYIlzoHILF5kBJ1vxPpC3JjEZ58wBvdGOKFtCJwaDTYEVyOI8Y2KK4poTWzqPw7LI2Mxo9wxbe6MYWbMRiw8sVbegdsXj9XgGdGICp6aTqdt4BzRauomUAadHBSLzBxrbrFTuTyJkrPJ3EFuZcqxlff7HKJ3uvAj4xFKVFYXjCxmdAswURQqC8eWBBowUAWJEUgVCDlhMDW7DK1qk9MEXpc6+AcxcnBufGpKpW3njE5q9tYBw9w5Z57V+4llZDWJMRjXJODGyBqlrNyI4LRVSIwev3CvjEsCQ+DGFBOpzlxMAWwPXEX7zAxABM1RnquoYwPGH1dFjMTwkhcLbVjNU36ODrSQGfGLQaQkFqJPewYQty5soAQg1arEiKWPC/uzYzGg4Bfhhh89Y5OIGeYYtPppEATgwApqaTajuHMGG1yw6FqUT5lQGsyYiGVvPRE9vmsiYjCkRcgGbz55rqLprnnhl3cWLAVDHHahc438GtCtjchiesqOsaWnDh2SXcqMfyxHBODGzezraaYdBqsDLZN0cRc2LAh1V+LkCz+TjbaoZDLGz/wvVKsqJR2WKGnZdJs3mobDVjZUoEgnRzt17xBE4MmGpVkBxp5DlfNi9nrgyACHO22r6RtZnRGLHYUNc17MHImD+y2R2obhvEGh/VFwBODFcVpUdxYmDzcubKAJYnhiPcqF/01yjJjJn6Wtw3ic2hoWcE41Y7Vqf7ZkUSwInhqtXpUWjpH0P/6KTsUJiC2R0ClS3mBfVHmkladDDiw4NwhjutsjmcvbqxzTeFZ4ATw1VcZ2DzUd89jBGLza36AgAQEUoyo3nEwOZU1WpGZLAeWbEhPrsnJwangtRIaOjDbeeMzaT86oltMW5/rbWZ0WjtH0fP0ITbX4v5r7OtZqxOjwLRwpdGLxYnBqfQIB2WJYbziIHdUMWVAcSHByEtOtjtr8UN9dhcRi021HcP+2xjmwsnhmsUpUehqs3MvfLZrMqv9GNtRrRHnt7yUyIRpNPgdDMnBjaz6vZBOAR8uiIJ4MQwTVF6FMxjVjT38VGf7KPazeNo7R/H+mz3p5GAqfNAitKjcKq5zyNfj/kfV+G50Ec9klw4MVxjNReg2Q2cbJz6Bb4hJ9ZjX3NDTiwudAxhcJwb6rGPqmo1IyMmBLFhQT69LyeGayxLDEeIQcv7GdiMTjb2IzJYjxVJnmtLUJoTA4f48DQ4xq51ttXs8/oC4GZiIKJ7ieg8ETmIqOQG1z1FRD1EVHPd6z8hootEdI6I9hOR7/8LXMPVaZVXJrGZnGzqw7qsGGgW0ThvNsUZ0TBoNTjZxImBTdc9NIHOwYmrMxm+5O6IoQbAPQDK5rjuGQA7Znj9EIB8IUQhgHoA33QzHrcVpUehtmMIFht3WmUf6hqcQHPfGDbkeKa+4GLUa1GUHnV1mooxl7M+PLHtem4lBiFErRCibh7XlQH4yCOREOJtIYTN+eEJAGnuxOMJRelRmLQ7UNvJPWzYh042eb6+4FKaE4Pq9kE+uIdNc7bVDJ2GsCpl4Wd+uEtJNYaHAbw52yeJ6FEiKieicpPJ5LUgXEd9nuUdqewaJxr7EW7UYWWy59+kpdmxU3UG3s/ArnG2xYyVyREw6n3TUfVacyYGIjpMRDUz/NnrqSCI6HEANgB/mO0aIcQTQogSIURJfHy8p279EUkRRiSEB6GqbdBr92Dqc7KxD+uzYhZ1MM9cijOjoNcSTjZynYFNsTsEqtsHpUwjAYBurguEENu9GQARPQRgN4BbhQJ2lhERd1pl0/QMTaCxdxR/sz7dK18/xKBDYVoUTnCdgTldNo1gxGKTlhikTiUR0Q4A/whgjxBCMbvK1mREo6l3lDutMgC4umKoNNvz9QWX0uypOsOoxTb3xczvVTinFd0588Md7i5XvZuI2gBsBPA6ER10vp5CRG9cc91zAD4AsJyI2ojoEeen/gNAOIBDRHSWiH7jTjyeUuz8n1HJdQYG4ERjH8KCdF4tAm7IiYXdIbhvEgMAVLQMIDpEj+y4UCn3n3Mq6UaEEPsB7J/h9Q4AO6/5+P5Z/v1cd+7vLYVpUdBpCBUtA7h1ZaLscJhkJ5v6UZIVDZ3WewPstZnR0GoIJxr7cPMy79XQmDpUtJixxkM9uRZDSauSFCPYoEVeSgQ/vTGYhi241DPi1WkkYKq7b0FqJG90YzCPTeJSz4jbZ364gxPDLIozolHVOgib3SE7FCbRKecvak9vbJvJhpxYnGszY2yS6wyBzNV5QVZ9AeDEMKs1GVEYt9pxkQ9rD2gnm/oQYtAiP9X73S1Lc2JgtQtUXOEVcYGs8soANASsTuPEoDjFGVPDuAouQAe0E419WJsZDb0X6wsuJZnR0NCHu6xZYKpoMWNFUgRCg9wqAbuFE8Ms0qKDkRAedHXZGAs8/aOTqO8e8UobjJmEG/VTdQbe6Baw7A6Bs61mqfUFgBPDrIgIxRnRqGjhYX2gOnW1P5L36wsupTmxONtqxoSVmzgGovruYYxYbCjOlNpomhPDjRRnRqGlfwymYYvsUJgEJxr7YdRrUJDquzdpaXYMJu0OnsIMUK7/766pbFk4MdwA1xkC24nGPpRkxsCg893bpCQrBhqaSkos8FRcMSM21ICMmBCpcXBiuIH81EjotcSJIQD1jlhwsWsYG5f4pr7gEhk8VWc4fqnXp/dlylDZMoDiTHkb21w4MdyAUa/FqpRIVPLywYDzvvMX8+bcOJ/fe1NuHCpbzXw+Q4DpH51EY++o9GkkgBPDnIozolHVZoaVN7oFlPcv9SLCqPPJ/oXrbc6Ng90hrm6uY4Gh8mp9QW7hGeDEMKfizChYbA5c6BiSHQrzESEEjjX04qYlcV45f2EuxZnRCNJpcIynkwJKRcsAdBpCocSNbS6cGObgWk/MdYbA0dw3ho7BCWxe6vtpJGBqCnN9dgyONXBiCCQVV8zIS4lAsMH3J7ZdjxPDHJIjg5EcaeT9DAHkWMPU0bEy6gsum3Pj0NAzgu6hCWkxMN+x2R2oajMror4AcGKYl+KMaN4BHUCOXepFalQwMmPlLRnc5ExK7/N0UkC42DWMsUm71MZ51+LEMA9rMqLQbh7np7cAYHcIHL/ch825cVKXDOYlRyAm1MDTSQGiUiEb21w4MczD1ToDjxr8XnX7IIYnbNLqCy4aDeGmJbE4dqkXCjgKnXlZRYsZCeFBSIsOlh0KAE4M87IqJRIGnYYL0AHAVV+4yccb22ayOTcOPc6Dgph/q2gZQLHEE9uux4lhHgw6DQpSI/lEtwBw7FIv8pIjEBsWJDuUq3UGXrbq33pHLLjSNya9cd61ODHMU0lmNGrahzA+yV0v/dXYpA0VV8zSp5Fc0mNCkBkbwnUGP1fePLWRcW2m77r4zoUTwzyV5kx1vazk6SS/dbp5AJN2h9RlqtfbnBuHE419vPPej51o7EewXovCNN/vsp8NJ4Z5utr1ktsU+K1jDSYYtBqsy1LOk9vm3DiMTtpR1cr7aPyVL08JnC/lRKJwEUY9VqVE4mQjH7vor45dmnqDKmHnqcvGJbEgAo7ydJJfMo9Noq572KeHQc0HJ4YFKM2OQSWfruWXekcsqO0cUkx9wSUqxIDC1Eje6OanTjb1Q4ipk/uUhBPDAmzIicWkzYGzPKz3O8cvT40ElVRfcOE23P7rZGM/gnQaRdUXAE4MC7IuOwZE4MPa/dCxBpO0Nttz4Tbc/stVXwjSKWf6EuDEsCCRwXrkJUfgBNcZ/IrsNttzKc6MhlGv4TqDnxkcs6K2awil2cqaRgI4MSxYaXYsKloGYLFxncFfXOoZQcfgBLYuj5cdyoyMei025MSirN4kOxTmQaebp+oLSis8A5wYFqw0JwYWmwPn2gZlh8I85IjzF+7Ny5SZGABg67J4NPaOoqVvTHYozENONPbBoNNgdbpydjy7cGJYoFJnneHEZZ5O8hdH6k1YmhCG1ChlNDCbyVZn0jpS3yM5EuYpJ5v6sSY9Cka9suoLgJuJgYjuJaLzROQgopIbXPcUEfUQUc0sn/8GEQkiUt6SkOtEhRiwPDEcJ7kQ6BfGJm042dh/9RevUmXHhSI9Jvjq6Iap29CEFec7BrFBYctUXdwdMdQAuAdA2RzXPQNgx0yfIKJ0ALcBaHEzFp/ZkBOLM1cGMGnjNgVqd6KxD5N2h2LrCy5EhG3LEnD8ch/Xt/xAeXM/HGJqalqJ3EoMQohaIUTdPK4rAzDbI/YvAPwDANU0nd+QE4Nxqx3V7byfQe2O1JkQrNcqqg3GbLYui8fYpB1nmrlfl9qdbOyHQatRzME815NaYyCiPQDahRBVMuNYqPXO5WUneD+D6h2pN2HjklhFzvNeb+OSWOi1hPd4Okn1TjT2oUih9QVgHomBiA4TUc0Mf/a6c2MiCgHwOIDvzvP6R4monIjKTSa5b4yY0Kk6A+9nULfm3lE0940pvr7gEhqkw7qsGByp48SgZsMTVtR0DCl2GgmYR2IQQmwXQuTP8OdVN++9BEA2gCoiagaQBqCCiJJmieMJIUSJEKIkPl7+G7k0JwZnrgxwO2QVcxVy1ZIYAGDb8njUdQ+jc3BcdihskcqvDMDuEIrc2OYibSpJCFEthEgQQmQJIbIAtAEoFkJ0yYppIUqzYzE2aUdNO+9nUKsj9SZkxoYgKy5UdijztnVZAgDwZjcVO9nYD72WFHVi2/XcXa56NxG1AdgI4HUiOuh8PYWI3rjmuucAfABgORG1EdEj7txXCVzDQK4zqNOE1Y4PLveparQAAMsSw5AUYeRlqyp2sqkPhWlRCDHoZIcyK3dXJe0XQqQJIYKEEIlCiNudr3cIIXZec939QohkIYTeef2+Gb5WlhBCNc1g4sKCkJsQhpNNXGdQo/LmAYxb7dim8GWq1yMibF0Wj6MNvbDxNKbqjFpsONc2iNJs5dYXAN757JaNObE41dTP+xlU6L26Hhi0GsVuMLqRbcvjMTxhQyW3f1edU039sDuE4n/uODG4YcvSOIxN2lHB50CrzpF6E9Znxyh6OD+bm3KnusDy6iT1KWswIUinwXoeMfivjUtiodMQFwJVpt08joaeEdXVF1wig/UozojiOoMKldWbUJqj/H0znBjcEG7UozgjGmUN/AZVE1ciV3objBvZuiwe1e2D6B2xyA6FzVO7eRyXTaO4WWHHx86EE4Obbl4Wh5r2IX6Dqsh7dT1IiTRiaUKY7FAWbdtyXraqNmUqaO/uwonBTa7/ycf4dC1VsNjsONbQi63L40GkvNPa5isvOQJxYUF49yK34VaLsnoTkiLU8UDCicFN+SmRiA7R83SSSpxo7MfopB3bVybKDsUtGg3h1hUJOFJn4lVxKmCzO/D+pV7cvCxOFQ8knBjcpNEQNi+dWlcuhGoaxAaswxe6EazXYlOu8ud553JbXiKGLTac4rNBFK+qbRBDEzZVTCMBnBg84ualcTANW1DbOSw7FHYDQggcru3GlqVxil8VMh+bcuNg1GtwuLZbdihsDmX1JhABm1XyQMKJwQNcTwE8naRs5zuG0Dk4ge156p5Gcgk2aLE5Nx6HLnTzaFXhyhpMKEyLQlSIQXYo88KJwQMSI4xYkRTOK0QU7tCFbhABt6xIkB2Kx9yWl4B28ziPVhVscMyKqlYztqpgmaoLJwYPuXlZPMqbBzA2aZMdCpvF4dpurM2IRlxYkOxQPOaWFYkgAk8nKdj7l3vhEOpYpurCicFDtiyNw6TdgZPcbVWROszjON8x5DfTSC7x4UFYkx6FQxc4MShVWb0J4UYditKV22b7epwYPGRdVgyMeg23KVCod5xP1GpfpjqT7XmJqG4f5MN7FEgIgbJ6EzYtiYNOq55ft+qJVOGMei1Ks2O5AK1Qb1/oRk5cKHJVsLlooT7uHAW9U8ub3ZTmsmkEHYMTqppGAjgxeNTNy+LRaBpF28CY7FDYNYYnrDjR2Od300guS+LDkBUbwnUGBTpSP9URYYuKCs8AJwaPcjXHOsrtMRSlrL4XVrvwy2kkYOrwnu0rE3H8Uh9GLLz4QUmONpiQExeK9JgQ2aEsCCcGD8pNCENypJGXrSrM4dpuRIdMtar2V9vzEjFpd+Ao/+wpxoTVjhONfaqbRgI4MXgUEWHb8qn2GBNWu+xwGKZ61Lx7sQe3rEhUVfFvoUoyoxEVoschnk5SjOOXezFhdaiyvbv/vlMk2VmQjBGLDe/x6VqKcLp5AIPjVtyW5z+b2mai02pwy/IEvHuxh8+CVojXqjoRGazHpiXqqi8AnBg8bmNOLGJDDXjtXIfsUBimppEMWg22LFXfU9tCbc9LhHnMijNX+KhZ2Sasdrx9oRu3r0qEQae+X7Pqi1jhdFoNduQn4Z3aHt4FLZkQAm9f6MJNubEIDVLf2c4LdfOyeBi0Ghw8z9NJsr1XZ8KIxYbdhSmyQ1kUTgxecOfqFIxb7byuXLKqtkG09o9jV0Gy7FB8IixIh63L4/F6dQccDm6qJ9Nr5zoQE2rATUtiZYeyKJwYvGBdVgwSwoN4OkmyA1UdMGg1+PiqJNmh+Mydq1PQPWTB6WZuzSLL2KQN79T24I78JNUueFBn1Aqn1RB2FiTjr3UmDE9YZYcTkOwOgdfOdWDr8nhEButlh+Mz21cmIFivxV+q+KFElncv9mDcalftNBLAicFr7lydjEmbg5ubSXK6uR/dQxbsWa3eN+dihBh02J6XiDdrumDl1UlSHKjqQHx4ENZnx8gOZdE4MXjJmvRopEYF47VznbJDCUgHqjoQrNfi1pX+vUx1JncWJqN/dBLHL/fJDiXgDE9Y8dc6E3YVJEOrUf7ZzrPhxOAlGg1hV2EyjjaYMDjG00m+ZLU78GZNF7bnJSLE4P+rka63dXk8wo06HODpJJ87XNuNSZsDuwvVveCBE4MX7S5MhtUucPB8l+xQAsr7l3rRPzoZcNNILkE6LW5flYSDNV28A9/HXqvqRHKkEcUZ0bJDcQsnBi8qSI1ERkwIDvDqJJ86UNWJcKMONy9T345TT9mzOgXDFhufD+JDg2NWlDWYsLswGRoVTyMBbiYGIrqXiM4TkYOISm5w3VNE1ENENTN87stEVOf8Oj92Jx6lISLcuToZxy/3oW/EIjucgDBhtePt813YsSoJQTqt7HCkuWlJLGJCDTyd5EMHz3fBaheqXo3k4u6IoQbAPQDK5rjuGQA7rn+RiD4GYC+AQiHEKgA/dTMexdldmAK7Q+DNGp5O8oUj9SYMW2zYU6T+N6c7dFoNdhbwDnxfOnCuAxkxIShMi5QditvcSgxCiFohRN08risDMNOOmy8B+DchhMV5nd9tFV6RFI4l8aH85OYjf6nqQGyoARtz1Lnj1JPuLJzagc9Lpr2vb8SC45f7sKswGUTqnkYC5NcYlgHYQkQniegIEa2THI/HERH2rE7FqeZ+dJj5TF5vGrXY8E5tN3YWJKt2x6knrcuKQVKEEQeqeMm0t71e3Qm7Q+BOP5hGAuaRGIjoMBHVzPBnrwfurwMQDWADgP8D4EWaJd0S0aNEVE5E5SaTugpqd61JgRDg3ahedri2GxNWR8BPI7loNITdhck4Ut/DS6a9bH9lO1YkhSMvJUJ2KB4xZ2IQQmwXQuTP8OdVD9y/DcDLYsopAA4AMy4lEUI8IYQoEUKUxMerq4VyZmwo1mZGY39FO4Tg5mbecqCqA8mRRqxV+VJBT7pzdQovmfay5t5RVLaYcfeaVNmheIzs8fYrAG4BACJaBsAAwC8PTL5rTSrquodR2zksOxS/ZBq24L06E/asTlH9UkFPKkyLRHZcKP5c0SY7FL+1v7IdRPCrkaq7y1XvJqI2ABsBvE5EB52vpxDRG9dc9xyADwAsJ6I2InrE+amnAOQ4l7E+D+Ah4aeP1LsLkqHTEF452y47FL/0SmU7bA6Be0vSZYeiKESEe0vScKqpH029o7LD8TtCCLxyth0bc2KRHBksOxyPcXdV0n4hRJoQIkgIkSiEuN35eocQYuc1190vhEgWQuid1+9zvj4phHjAOTVVLIR4171vR7miQw3YtjwBr55th5175XuUEAIvlLdibWY0chPCZIejOJ8oToOGgD+Vt8oOxe9UtppxpW8Md/nRNBIgfyopoNy9JhXdQxZ8wM3NPKqy1YxLPSO4ryRNdiiKlBhhxLblCXipoo3Pg/awVyrbEaTT4I58/zrzgxODD926MgHhQTrsr+TpJE/6U3krQgxa7PKTpYLecF9JOrqHLDja4JclPCmsdgcOVHXgtrxEhBv968wPTgw+ZNRrsbMgGW/VdGJ8kpubecLYpA0HqjqxqyAZYQFwrvNi3bIiAbGhBrxwmqeTPKWs3oSBMatfrUZy4cTgY3etScXopB2Hank3qie8Ud2FEYsN963jovONGHQa3L0mFYdru7lvl4fsr2xHdIgeNy9T1/L5+eDE4GOl2TFIjjTiFZ5O8ogXy1uRExeKkkzeuzCX+9alw+YQPJXpAUMTVhy60I07V6dA74e77P3vO1I4jYawtygVR+pN6OUnN7c09Y7iVFM/7i1J94v+NN62LDEcRelReLG8lTdauumtmi5YbA6/W43kwolBgnuKU6cOq+cWGW75U3krtBrCJ4r9883pDfeVpKO+ewRVbYOyQ1G1VyrbkRUbgjXpUbJD8QpODBIsSwxHXnIEXuYh/aLZ7A68VNGGbcvikRBhlB2OauxenQyjXoMXeU/DonWYx/FBYx/uWpPqtyNVTgyS3FuShnNtgzjfwU9ui1HWYEL3kIV3Oi9QhFGPnfnJOHC2g1fGLZIrqX6i2H/3zXBikOTuNakI0mnw/Cl+cluMF063Ii7MgFtXJsgORXXuW5eOYYsNb1RzO+6FsjsEXjjdii1L45EeEyI7HK/hxCBJVIgBuwqS8UplO5+wtUBa17sDAAASgUlEQVTt5nEcru3BJ4rT/HJFiLeVZscgJy4Uv/2gmYvQC/ReXQ86Byfw6fX+PVLld5VE95dmYNhiw2t8kMqC/PaDZgDAZ2/KkhmGahERPrcpC1Vtg6hoMcsOR1WeO9WC+PAg3LoyUXYoXsWJQaKSzGgsTQjDH0+1yA5FNcYmbXj+VCtuX5WI1Cj/6Wbpa58oTkO4UYen32+SHYpqdA6O492LPbivxP9Hqv793SkcEeH+9Rk422rGhY4h2eGowssV7Rgct+LhTdmyQ1G10CAd/mZdOt6s6eIjZ+fphdOtcAjgb9ZlyA7F6zgxSHZPcSoMOg2e41HDnBwOgaffb0JhWiTW8k5nt312YxaEEPjdiSuyQ1G8D4vOcX5ddHbhxCAZF6Hn7+ilXlw2jeLzm7L8dv24L6XHhODjeUn448kWXro6hyP1rqKz/48WAE4MivBpVxH6HBehb+SpY02IDw/CrgJur+0pD2/OxuC4FS9X8tGfN/LHk62ICwvC9jz/Ljq7cGJQgBLnyWN/PMnTSbO51DOCI/UmPLghEwYd/9h6yrqsaKxKicDT7/PS1dlMFZ27A6Lo7BIY36XCcRF6bs8cb4JBp8GnSwNjKO8rRISHN2XjUs8IH+IzixdPtwVM0dmFE4NCfIKL0LMaHLPipTPtuKsoBXFhQbLD8Tu7VycjLiyIl67OYKro3IItS+OQEev/RWcXTgwKERViwO6CZOyvbMfQhFV2OIry/OkWjFvt+DwvUfWKIJ0WD2zIwF/rTLhsGpEdjqIcru1GRwAVnV04MSjI5zdlY8Riw4t8/OJVNrsDv/3gCjbkxGBlcoTscPzWZ0ozYdBqeNRwnX1Hm5AaFYzbAqTo7MKJQUEK0iKxPjsGT7/fDJvdITscRXizpgvt5nF8YXOO7FD8Wnx4EPYWpeClM+0wj03KDkcRzrWZcaq5H5/flAVdgBSdXQLru1WBRzZno908joPn+UxoIQSePNaE7LhQ3LKCu6h62yNbsjFutXOLFqd9x5oQFqTDpwLwPHFODAqzfWUiMmJCsO9Yo+xQpKtoGUBVqxkPb86GRsMb2rxtRVIEtiyNw7PHmzFpC+wRa+fgOF4/14lPrUtHuFEvOxyf48SgMFoN4eFNWahoMaOiZUB2OFI9ebQJkcF6PrrThx7enI3uIQterw7sY2efPX4FDiHwuQDt4MuJQYHuLUlHuFGHfccCtxDY0jeGg+e78JnSDIQYdLLDCRhbl8YjNyEMTx5tCtgNb6MWG/548gp25CcFRF+kmXBiUKDQIB0+vT4Db9V0oW1gTHY4Ujx9vAlaDeGhAH1ik0WjITyyORvnO4ZwsqlfdjhSvFTRhqEJGx7ZHLjLozkxKJTrF+Kzx5ulxiHD4LgVL55uxZ2FKUiMMMoOJ+DcvSYVMaEGPHk08EasDofAU8eaUJQeheKMwO3gy4lBoVKignFHfhKeP9WKEUtgdV194XQLRifteDiAn9hkMuq1eGBDJt652I3GANvw9s7FHjT3jeGRzdkB3cHXrcRARPcS0XkichBRyQ2ue4qIeoio5rrXi4joBBGdJaJyIlrvTjz+5gtbcjBsseFP5YGz4c1md+CZ95uxIScG+amRssMJWA9uyIReo8HT7zfLDsWn9h1rRKrzoSyQuTtiqAFwD4CyOa57BsCOGV7/MYDvCSGKAHzX+TFzKkqPQklmNJ4oa8RogIwa3qzpQsfgBG9ok8y14e3PZ9oCZsPbB5f7cKKxHw/dlBlwG9qu59Z3L4SoFULUzeO6MgAzVbIEAFefg0gAgb1Gbgbf3LkCnYMT+OXhetmheJ3DIfCbI5d5Q5tCuDa8PRMAdS6LzY7HX6lGekwwHtyQJTsc6WSnxa8B+AkRtQL4KYBvSo5HcdZmxuD+9Rl46v1mv2/J/crZdpzvGMKXb8nlDW0KsCIpAjtWJeGJskb0DE3IDsernjjSiEbTKL6/Nx/BBq3scKSbMzEQ0WEiqpnhz14P3P9LAP63ECIdwP8GsO8GcTzqrEOUm0wmD9xaPR7bsQJRwXp8a3817A7/XFs+PmnHTw7WoSA1EncV8YY2pXjsjhWw2h34+SH/HbE2947iV3+9hF0FyfjYch6pAvNIDEKI7UKI/Bn+vOqB+z8E4GXn3/8EYNbisxDiCSFEiRCiJD4+3gO3Vo/IED2+vXslzraa/fa8hn3HGtE5OIHHd63k0YKCZMWF4sENWXixvBUXu/xvxCqEwHderYFBq8F378yTHY5iyJ5K6gCw1fn3WwA0SIxF0e4qSsWm3Fj86K2L6Bn2r2G9adiCX793GbflJWJDTqzscNh1vnJrLsKNevzwjYuyQ/G418514mhDL77x8WW8Z+Ya7i5XvZuI2gBsBPA6ER10vp5CRG9cc91zAD4AsJyI2ojoEeen/hbAz4ioCsAPATzqTjz+jIjwg735sFgd+OfXamWH41G/OFwPi82Bb96xQnYobAZRIQZ8+ZZclNWbcKTef6ZxB8et+P5rF1CQGokHN2bJDkdR3F2VtF8IkSaECBJCJAohbne+3iGE2HnNdfcLIZKFEHrn9fucrx8TQqwVQqwWQpQKIc649+34t5z4MHxp2xL8paoDZX7yBq3vHsbzp1rwwIZM5MSHyQ6HzeKzG7OQGRuCH75e6zd1rp8erEPfiAU/vLsAWp6+nEb2VBJboC9tW4LsuFB859UaTFjtssNx2w/fqEVokA5fuXWp7FDYDRh0Gjy2YwXquofxoh9suKxsGcDvT17BZzdmoSCNN1JejxODyhj1WvzL3fm40jeGXx5Wd0nmaIMJ79WZ8OVbchETapAdDpvDjvwkrMuKxs/erld1m5ZJmwOPvVSNpAgjvnH7ctnhKBInBhW6aUkcPlWSjv8+2ojzHYOyw1kU89gkvvNKDdJjgrmDqkoQER7flYfeEQt+cOCCattyP1F2GXXdw/jnu/IRFsQt3WfCiUGlvrVzJaJDDHjspWrVnQ89aXPgi787gw7zBH5+XxGCdLyhSC2K0qPwdx/LxQvlrfjvo+o7ZfCyaQT//s4l7CpMxq0rE2WHo1icGFQqMkSP7+1Zher2QVW1LBBC4JsvV+NkUz9+/MlCrMuKkR0SW6Cv37YMuwqT8a9vXsRbNV2yw5k3h2PqZy/YoMU/3blKdjiKxolBxXYWJGH7ygT87O16tPar40Cf/3zvMl6qaMNXb12Ku9bwDmc10mgIP7t3NYrSo/C1FypR1WqWHdK8PH+6Faea+vH4zpWIDw+SHY6icWJQMSLC9/fmQ0PAt/ZXK37O90BVB35ysA53FaXga9t5FZKaGfVa/PdnSxAXFoQv/LYc7eZx2SHdUPfQBP71zVpszInFvSVpssNRPE4MKpcSFYx/vGMFjjb0Yn9lu+xwZlXRMoC//1MV1mVF40efLAzoQ1D8RVxYEJ7+3DpMTNrxyDOnMTxhlR3SrP7vq+cxaXPgh/cU8M/ePHBi8AMPlGaiOCMK33/tAroV2AWzbWAMj/62HMmRRvzXgyVcbPYjSxPD8esH1qKhZwRfea5SkZvfDlR14K3zXfjq9qXIjguVHY4qcGLwAxoN4af3rsaE1Y5/fOmcoqaURiw2fOHZclhsDux7aB3vV/BDm5fG4ft7V+GvdSb88A1ltWvpHprAt1+pQVF6FB7dwoc/zRcnBj+REx+Gb+1ciffqTPijQjqw2h0CX32uEg09I/jPzxQjN4FbXvirz5Rm4nM3ZWHfsSbFdAAWQuAf/nwOFpsdP79vdcCfyrYQ/F/KjzxQmoktS+Pwz6/Vorl3VHY4+NFbF/HOxR7805152LI0sFqlB6Jv71qJrcvi8Z1XanD8cq/scPCHky04Um/Ct3au5D5cC8SJwY9oNIQff7IQOi3h7/9UJXW+94XTLXiirBEPbczkzpUBQqfV4FefXoPsuFB86fcVaJL4cNLcO4p/eb0WW5bG4YHSTGlxqBUnBj+THBmMH+zNx5krA/ivsstSYjjR2IfH99dgy9I4fGc3H34SSCKMeux7aB00BDzyzGkMjvl+pZLdIfD1F89Cr516UOKDnxaOE4Mf2luUgl0FyfjFoXqfnxN9pN6Eh585jczYEPzHp4t5XjcAZcSG4L8eLEHrwBg+9cQHPl8p95sjl1HRYsYP7spHcmSwT+/tL/hd64eICD+4Kx9RIQb87W/LcdZHO1NfrmjDI8+cRlZsKJ772w2IDNb75L5MedZnx2DfQ+vQ2j+Ge/7zOC71jHj9nkIIPPN+E35xqB67CpKxZ3WK1+/przgx+KmYUAOeemgdAODe3xzHvmNNXlvGKoTAr9+7jK+/WIXSnBi88MUNSOBjEgPezcvi8cIXN8Jis+OTvzmOM1cGvHavoQkr/ucfKvBPBy7g5mXxvJHNTZwY/FhBWiTe+MoWbFuegB+8dgFf/N0Zj8/52h0C3ztwAT966yL2rE7B059bj3AjjxTYlPzUSLz0pZsQFazHZ548gUMXuj1+j+q2Qez+92N4+0I3vrVzBZ78bAmPVt3EicHPRYbo8cSDa/HtXSvx7sUe7PrVURy/1OuR0UPbwBi+8OxpPHO8GV/YnI1ffqoIBh3/SLHpMmND8ecv3YTlieH44u/K8YtD9R45fXB80o59x5rwiV8fh83uwItf3IBHb17CxWYPICXtkp2vkpISUV5eLjsM1alsGcDf/bES7eZxpEUHY8/qFOwtSsXypPAFfR2LzY4njzbhV+82gEB47I4VfNgOm9OoxYZvvlyNv1R1ICMmBN/bswofW5GwoK9hsztw7FIv/nK2AwfPd2F00o5bViTgZ/euRjTvqp8TEZ0RQpTMeR0nhsAyarHh4PkuvHq2A8cu9cLuEFiRFI5NuXGIDw9CbKgBceFBiAsNQnSoHhHBeoQZdFefwo419OK7r9agsXcUO1Yl4Tt35iE1ild+sPl7/1IvvvNqDRpNo/h4XiK+e2ce0qJDAABWuwPDEzYMjVvRN2pB78gkekcs6BuZRPvAOA7XdqNvdBIRRh12FiRjT1EKNubEcj1hnjgxsDn1jljw+rlOvHq2HRc6hzBhnfkkOA0B4UY9Qg1adAxOIDN26mlv2/KFPe0x5jJpc+DJY4341TuX4BAC0SEGDE1YMTY5+xRTZLAem3PjsKcoBduWx3MzxkXgxMAWRAiBsUk7ekc+fEobHLNiaMKKoXErhiZsGBy3YlliOD6/KQtGPb8pmfvazeP4zXuXYbHZEWGcGqFGGHUIN+oRE2ZAfFgQYsMMiA0N4vqVB8w3MfBJ2AzA1N6H0CAdQoN0yIzl1sTMN1KjgvGDu/Jlh8GuwymYMcbYNJwYGGOMTcOJgTHG2DScGBhjjE3DiYExxtg0nBgYY4xNw4mBMcbYNJwYGGOMTaPKnc9EZAJwRXYcixAHQP4p6b7H33dg4e9buTKFEPFzXaTKxKBWRFQ+n+3o/oa/78DC37f68VQSY4yxaTgxMMYYm4YTg289ITsASfj7Diz8fasc1xgYY4xNwyMGxhhj03BikISIvkFEgojiZMfiC0T0EyK6SETniGg/EUXJjslbiGgHEdUR0SUiekx2PL5AROlE9FciqiWi80T0Vdkx+RIRaYmokohekx2LJ3BikICI0gHcBqBFdiw+dAhAvhCiEEA9gG9KjscriEgL4P8BuANAHoD7iShPblQ+YQPw90KIlQA2APhfAfJ9u3wVQK3sIDyFE4McvwDwDwACpsAjhHhbCGFzfngCQJrMeLxoPYBLQohGIcQkgOcB7JUck9cJITqFEBXOvw9j6pdkqtyofIOI0gDsAvCk7Fg8hRODjxHRHgDtQogq2bFI9DCAN2UH4SWpAFqv+bgNAfIL0oWIsgCsAXBSbiQ+80tMPeg5ZAfiKXzmsxcQ0WEASTN86nEA3wLwcd9G5Bs3+r6FEK86r3kcU9MOf/BlbD5EM7wWMCNDIgoD8BKArwkhhmTH421EtBtAjxDiDBFtkx2Pp3Bi8AIhxPaZXieiAgDZAKqICJiaTqkgovVCiC4fhugVs33fLkT0EIDdAG4V/rtOug1A+jUfpwHokBSLTxGRHlNJ4Q9CiJdlx+MjmwDsIaKdAIwAIojo90KIByTH5RbexyARETUDKBFCKL3xltuIaAeAnwPYKoQwyY7HW4hIh6ni+q0A2gGcBvBpIcR5qYF5GU096TwLoF8I8TXZ8cjgHDF8QwixW3Ys7uIaA/OV/wAQDuAQEZ0lot/IDsgbnAX2vwNwEFMF2Bf9PSk4bQLwIIBbnP9/zzqfopkK8YiBMcbYNDxiYIwxNg0nBsYYY9NwYmCMMTYNJwbGGGPTcGJgjDE2DScGxhhj03BiYIwxNg0nBsYYY9P8f3D9TswVzU51AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import scipy as sp\n",
    "import scipy.linalg\n",
    "U = np.linalg.matrix_power(AM, Nt)\n",
    "H = 1j*h*sp.linalg.logm(U)/T\n",
    "assert np.allclose(H, H.T.conj())\n",
    "#assert np.allclose(H.imag, 0)\n",
    "#H = H.real\n",
    "plt.plot(x, np.diag(H))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0xa1a206310>]"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAHtpJREFUeJzt3Xl4VPWh//H3N/sCWUgChARIIGEJOwQQcENF0Vp3rLbutqhXbf3Z5ept7a33Z59qa2ttr1eLK637UpVqBRE3qiwB2cK+BEggIRvZ95nv/SNjn1wLEgiTMzPn83qePJNzZsj5HJj55PCd75xjrLWIiEjoC3M6gIiI9A4VvoiIS6jwRURcQoUvIuISKnwREZdQ4YuIuIQKX0TEJVT4IiIuocIXEXGJCKcDdJWammqzsrKcjiEiElTWrl1baa1NO9bjAqrws7KyWLNmjdMxRESCijFmX3cepyEdERGXUOGLiLiECl9ExCVU+CIiLqHCFxFxCRW+iIhLqPBFRFxChS8i4hBrLdvL6nnik918vqvS79sLqA9eiYiEuqa2DpbvrOTj7eV8vL2C0toWAG49Yzgzc1L9um0VvohIL9lVXs/1zxRwoKaZPtERnJqTyg/OTuOMkWmkJ8b6ffsqfBGRXlCwt5rvLlxDZHgYz904lZnDU4mK6N1RdRW+iIifLS4s4wcvryMjKZaFN01jcL84R3Ko8EVE/OgvK/by80WbmTg4iaevn0q/+CjHsqjwRUROEmsth+pa2XSglk0lNawrrmH5zkrOGd2fP149mdiocEfzqfBFRE6CFburuOuVdRyqawUgzMCIAX25Y3YOd52TS0S487PgVfgiIifBo8t2YDD84pt5jMtMIi89wfEj+q9S4YuI9ND+qiZW7qnmR+eO4IZZ2U7HOSrn/48hIhLkXl9bjDFw+ZRMp6N8LRW+iEgPeLyW19eWcHpu73x4qidU+CIiPfDZrkoO1rZwZf5gp6MckwpfRKQHXltbQlJcJOfk9Xc6yjGp8EVETlBNUxtLNpdxycQMoiMCa0bOkajwRURO0KINB2nr8DIvP7DfrP2SCl9E5AS9tqaEvPQExgxKdDpKt6jwRUROwJaDdWw6UMuVQXJ0Dyep8I0xzxhjyo0xhV3W9TPGLDXG7PTdJp+MbYmIBILX1hYTFR7GxRMznI7SbSfrCP85YO5X1t0DLLPW5gLLfMsiIkGvrcPLW+sOMGfMAJIdPPvl8TophW+t/RSo/srqi4GFvu8XApecjG2JiDjtzyv2cripnXkB/snar/LnGP4Aa20pgO828Cepiogcwxf7D/Pge9uYkzeAM0akOR3nuDj+pq0xZr4xZo0xZk1FRYXTcUREjupwYxt3vPAF6UkxPHzFBIwxTkc6Lv4s/EPGmHQA3235kR5krV1grc231uanpQXXb0sRcQ+v13L3q+upbGjjsW9PJjEu0ulIx82fhb8IuN73/fXA237cloiIXz3x6W4+2l7BfReOZnxmktNxTsjJmpb5ErACGGmMKTHG3Aw8CMwxxuwE5viWRUSCzqo9VTy8ZDsXjk/nmlOGOh3nhJ2UC6BYa68+yl1nn4yfLyLilKLKRu58aR1DU+L51WXjgm7cvitd8UpE5CjWF9dw03MFADx+zWT6xgTfuH1Xjs/SEREJRB9uO8TVC1YSHx3OG7fNZNTABKcj9ZiO8EVEvuKVgv38x5uFjE7vy7M3TCOtb7TTkU4KFb6IuFpLu4dDdS2U1rZQVtvCuv2HWbhiH6flpvL4NVPoEx06NRk6eyIicpy+/9I6Fm04+C/r503J5JeXjiMqIrRGvVX4IuJKxdVNLNpwkLljBnJO3gDSE2MYmBjDwIQY4kPoqL6r0NwrEZFjeKWgmDADP/9mHoOSYp2O0ytC6/8rIiLd0OHx8traYs4YkeaasgcVvoi40Cc7KjhU18pV04Y4HaVXqfBFxHVeWl1Map9ozhrlrrO2q/BFxFUO1bXw0fZyrpiSSWS4uyrQXXsrIq73+toSPF7Lt6YOdjpKr1Phi4hreL2WV9cUc8qwfmSnxjsdp9ep8EXENVbuqWJfVRNXTXXXm7VfUuGLiGu8XFBMQkwEc8cOdDqKI1T4IuIKhxvbWFxYxqWTMoiJDHc6jiNU+CLiCs+v3Eebx+u6ufdd6dQKIhLSOjxeHlq8jSeXFzF7ZBqj04P/vPYnSoUvIiGrqqGVO15cx4o9VVw/Yyg//Uae05EcpcIXkZC0obiG255fS1VjGw/Pm8AVUzKdjuQ4Fb6IhJQDNc38ecVenv1sL2l9onnjtpmMzUh0OlZAUOGLSNCz1rKqqJrnPtvL+1vKMMZwwbh07r9oDP3io5yOFzBU+CISdNo9XnaVN7D5YB2FB2pZsbuK7YfqSYqL5JYzhnPNKUPJcNFpj7tLhS8iAaPD46W2uZ3DTe3UNrdxuLGdioZWDtW1UF7fSnldCwdrWthV0UBbhxeA2MhwxmYk8NDl47h4onvn2HeHCl9EHGet5c6X1vHOxtKjPiYlPor+CTEMSIjmtNxU8gYlMGZQItmp8YSHmV5MG7xU+CLiuFfXFPPOxlKuzM8kLz2B5PgoEmMjSY6LIq1vNKl9okPuguJOUOGLiKPKalt44J2tTM/ux4OXjSdMR+t+o1+ZIuIYay0/e2sT7V4vD12usvc3Fb6IOGbRhoN8sLWcH84ZSZYLz0/f21T4IuKIqoZW7v/bFiYMTuKmU7OdjuMKKnwRccR/LtpMfUs7v7livGbZ9BK/v2lrjNkL1AMeoMNam+/vbYpIYKpvaafwQB3/2FXBOxtLuXvOCEYM6Ot0LNforVk6s621lb20LREJIHsqGvjjh7vYUFLDnorGf66fMSyFW88Y7mAy99G0TBHxm4r6Vq59ejV1ze2cMjyFSydmMC4zkXEZiaT0iXY6nuv0RuFb4H1jjAX+ZK1d0AvbFBGHtbR7mP+XNVQ1tvLqLTMYn5nkdCTX643Cn2WtPWiM6Q8sNcZss9Z++uWdxpj5wHyAIUPce+kxkVDi9Vp+9NoG1u2v4YlrJqvsA4TfZ+lYaw/6bsuBN4FpX7l/gbU231qbn5aW5u84ItILHvlgB+9sLOWe80cxd2y603HEx6+Fb4yJN8b0/fJ74Fyg0J/bFBFnvbG2hD9+uItv5Q/mltOHOR1HuvD3kM4A4E1jzJfbetFau9jP2xQRh2w5WMc9f93IzOEp/P9LxuJ77UuA8GvhW2v3ABP8uQ0RCRy/em8r8dER/M93JuvslgFI/yIiclL8Y2cly3dWcsfsHJLidFnBQKTCF5Ee83otDy3eRkZSLNfOGOp0HDkKFb6I9Ni7m0rZdKCWu+eMIDpClxgMVCp8EemRdo+Xh9/fzqiBfblkUobTceRrqPBFpEdeXr2ffVVN/PvcUTrrZYBT4YvICWts7eDRZTuZnt2PM0fqg5OBToUvIifsqeVFVDa0cc/5ozTnPgjobJki0m3WWnaWN7ByTxUr91SxbGs5c8cMZNKQZKejSTeo8EXkmNo6vPzXO5t5b1MZVY1tAAxKjOHC8YP4ydyRDqeT7lLhi8jXavd4uf3FL1i65RAXTRjEqbmpzBiWQmZyrIZxgowKX0SOqt3j5c4X17F0yyHuv2gM18/McjqS9IDetBWRI+rweLnr5fUs3lzGfRfmqexDgApfRP5Fh8fL/3t1A+9uKuWnF4zm5lOznY4kJ4GGdETkn1o7PHy0rZyFn+9jxZ4q7jl/FN/TOe1DhgpfxOWstWw6UMsba0t4e8NBapra6d83mv+6eAzXzchyOp6cRCp8ERcrrW3mJ69vZPnOSqIiwjhvzEAun5zBqTmpRIRrxDfUqPBFXMhay9vrD3Lf24V0eCw/+8Zo5uUPJjE20ulo4kcqfBGXqW5s42dvbeLvm8qYMjSZ386bQFZqvNOxpBeo8EVcwFrLjkMNLN1SxsIV+6hpauPf545i/unDdIZLF1Hhi4Qoay2ri6pZsvkQS7eWUVzdDMDUrGTuv3EaeYMSHE4ovU2FLxKCWjs8/OzNQl5bW0JURBizhqdw2xk5nD26PwMSYpyOJw5R4YuEmIr6Vm57fi1r9h3mzrNyuPWM4cRH66UuKnyRkLL5YC3fW7iG6qY2Hvv2ZL4xPt3pSBJAVPgiIeK9TaXc/eoGkuIief3WmYzNSHQ6kgQYFb5ICPh8dyX/9uIXTBqcxBPXTqF/X43Ty79S4YsEucONbdz9ygayU+J5/rvTiYvSy1qOTJ+dFgli1lru/esmqhpb+cPVk1T28rVU+CJB7JWCYhZvLuNH547UmL0ckwpfJEjtrmjg/r9tYVZOCt87TacwlmNT4YsEobaOzqtRRUeG8dt5EwnT6RGkGzTgJxJE6lva2VZWzysFxWw6UMufrp3CwETNyJHu8XvhG2PmAo8C4cBT1toH/b1NkVDQ7vGysaSW1UXVbCypYUtpHfuqmv55/w0zszhvzEAHE0qw8WvhG2PCgceAOUAJUGCMWWSt3eLP7YoEqwM1zbyxtoTVRdWs3XeY5nYPAFkpcYwZlMC8KZnkDUpgdHoC6YmxDqeVYOPvI/xpwC5r7R4AY8zLwMWACl/kK5bvrOCOF9dR19LOyAF9+dbUwUzP7sfU7H6k9ol2Op6EAH8XfgZQ3GW5BJje9QHGmPnAfIAhQ4b4OY5I4LHW8tTyIn713lZy+vfhrdtnka0Lkogf+LvwjzR1wP6fBWsXAAsA8vPz7REeLxKymts83PPXjby9/iBzxwzkt1dO0JktxW/8/cwqAQZ3Wc4EDvp5myJBobS2me8uXMOW0jp+dO4Ibp+dgzGaXin+4+/CLwByjTHZwAHgKuDbft6mSMDbU9HAtU+vpra5naeuy+fs0QOcjiQu4NfCt9Z2GGPuAJbQOS3zGWvtZn9uUyTQbT5Yy/XPrMZaeHn+KTolgvQavw8WWmv/Dvzd39sRCQZr9lZz43MF9I2O4C/fnc7wtD5ORxIX0btDIr3k4+3l3Pr8WgYlxvKX704nI0nz6KV3qfBF/Ky0tpnHP97Ni6v2M3JgXxbeNE3z6sURKnwRPzlQ08zjH+/i1YISvNYyLz+Tey8YTUJMpNPRxKVU+CInWVltC3/4cCevren8zOG8/MH825nDyUyOcziZuJ0KX+QkqW1q5/FPdvPsZ0V4reVbUwdz25k5GquXgKHCF+mhlnYPz32+l//5aBf1rR1cMjGDu+eMYHA/HdFLYFHhi/RAWW0LVz+5kqLKRs4cmcZPzhtF3qAEp2OJHJEKX+QElde38O0nV1JR38pfbp7GablpTkcS+VoqfJETUNXQyjVPraKsroWFN01jalY/pyOJHJOuaStynGqa2rjm6dXsq2riqevzVfYSNFT4Isehtrmd655Zze7yBp68Lp+Zw1OdjiTSbRrSEemirqWdDcU1rNtfw7r9h9lb1URTWwfNbR5a2r20ebxEhhueuGYKp4/QmL0EFxW+uFpzm4eVRVV8sr2Cz3ZVsquiAWvBGMjt34e8QQn0iYogNiqc2Khw4iLDmTE8hXwN40gQUuGL63R4vLxUUMzSLYdYtaeK1g4vMZFhTM9O4aIJg5g0JJnxgxN1CgQJOSp8cZ0nPtnNw+/vYFhaPN+ZPpQzRqYxPbsfMZHhTkcT8SsVvrjKzkP1/GHZLr4xPp3Hvj3Z6TgivUqzdMQ1PF7Lj1/fSHx0OPdfNMbpOCK9Tkf44hrPflbE+uIaHr1qos5HL66kI3xxhaLKRn6zZDvnjB7ARRMGOR1HxBEqfAl5Xq/l39/YSFREGL+8dCzGGKcjiThChS8h74VV+1hdVM19F+YxICHG6TgijtEYvoSkxtYOlm0r550NB/loezmn5aYyb0qm07FEHKXCl5BhrWXZ1nLe+KKED7eV09rhpX/faL4zfSh3npWjoRxxPRW+hIQv9h/ml+9uZe2+w6T1jeaqqYP5xvhB5A9NJixMRS8CKnwJcvuqGvn14u28u6mUtL7R/OqyccybkklEuN6eEvkqFb4ErYWf7+WBd7cQERbGXefk8r3ThhEfrae0yNHo1SFBqaapjYcWb2Nadj8euXIi/TX7RuSY9P9eCUovrNpPU5uH+y7MU9mLdJMKX4JOa4eHZz/by+kj0hg1MMHpOCJBQ4UvQeetdQeobGjlltOHOR1FJKj4rfCNMb8wxhwwxqz3fV3gr22Je3i9lieXF5GXnsDM4SlOxxEJKv5+0/YRa+3Dft6GuMhH28vZVd7Ao1dN1AepRI6ThnQkqCz4dA+DEmO4YFy601FEgo6/C/8OY8xGY8wzxphkP29LQtyG4hpWFVVz06nZROqDVSLHrUevGmPMB8aYwiN8XQw8DgwHJgKlwG+P8jPmG2PWGGPWVFRU9CSOhLgFy/fQNyaCq6YNcTqKSFDq0Ri+tfac7jzOGPMk8M5RfsYCYAFAfn6+7UkeCU0dHi+7Khp4b1Mp808fTh99mlbkhPjtlWOMSbfWlvoWLwUK/bUtCS2ri6p58L2tVDa0cbipjfqWDgAiww03zMxyNpxIEPPnodKvjTETAQvsBW7x47YkROw4VM/NCwtIiIlkalYySXFRJMZGkhwXyZSh/RiYqE/VipwovxW+tfZaf/1sCU2H6lq48dkCYiLDeeWWU8hMjnM6kkhI0VQHCQgNrR3c9FwBh5vaePaGqSp7ET/Qu1/iuHaPl9tf+IJtZfU8dX0+YzMSnY4kEpJ0hC+O8not971VyCc7KnjgkrHMHtnf6UgiIUtH+OKYjSU1/OeizazbX8Mds3O4WvPrRfxKhS+9rrqxjd8s2cbLBcWkxEfz8LwJXD45w+lYIiFPhS+9pqXdw8ur9/O7pTtoavNw86xsvn9OLgkxkU5HE3EFFb74XVNbBy+u2s+CT/dQXt/KrJwUfvHNMeQO6Ot0NBFXUeGL39Q2tfPnFXt55rMiDje1M2NYCo98ayIzh6fo1MYiDlDhy0nT2uFh3f4aPt9VyWe7q9hQXEOH13LWqP7cPjuHKUN1wlQRJ6nwpccO1jTz8PvbeW9TGc3tHsIMjMtI5HunD+Ob4weRN0jXnRUJBCp8OWENrR088fFunly+BwvMm5LJ6SPSOGVYComxeiNWJNCo8OW4ebyW19YU8/D7O6hsaOXiiYP48XkjdToEkQCnwpfjsnZfNfe9tZktpXVMGZrMk9dNYdIQjc2LBAMVvnRLRX0rD763jTe+KGFgQgx/vHoSF45P12wbkSCiwpevVdvUzhtflPDI0h20dHi47czh3DE7h3hddUok6OhV62L1Le3Ut3Tg8drOL2tpavWw6UAt6/Yf5ov9h9ld0QjAabmp/OKiMQxP6+NwahE5USp8l+nwePlkRwWvFBTz4bZyOrxHvoxwclwkk4Ykc+mkDKZlpzA1K1nDNyJBToXvEgdqmnlh5T5eX1tCeX0rqX2iuHFWFsPT+hAWZogIM4SHGaLCwxiVnkBWSpwKXiTEqPBdoLG1g0se+4yqhlbOHNmfK/MHc/bo/kSG63IIIm6iwneBp5YXUVHfymu3zmBqVj+n44iIQ3SIF+KqGlpZ8OluzhszQGUv4nIq/BD33x/torndw4/PG+l0FBFxmAo/hBVXN/HCyv3MmzKYnP4697yI26nwQ9gjH+zAGLhrTq7TUUQkAKjwQ9S2sjreXHeAG2ZmkZ4Y63QcEQkAKvwQ9ZvF2+kbHcFtZw53OoqIBAhNywwB1lpa2r00tnXQ3Oah8EAty7aV85O5I0mKi3I6nogECBV+EKlubGPF7ir2VjWyr6qRfVVN7Ktq4lB9C/YrZ0gYkBDNjTOznQkqIgFJhR/gKupbWbK5jPcKS1m5pxqP79w3aX2jGdovjpk5KQxKjCU+OoK4qHBio8KJj4pgytBkYqPCHU4vIoFEhR+ASg438f7mQyzeXEbB3mqshWGp8dx6xjDm5A0kt38fnZ5YRI6bWsNB1loaWjuobW6nsqGN5TsqWLKljMIDdQCMGNCH75+VywXj0hkxoI9OZiYiPdKjwjfGzAN+AYwGpllr13S5717gZsADfN9au6Qn2wpW/9hZybJth6hpaqemqY3DTe3UNnd+X+c7F31Xk4Ykcc/5ozhvzECyU+MdSi0ioainR/iFwGXAn7quNMbkAVcBY4BBwAfGmBHWWk8Ptxc02jq8/GbJNp5cXkRsZDgpfaJIioskOS6KzORYkuIiSYzt+hXFpCFJDEiIcTq6iISoHhW+tXYrcKShhouBl621rUCRMWYXMA1Y0ZPtBYvi6ibufGkd64truG7GUP7jgtHEROoNVBFxlr/G8DOAlV2WS3zr/oUxZj4wH2DIkCF+itN7FheW8uPXNwLw+Hcmc/64dIcTiYh0OmbhG2M+AAYe4a6fWmvfPtofO8K6I15Lz1q7AFgAkJ+ff+Tr7QWJPyzbye+W7mBCZiJ/vHoyQ1LinI4kIvJPxyx8a+05J/BzS4DBXZYzgYMn8HOCxlPL9/C7pTu4bFIGD14+nqgInbVCRAKLv1ppEXCVMSbaGJMN5AKr/bQtx71aUMwD727l/LED+fUVKnsRCUw9aiZjzKXGmBJgBvCuMWYJgLV2M/AqsAVYDNweqjN0/r6plHv+upHTclP5/VUTidB1YkUkQPV0ls6bwJtHue+XwC978vMD3Sc7KvjBy+uYNCSZP107hegIzcQRkcClT9oep3aPl4K91Xy4tZznV+0jt39fnrlhKnFR+qsUkcCmluoGay1/21jKks1lfLq9gvrWDqLCwzgtN5UHLx9PYmyk0xFFRI5JhX8M1loeeHcrT/+jiLS+0VwwLp2zRvfn1JxUncBMRIKKGusYHl22k6f/UcQNM7P4+YV5hIXpBGYiEpw0peRrPLV8D7//YCdXTMlU2YtI0FPhH8UrBfv/Obf+wcvGqexFJOi5fkinvL6FstrOSwR6rcUCWw7Wcd/bhZwxIk1z60UkZLiy8Fs7PHywpZxX1xSzfGcF3iOcwWdaVj+euEZz60UkdLim8FvaPWworuG9wjLeWn+AmqZ20hNjuH12DuMzkwgzEGYMGIgMCyM/K1mnNBaRkBKyhd/a4WHlnmpWF1VRUHSY9cU1tHm8REWEcW7eAK7MH8ysnFTCNTYvIi4RcoXf0u7hpdX7+dMneyirayE8zDA2I5EbZmUxLasfU7P76YNSIuJKIVP4ja0dvLBqHws+LaKyoZVpWf144JKxzBieog9IiYgQIoX/4bZD/PDVDRxuamdWTgr/fdYkThmW4nQsEZGAEhKFn5USz6Qhydw+O4cpQ5OdjiMiEpBCovCHpfXhmRumOh1DRCSg6RNFIiIuocIXEXEJFb6IiEuo8EVEXEKFLyLiEip8ERGXUOGLiLiECl9ExCWMtUc4GbxDjDEVwL4T/OOpQOVJjOOEYN8H5XdesO+D8p+YodbatGM9KKAKvyeMMWustflO5+iJYN8H5XdesO+D8vuXhnRERFxChS8i4hKhVPgLnA5wEgT7Pii/84J9H5Tfj0JmDF9ERL5eKB3hi4jI1wiJwjfGzDXGbDfG7DLG3ON0nu4wxjxjjCk3xhR2WdfPGLPUGLPTdxuQV3Mxxgw2xnxkjNlqjNlsjPmBb31Q5AcwxsQYY1YbYzb49uF+3/psY8wq3z68YoyJcjrr1zHGhBtj1hlj3vEtB01+Y8xeY8wmY8x6Y8wa37qgeQ4BGGOSjDGvG2O2+V4PMwJ5H4K+8I0x4cBjwPlAHnC1MSbP2VTd8hww9yvr7gGWWWtzgWW+5UDUAfzQWjsaOAW43fd3Hiz5AVqBs6y1E4CJwFxjzCnAQ8Ajvn04DNzsYMbu+AGwtctysOWfba2d2GUqYzA9hwAeBRZba0cBE+j8twjcfbDWBvUXMANY0mX5XuBep3N1M3sWUNhleTuQ7vs+HdjudMZu7sfbwJwgzh8HfAFMp/NDMxG+9f/nuRVoX0AmnYVyFvAOYIIs/14g9SvrguY5BCQARfjeCw2GfQj6I3wgAyjuslziWxeMBlhrSwF8t/0dznNMxpgsYBKwiiDL7xsOWQ+UA0uB3UCNtbbD95BAfy79HvgJ4PUtpxBc+S3wvjFmrTFmvm9dMD2HhgEVwLO+YbWnjDHxBPA+hELhmyOs09SjXmCM6QO8Adxlra1zOs/xstZ6rLUT6TxSngaMPtLDejdV9xhjLgTKrbVru64+wkMDMr/PLGvtZDqHY283xpzudKDjFAFMBh631k4CGgmk4ZsjCIXCLwEGd1nOBA46lKWnDhlj0gF8t+UO5zkqY0wknWX/grX2r77VQZO/K2ttDfAxne9HJBljInx3BfJzaRZwkTFmL/AyncM6vyd48mOtPei7LQfepPOXbjA9h0qAEmvtKt/y63T+AgjYfQiFwi8Acn2zE6KAq4BFDmc6UYuA633fX0/n2HjAMcYY4Glgq7X2d13uCor8AMaYNGNMku/7WOAcOt9w+wi4wvewgN0Ha+291tpMa20Wnc/5D6213yFI8htj4o0xfb/8HjgXKCSInkPW2jKg2Bgz0rfqbGALgbwPTr+JcJLePLkA2EHnGOxPnc7TzcwvAaVAO51HCjfTOQa7DNjpu+3ndM6jZD+VzqGCjcB639cFwZLftw/jgXW+fSgEfu5bPwxYDewCXgOinc7ajX05E3gnmPL7cm7wfW3+8nUbTM8hX96JwBrf8+gtIDmQ90GftBURcYlQGNIREZFuUOGLiLiECl9ExCVU+CIiLqHCFxFxCRW+iIhLqPBFRFxChS8i4hL/C6nIWp9RkkQ+AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(np.linalg.eigvalsh(H))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.00000000e+00+0.00000000e+00j, -9.88415498e-15-7.49400542e-17j,\n",
       "        -1.11022302e-16+6.38378239e-16j, -8.19939154e-15+4.76210391e-18j],\n",
       "       [-9.88415498e-15+7.49400542e-17j,  1.00000000e+00+0.00000000e+00j,\n",
       "        -8.04951143e-15+4.76210391e-18j, -1.11022302e-16-6.38378239e-16j],\n",
       "       [-1.11022302e-16-6.38378239e-16j, -8.04951143e-15-4.76210391e-18j,\n",
       "         1.00000000e+00+0.00000000e+00j, -9.89367919e-15+7.49400542e-17j],\n",
       "       [-8.19939154e-15-4.76210391e-18j, -1.11022302e-16+6.38378239e-16j,\n",
       "        -9.89367919e-15-7.49400542e-17j,  1.00000000e+00+0.00000000e+00j]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U = np.linalg.matrix_power(M, 2)\n",
    "(U.conj().T.dot(U))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.35355339-3.53553391e-01j,  0.5       +1.53869061e-15j,\n",
       "        -0.35355339+3.53553391e-01j,  0.5       +1.88411095e-16j],\n",
       "       [ 0.5       +1.53869061e-15j,  0.35355339-3.53553391e-01j,\n",
       "         0.5       +1.53869061e-15j, -0.35355339+3.53553391e-01j],\n",
       "       [-0.35355339+3.53553391e-01j,  0.5       +1.53869061e-15j,\n",
       "         0.35355339-3.53553391e-01j,  0.5       +1.53869061e-15j],\n",
       "       [ 0.5       +1.88411095e-16j, -0.35355339+3.53553391e-01j,\n",
       "         0.5       +1.53869061e-15j,  0.35355339-3.53553391e-01j]])"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Harmonic Oscillator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    },
    {
     "ename": "AssertionError",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAssertionError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-8-1846013a711c>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     21\u001b[0m \u001b[0mM\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdx\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexp\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1j\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0mh\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mm\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx_\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0my_\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;36m2.0\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0mdt\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mV\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx_\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0my_\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0mB\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     22\u001b[0m \u001b[0mI\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mM\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mconj\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mT\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mM\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 23\u001b[0;31m \u001b[0;32massert\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mallclose\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0meye\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mNx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mI\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mAssertionError\u001b[0m: "
     ]
    }
   ],
   "source": [
    "%pylab inline --no-import-all\n",
    "\n",
    "Nx = 64\n",
    "L = 10.0\n",
    "w = m = h = 1.0\n",
    "dx = L/Nx\n",
    "x = np.arange(Nx)*dx - L/2\n",
    "\n",
    "n_t = 20\n",
    "Nt = 2*n_t+1\n",
    "dt = m*L*dx/2/np.pi/Nt/h\n",
    "T = Nt*dt\n",
    "\n",
    "def V(x):\n",
    "    return m*(w*x)**2/2.0\n",
    "\n",
    "x_ = x[:,None]\n",
    "y_ = x[None,:]\n",
    "\n",
    "B = np.sqrt(2j*Nt*np.pi*h*dt/m)\n",
    "M = dx*np.exp(1j/h*(m*(x_-y_)**2/2.0/dt - V(x_+y_)))/B\n",
    "I = M.conj().T.dot(M)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:_gpe]",
   "language": "python",
   "name": "conda-env-_gpe-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  },
  "nikola": {
   "category": "",
   "date": "2018-08-23 16:14:34 UTC-07:00",
   "description": "",
   "link": "",
   "slug": "path-integrals",
   "tags": "",
   "title": "Path Integrals",
   "type": "text"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": true,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}