{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ff8465b6",
   "metadata": {},
   "source": [
    "| [**Lattice QCD I**](https://homepages.uni-regensburg.de/~lec17310/teaching/wise2122/lqft.html) | - | [**Christoph Lehner**](https://go.ur.de/lehner) |\n",
    "|:-|--|-:|\n",
    "\n",
    "# Chapter 10: fermions on a lattice\n",
    "\n",
    "Adding fermions to our so-far purely bosonic theories requires several modifications that we will now discuss."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "805d159d",
   "metadata": {},
   "source": [
    "## Fermionic path integral\n",
    "\n",
    "The quantization of bosonic fields $\\hat{q}_a$ and their conjugate momenta fields $\\hat{p}_a$ required the imposition of commutation relations\n",
    "$$\n",
    "\\begin{align}\n",
    "[\\hat{q}_a,\\hat{p}_b] &= i \\delta_{ab} \\,, \\\\\n",
    "[\\hat{q}_a,\\hat{q}_b] &= 0 \\,, \\\\\n",
    "[\\hat{p}_a,\\hat{p}_b] &= 0\n",
    "\\end{align}\n",
    "$$\n",
    "extending trivially the derivation of the path integral for the one-dimensional Hamiltonian systems.\n",
    "\n",
    "The Fermi exclusion principle, however, is incompatible with the above quantization rules.  In particular,\n",
    "we need $\\hat{q}_a^2 = 0$ corresponding to the limitation that fermions occupy each possible state with occupation number 0 or 1.  This can be achieved modifying the above relations to **anticommutator** relations\n",
    "$$\n",
    "\\begin{align}\\label{eqn:full}\n",
    "\\{ \\hat{q}_a,\\hat{p}_b \\} &= i \\delta_{ab} \\,, \\\\\n",
    "\\{\\hat{q}_a,\\hat{q}_b\\} &= 0 \\,, \\\\\n",
    "\\{\\hat{p}_a,\\hat{p}_b\\} &= 0 \\,,\n",
    "\\end{align}\n",
    "$$\n",
    "where $\\{A,B\\} = AB + BA$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d3672560",
   "metadata": {},
   "source": [
    "Let us first consider the case of a single fermionic degree of freedom $\\hat{q}$ which, one could write in matrix form as\n",
    "$$\n",
    "\\begin{align}\n",
    "\\hat{q} = \\begin{pmatrix} 0 & i \\\\ 0 & 0 \\end{pmatrix}\n",
    "\\end{align}\n",
    "$$\n",
    "which would indeed satisfy $\\hat{q}^2 = 0$.  We can complete this variable by a canonical momentum\n",
    "$$\n",
    "\\begin{align}\n",
    "\\hat{p} = \\begin{pmatrix} 0 & 0 \\\\ 1 & 0 \\end{pmatrix}\n",
    "\\end{align}\n",
    "$$\n",
    "for which the above anticommutator relations hold.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac5ac65c",
   "metadata": {},
   "source": [
    "Because of $\\hat{q}^2=0$ and $\\hat{p}^2=0$ we can always consider vectors and dual vectors\n",
    "$$\n",
    "\\begin{align}\n",
    "\\vert 0 \\rangle &= \\hat{q} \\vert \\Omega \\rangle \\,, \\\\\n",
    "\\langle 0 \\vert &= \\langle \\Omega \\vert \\hat{p} \n",
    "\\end{align}\n",
    "$$\n",
    "which satisfy\n",
    "$$\n",
    "\\begin{align}\n",
    "\\hat{q} \\vert 0 \\rangle = \\langle 0 \\vert \\hat{p} = 0\n",
    "\\end{align}\n",
    "$$\n",
    "based on any $\\vert \\Omega \\rangle$ and $\\langle \\Omega \\vert$ for which $\\vert 0 \\rangle$ and $\\langle 0 \\vert$ do not vanish.  In our example, this could be\n",
    "$$\n",
    "\\begin{align}\n",
    "\\vert 0 \\rangle &= \\begin{pmatrix} 1 \\\\ 0 \\end{pmatrix} \\,, \\\\\n",
    "\\langle 0 \\vert &= \\begin{pmatrix} 1 & 0 \\end{pmatrix} \\,.\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00981372",
   "metadata": {},
   "source": [
    "In our derivation of the bosonic path integral, we were next using completeness relations\n",
    "$$\n",
    "\\begin{align}\n",
    "1 = \\int dq \\vert q \\rangle \\langle q \\vert\n",
    "\\end{align}\n",
    "$$\n",
    "for which eigenvectors $\\vert q \\rangle$ of $\\hat{q}$ were needed.  The fermionic operators $\\hat{q}$, however, have no non-zero eigenvalues and no basis of eigenvectors in the usual sense.  We therefore will need to devise an alternative approch.\n",
    "\n",
    "We first note that in our example\n",
    "$$\n",
    "\\begin{align}\n",
    "\\hat{p} \\vert 0 \\rangle = \\begin{pmatrix} 0 \\\\ 1 \\end{pmatrix}\n",
    "\\end{align}\n",
    "$$\n",
    "together with $\\vert 0 \\rangle$ forms a basis of the space in which $\\hat{q}$ and $\\hat{p}$ act.  To complete the discussion along the lines of the bosonic path integral, we follow the nice exposition of chapter 9.5 of Steven Weinberg's book The Quantum Theory of Fields I."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6f6b875",
   "metadata": {},
   "source": [
    "Let us consider the full list of operators $\\hat{q}_a$ and $\\hat{p}_a$ endowed with anticommutation rules of Eq. \\eqref{eqn:full} and an additional algebra of **Grassmann variables** $q_a$ and $p_q$ that satisfy\n",
    "$$\n",
    "\\begin{align}\n",
    "\\{ q_a, q_b \\} = \\{ q_a, p_b \\} = \\{ q_a, \\hat{q}_b \\} = \\{ q_a, \\hat{p}_b \\} =\n",
    "\\{ p_a, p_b \\} = \\{ p_a, \\hat{q}_b \\} = \\{ p_a, \\hat{p}_b \\} = 0 \\,.\n",
    "\\end{align}\n",
    "$$\n",
    "The state\n",
    "$$\n",
    "\\begin{align}\n",
    "\\vert q \\rangle = \\exp\\left( - i \\sum_a \\hat{p}_a q_a \\right) \\vert 0 \\rangle\n",
    "\\end{align}\n",
    "$$\n",
    "then satisfies\n",
    "$$\n",
    "\\begin{align} \\label{eqn:reigen}\n",
    "\\hat{q}_a \\vert q \\rangle = q_a \\vert q \\rangle\n",
    "\\end{align}\n",
    "$$\n",
    "with state $\\vert 0 \\rangle$ for which $\\hat{q}_a \\vert 0 \\rangle = 0$ and $\\langle 0 \\vert$ for which $\\langle 0 \\vert \\hat{p}_a = 0$.  We can then also define a state\n",
    "$$\n",
    "\\begin{align}\n",
    "\\langle q \\vert = \\langle 0 \\vert \\left(\\prod_{a} \\hat{q}_a \\right) \\exp\\left( - i \\sum_a q_a \\hat{p}_a \\right)\n",
    "\\end{align}\n",
    "$$\n",
    "that satisfies\n",
    "$$\n",
    "\\begin{align}\\label{eqn:leigen}\n",
    "\\langle q \\vert \\hat{q}_a = \\langle q \\vert q_a \\,.\n",
    "\\end{align}\n",
    "$$\n",
    "Since the order of factors is now important, we define the canonical and reverse canonical product explicitly as\n",
    "$$\n",
    "\\begin{align}\n",
    "\\prod_a \\hat{q}_a &= \\hat{q}_1 \\hat{q}_2 \\cdots \\hat{q}_n \\,, \\\\\n",
    "\\tilde{\\prod}_a \\hat{q}_a &= \\hat{q}_n \\hat{q}_{n-1} \\cdots \\hat{q}_1 \\,.\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "***Homework***: Show that equations \\eqref{eqn:reigen} and \\eqref{eqn:leigen} hold."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58621ea2",
   "metadata": {},
   "source": [
    "At this point we have again a concept similar to eigenvectors and eigenvalues for the bosonic degrees of freedom.  We lack, however, a completeness relation for the corresponding states.\n",
    "\n",
    "One first can show that\n",
    "$$\n",
    "\\begin{align}\n",
    "\\langle q' \\vert q \\rangle = \\prod_a (q_a - q'_a)\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "***Homework:*** Show this equation.\n",
    "\n",
    "Before completing the discussion of completeness relations, we shall see that this product can play the role of a delta function similar to the bosonic case."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03a59f09",
   "metadata": {},
   "source": [
    "We can similarly show that\n",
    "$$\n",
    "\\begin{align}\n",
    "\\hat{p}_a \\vert p \\rangle &= p_a \\vert p \\rangle \\,, \\\\\n",
    "\\langle p \\vert \\hat{p}_a &= \\langle p \\vert p_a \\,, \\\\\n",
    "\\vert p \\rangle &= \\exp \\left(-i\\sum_a \\hat{q}_a p_a \\right) \\left(\\prod_a \\hat{p}_a \\right) \\vert 0 \\rangle \\,, \\\\\n",
    "\\langle p \\vert &= \\langle 0 \\vert \\exp \\left(-i\\sum_a p_a \\hat{q}_a \\right) \\,, \\\\\n",
    "\\langle p' \\vert p \\rangle &= \\prod_a (p'_a - p_a) \\,.\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cefa526e",
   "metadata": {},
   "source": [
    "Next, we compute the inner product between two such states\n",
    "$$\n",
    "\\begin{align}\n",
    "\\langle q \\vert p \\rangle &=\n",
    "\\langle q \\vert \\exp \\left(-i\\sum_a \\hat{q}_a p_a \\right) \\left(\\prod_a \\hat{p}_a \\right) \\vert 0 \\rangle \\\\\n",
    "&=\\langle q \\vert \\exp \\left(-i\\sum_a q_a p_a \\right) \\left(\\prod_a \\hat{p}_a \\right) \\vert 0 \\rangle \\\\\n",
    "&=\\exp \\left(-i\\sum_a q_a p_a \\right) \\langle q \\vert  \\left(\\prod_a \\hat{p}_a \\right) \\vert 0 \\rangle \\\\\n",
    "&=\\exp \\left(-i\\sum_a q_a p_a \\right) \\langle 0 \\vert \\left(\\prod_{a} \\hat{q}_a \\right) \\exp\\left( - i \\sum_a q_a \\hat{p}_a \\right)  \\left(\\prod_a \\hat{p}_a \\right) \\vert 0 \\rangle \\\\\n",
    "&=\\exp \\left(-i\\sum_a q_a p_a \\right) \\langle 0 \\vert \\left(\\prod_{a} \\hat{q}_a \\right) \\left(\\prod_{a}\\exp\\left( - i q_a \\hat{p}_a \\right) \\hat{p}_a \\right) \\vert 0 \\rangle \\\\\n",
    "&=\\exp \\left(-i\\sum_a q_a p_a \\right) \\langle 0 \\vert \\left(\\prod_{a} \\hat{q}_a \\right) \\left(\\prod_{a}\\left(\\hat{p}_a - i q_a \\hat{p}_a^2 \\right) \\right) \\vert 0 \\rangle \\\\\n",
    "&=\\exp \\left(-i\\sum_a q_a p_a \\right) \\langle 0 \\vert \\left(\\prod_{b} \\hat{q}_b \\right)  \\left(\\prod_a \\hat{p}_a \\right) \\vert 0 \\rangle \\\\\n",
    "&=\\chi_n \\exp \\left(-i\\sum_a q_a p_a \\right) =\\chi_n \\exp \\left(i\\sum_a p_a q_a \\right)\n",
    "\\end{align}\n",
    "$$\n",
    "with\n",
    "$$\n",
    "\\begin{align}\n",
    "\\chi_n = i^n (-1)^{n(n-1)/2}\n",
    "\\end{align}\n",
    "$$\n",
    "since $\\hat{q}_a p_a$ commute with each other and where $n$ is the number of operators $\\hat{q}_a$ and $\\hat{p}_a$, respectively.\n",
    "Similarly, we find\n",
    "$$\n",
    "\\begin{align}\n",
    "\\langle p \\vert q \\rangle &= \\exp \\left(-i\\sum_a p_a q_a \\right) \\,.\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc04663a",
   "metadata": {},
   "source": [
    "Extending the $n=1$ example given at the beginning of this section, a basis of the space within which the $\\hat{q}_a$ and $\\hat{p}_a$ act is given by\n",
    "$$\n",
    "\\begin{align}\n",
    "\\vert a,b,\\ldots \\rangle &= \\hat{p}_a \\hat{p}_b \\cdots \\vert 0 \\rangle \\,, \\\\\n",
    "\\langle a,b,\\ldots \\vert &= \\langle 0 \\vert \\cdots (-i\\hat{q}_b)(-i\\hat{q}_a)\n",
    "\\end{align}\n",
    "$$\n",
    "with\n",
    "$$\n",
    "\\begin{align}\n",
    "\\hat{q}_a \\vert b,c,\\ldots \\rangle &= 0 \\,, \\\\\n",
    "\\hat{p}_a \\vert b,c,\\ldots \\rangle &= \\vert a,b,c,\\ldots \\rangle \\,, \\\\\n",
    "\\hat{q}_a \\vert a,b,c,\\ldots \\rangle &= i \\vert b,c,\\ldots \\rangle \\,, \\\\\n",
    "\\hat{p}_a \\vert a,b,c,\\ldots \\rangle &= 0 \\,, \\\\\n",
    "\\langle c,d,\\ldots \\vert a,b,\\ldots \\rangle &= \\begin{cases} 0 & \\textrm{if } \\{ c,d,\\ldots \\} \\neq \\{a,b,\\ldots \\} \\\\\n",
    "1 & \\textrm{if } c=a, d=b, \\ldots\n",
    "\\end{cases} \\,.\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fef1aaab",
   "metadata": {},
   "source": [
    "If we now consider again the state\n",
    "$$\n",
    "\\begin{align}\n",
    "\\vert q \\rangle &= \\exp\\left( - i \\sum_a \\hat{p}_a q_a \\right) \\vert 0 \\rangle \\\\\n",
    "&=\\vert 0 \\rangle + i \\sum_a q_a \\vert a \\rangle - \\frac12 \\sum_{a,b} q_a q_b \\vert a,b \\rangle + \\ldots\n",
    "\\end{align}\n",
    "$$\n",
    "we note that we can expand any state $\\vert f \\rangle$ in this basis by considering\n",
    "the coefficients of $q_a q_b \\cdots$ in $\\vert q \\rangle$ which we denote by $\\vert q \\rangle_{ab\\ldots}$.\n",
    "We therefore can write\n",
    "$$\n",
    "\\begin{align}\n",
    "\\vert f \\rangle &= f_0 \\vert q \\rangle_0 + \\sum_a f_a \\vert q \\rangle_{a} + \\sum_{a \\neq b} f_{ab} \\vert q \\rangle_{ab} + \\ldots\n",
    "\\end{align}\n",
    "$$\n",
    "In order to introduce a sum over states it will therefore be very useful to introduce an operation that projects out the coefficients of $q_a q_b \\cdots$ from such expressions.  This can be done using the **Berezin integration** over Grassmann variables.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd378a71",
   "metadata": {},
   "source": [
    "Any function $g(\\xi)$ of a Grassmann variable $\\xi$ can be written as\n",
    "$$\n",
    "\\begin{align}\n",
    "g(\\xi) = g_0 +  \\xi g_1\\,,\n",
    "\\end{align}\n",
    "$$\n",
    "where $g_1$ is a Grassmann variable as well if $g(\\xi)$ is a commuting object.  This must hold since the algebra has now higher powers of $\\xi$.  The order of $\\xi$ and $g_1$ is therefore important.  The Berezin integration is then defined as the operation\n",
    "$$\n",
    "\\begin{align}\n",
    "\\int d\\xi g(\\xi) = g_1 \\,.\n",
    "\\end{align}\n",
    "$$\n",
    "Consider now the case of two Grassmann variables $\\xi_1,\\xi_2$ and\n",
    "$$\n",
    "\\begin{align}\n",
    "g(\\xi) = g_0 +  \\xi_1 g_1 + \\xi_2 g_2 + \\xi_1 \\xi_2 g_3 \\,,\n",
    "\\end{align}\n",
    "$$\n",
    "Then the iterated Berezin integrals\n",
    "$$\n",
    "\\begin{align}\n",
    "\\int d\\xi_2 \\int d\\xi_1 g(\\xi) = \\int d\\xi_2 (g_1 + \\xi_2 g_3) = g_3\n",
    "\\end{align}\n",
    "$$\n",
    "project out the coefficient of $\\xi_1 \\xi_2$.  Note, however, that the order of integration was reversed.\n",
    "One can consistently define the infinitesimals $d\\xi_i$ to be Grassmann variables themself, which anticommute\n",
    "with all other Grassmann variables and themself."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c5bd38e",
   "metadata": {},
   "source": [
    "With this projection at hand, we can therefore write any state as\n",
    "$$\n",
    "\\begin{align}\n",
    "\\vert f \\rangle = \\int \\tilde{\\prod}_a dq_a f(q) \\vert q \\rangle\n",
    "\\end{align}\n",
    "$$\n",
    "for a given function $f(q)$ of all $q_a$.  After a lengthy but straightforward calculation, one finds\n",
    "$$\n",
    "\\begin{align}\n",
    "\\langle q' \\vert f \\rangle &= (-1)^n f(q')\n",
    "\\end{align}\n",
    "$$\n",
    "such that for any state\n",
    "$$\n",
    "\\begin{align}\n",
    "\\vert f \\rangle = (-1)^n \\int \\vert q \\rangle \\left(\\tilde{\\prod}_b dq_b \\right) \\langle q \\vert f \\rangle\n",
    "\\end{align}\n",
    "$$\n",
    "or our desired completeness relation\n",
    "$$\n",
    "\\begin{align}\n",
    "1 = \\int \\vert q \\rangle \\left(\\tilde{\\prod}_b (-dq_b) \\right) \\langle q \\vert  \\,.\n",
    "\\end{align}\n",
    "$$\n",
    "In the same way one can also show\n",
    "$$\n",
    "\\begin{align}\n",
    "1 = \\int \\vert p \\rangle \\left(\\tilde{\\prod}_b dp_b \\right) \\langle p \\vert  \\,.\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e4455dd7",
   "metadata": {},
   "source": [
    "With these relations we can now finally complete the derivation of the fermionic path integral.  The key ingredient is again to express the transfer matrix in a path integral way.  Let us consider without loss of generality a Hamiltonian\n",
    "$H(\\hat{p},\\hat{q})$, where all $\\hat{q}_a$ stand to the right of all $\\hat{p}_a$.\n",
    "\n",
    "We then find\n",
    "$$\n",
    "\\begin{align}\n",
    "\\langle q' \\vert e^{-i\\delta t H(\\hat{p},\\hat{q})} \\vert q \\rangle &= \\int \n",
    "\\langle q' \\vert p \\rangle  \\left(\\tilde{\\prod}_b dp_b \\right) \\langle p \\vert e^{-i\\delta t H(\\hat{p},\\hat{q})} \\vert q \\rangle \\\\\n",
    "&= \\int \n",
    "\\langle q' \\vert p \\rangle  \\left(\\tilde{\\prod}_b dp_b \\right) \\langle p \\vert q \\rangle e^{-i\\delta t H(p,q)}  \\\\\n",
    "&= \\chi_n \\int \n",
    " \\left(\\tilde{\\prod}_b dp_b \\right) \\exp\\left(i\\sum_a p_a(q'_a - q_a) -i\\delta t H(p,q) \\right)  \\,.\n",
    "\\end{align}\n",
    "$$\n",
    "We can write this as an operator equation by inserting further completeness relations\n",
    "$$\n",
    "\\begin{align}\n",
    "e^{-i\\delta t H(\\hat{p},\\hat{q})} &= \n",
    "\\chi_n \\int \\vert q' \\rangle \\left(\\tilde{\\prod}_a (-dq'_a) \\right) \\langle q' \\vert\n",
    "e^{-i\\delta t H(\\hat{p},\\hat{q})}\n",
    " \\vert q \\rangle \\left(\\tilde{\\prod}_b (-dq_b) \\right) \\langle q \\vert \\\\\n",
    " &=\\chi_n  \\int \\vert q' \\rangle \\left(\\tilde{\\prod}_a dq'_a \\right) \\left(\\tilde{\\prod}_b dp_b \\right)\n",
    " \\exp\\left(i\\sum_a p_a(q'_a - q_a) -i\\delta t H(p,q) \\right)\n",
    " \\left(\\tilde{\\prod}_a dq_a \\right) \\langle q \\vert \\,.\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "We again note the Legendre transformation that creates now a Berezin integral over the action.  The transition to Euclidean time $\\delta\\tau = i \\delta t$ goes through in the same way as in the bosonic case."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2bafb67a",
   "metadata": {},
   "source": [
    "## Euclidean quark action for QCD\n",
    "\n",
    "We remind ourselves of the real-time action\n",
    "$$\n",
    "\\begin{align}\n",
    "  {\\cal S} &= \\int d^4x \\bar\\psi (i \\gamma^\\mu \\partial_\\mu - m) \\psi\n",
    "\\end{align}\n",
    "$$\n",
    "with Minkowski gamma matrices in our chiral representation of chapter 5\n",
    "$$\n",
    "\\begin{align}\\label{eqn:defgamma}\n",
    "  \\gamma^0 = \\begin{pmatrix}\n",
    "    0 & \\mathbb{1} \\\\ \\mathbb{1} & 0\n",
    "  \\end{pmatrix},\\qquad \\gamma^i = \\begin{pmatrix} 0 & \\sigma_i \\\\\n",
    "    -\\sigma_i & 0\n",
    "  \\end{pmatrix} \\,.\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "The conjugate momentum to the field $\\psi$ is therefore given by\n",
    "$$\n",
    "\\begin{align}\n",
    "\\pi = \\bar\\psi i \\gamma^0 \\,.\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "We now introduce a Euclidean space formulation by\n",
    "considering the four vector in Euclidean space $\\tilde{x}$ that\n",
    "relates to the Minkowski vector $x$ by\n",
    "$$\n",
    "\\begin{align}\n",
    "(\\tilde{x}^1,\\tilde{x}^2,\\tilde{x}^3,\\tilde{x}^4)=(-x^1,-x^2,-x^3,i x^0)\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "We therefore have\n",
    "$$\n",
    "\\begin{align}\n",
    "\\partial_0 &= \\frac{\\partial}{\\partial x^0} = \\frac{\\partial \\tilde{x}^4}{\\partial x^0}\\frac{\\partial}{\\partial \\tilde{x}^4}\n",
    "=i \\tilde{\\partial}_4 \\,, &\n",
    "\\partial_i &= -\\tilde{\\partial}_i \\,,\n",
    "\\end{align}\n",
    "$$\n",
    "and\n",
    "$$\n",
    "\\begin{align}\n",
    " dx^0 &= -i d\\tilde{x}^4 \\,, &\n",
    " dx^1 dx^2 dx^3 &= \\vert (-1)^3 \\vert d\\tilde{x}^1 d\\tilde{x}^2 d\\tilde{x}^3 = d\\tilde{x}^1 d\\tilde{x}^2 d\\tilde{x}^3\n",
    "\\end{align}\n",
    "$$\n",
    "and find the Euclidean action\n",
    "$$\n",
    "\\begin{align}\n",
    "  {\\cal S}_E &= -i {\\cal S} \\\\\n",
    "  &= -\\int d\\tilde{x}^1 d\\tilde{x}^2 d\\tilde{x}^3 d\\tilde{x}^4 \\bar\\psi (- \\gamma_0 \\tilde{\\partial}_4 - i \\gamma^i \\tilde{\\partial}_i - m) \\psi \\\\\n",
    "  &= \\int d\\tilde{x}^1 d\\tilde{x}^2 d\\tilde{x}^3 d\\tilde{x}^4 \\bar\\psi (\\gamma^0 \\tilde{\\partial}_4 + i \\gamma^i \\tilde{\\partial}_i + m) \\psi \\\\\n",
    "  &= \\int d\\tilde{x}^1 d\\tilde{x}^2 d\\tilde{x}^3 d\\tilde{x}^4 \\bar\\psi (\\tilde{\\gamma}^4 \\tilde{\\partial}_4 + \\tilde{\\gamma}^i \\tilde{\\partial}_i + m) \\psi \\\\\n",
    "\\end{align}\n",
    "$$\n",
    "with\n",
    "$$\n",
    "\\begin{align}\n",
    "\\tilde{\\gamma}^4 &= \\gamma^0 \\,, &\n",
    "\\tilde{\\gamma}^i &= i \\gamma^i \\,.\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "b7ed52c1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import gpt as g"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "4bf16a29",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "gamma[1] =\n",
      " [[ 0.+0.j  0.+0.j  0.+0.j  0.+1.j]\n",
      " [ 0.+0.j  0.+0.j  0.+1.j  0.+0.j]\n",
      " [ 0.+0.j -0.-1.j  0.+0.j  0.+0.j]\n",
      " [-0.-1.j  0.+0.j  0.+0.j  0.+0.j]]\n",
      "gamma[2] =\n",
      " [[ 0.+0.j  0.+0.j  0.+0.j -1.+0.j]\n",
      " [ 0.+0.j  0.+0.j  1.+0.j  0.+0.j]\n",
      " [ 0.+0.j  1.+0.j  0.+0.j  0.+0.j]\n",
      " [-1.+0.j  0.+0.j  0.+0.j  0.+0.j]]\n",
      "gamma[3] =\n",
      " [[ 0.+0.j  0.+0.j  0.+1.j  0.+0.j]\n",
      " [ 0.+0.j  0.+0.j  0.+0.j -0.-1.j]\n",
      " [-0.-1.j  0.+0.j  0.+0.j  0.+0.j]\n",
      " [ 0.+0.j  0.+1.j  0.+0.j  0.+0.j]]\n",
      "gamma[4] =\n",
      " [[0.+0.j 0.+0.j 1.+0.j 0.+0.j]\n",
      " [0.+0.j 0.+0.j 0.+0.j 1.+0.j]\n",
      " [1.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
      " [0.+0.j 1.+0.j 0.+0.j 0.+0.j]]\n"
     ]
    }
   ],
   "source": [
    "print(\"gamma[1] =\\n\",g.gamma[0].tensor().array)\n",
    "print(\"gamma[2] =\\n\",g.gamma[1].tensor().array)\n",
    "print(\"gamma[3] =\\n\",g.gamma[2].tensor().array)\n",
    "print(\"gamma[4] =\\n\",g.gamma[3].tensor().array)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80805256",
   "metadata": {},
   "source": [
    "## Fermionic path integral is a determinant\n",
    "\n",
    "For a fermionic action of the kind\n",
    "$$\n",
    "\\begin{align}\n",
    "S_E = \\sum_{x,y,i,j} \\bar\\psi(x)_i D_{i,x;j,y}(U) \\psi(y)_j\n",
    "\\end{align}\n",
    "$$\n",
    "we find that the path integral formulation previously derived requires us to integrate\n",
    "over the fields $\\psi$ and the corresponding conjugate momenta $\\bar\\psi$ (see above) such that\n",
    "the fermionic part of the path integral yields\n",
    "$$\n",
    "\\begin{align}\n",
    "\\int d\\bar\\psi d\\psi e^{-S_E} = \\det(D(U)) \\,.\n",
    "\\end{align}\n",
    "$$\n",
    "This follows up to a sign from the definition of the\n",
    "Berezin integral.  We define the above integral such that\n",
    "the order of $d\\psi$ and $d\\bar\\psi$ is chosen such that this equation\n",
    "holds without additional sign."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e8eef5e",
   "metadata": {},
   "source": [
    "## Naive lattice fermion\n",
    "\n",
    "A naive lattice discretization of a fermion action would therefore correspond to\n",
    "the Euclidean action\n",
    "$$\n",
    "\\begin{align}\n",
    "S_E = \\sum_x \\bar\\psi(x) (\\gamma_\\mu D_\\mu + m) \\psi(x)\n",
    "\\end{align}\n",
    "$$\n",
    "with covariant derivative from chapter 6\n",
    "$$\n",
    "\\begin{align}\n",
    "D_\\mu \\psi(x) = \\frac1{2a}(C^+_\\mu - C^-_\\mu)\\psi(x) = \\frac1{2a}\\left(U_\\mu(x) \\psi(x + a \\hat{\\mu}) - U^\\dagger_\\mu(x - a\\hat{\\mu}) \\psi(x - a \\hat{\\mu})\\right)\\,.\n",
    "\\end{align}\n",
    "$$\n",
    "We will see in the next lecture that this action actually corresponds to more than\n",
    "one Fermion at the same time, which requires special attention."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "519ee56b",
   "metadata": {},
   "source": [
    "## Fermion doublers\n",
    "\n",
    "The study of the non-interacting Hamiltonian of a fermionic system on a lattice again proceeds similar to chapter 4.  The classical solutions of the non-interacting theory provide structure to the energy eigenstates of the system.\n",
    "Consider in Minkowski momentum space\n",
    "$$\n",
    "\\begin{align}\n",
    "\\psi(x) = e^{i x_\\mu p_\\mu} \\psi_0\n",
    "\\end{align}\n",
    "$$\n",
    "for which the equations of motion yield\n",
    "$$\n",
    "\\begin{align}\n",
    "0 &= \\frac12\\gamma_\\mu (\\psi(x + \\hat{\\mu}) - \\psi(x - \\hat{\\mu}))  + m \\psi(x) \\\\\n",
    "  &= \\left(\\frac12\\gamma_\\mu (e^{i p_\\mu} - e^{-i p_\\mu})  + m \\right) \\psi(x) \\\\\n",
    "  &= \\left(i\\gamma_\\mu \\sin(p_\\mu) + m \\right) \\psi(x) \\,.\n",
    "\\end{align}\n",
    "$$\n",
    "Note that the lattice momentum for such a naive fermion $\\sin(p_\\mu)$ behaves\n",
    "quite differently from the lattice momentum for a boson $2\\sin(p_\\mu/2)$, see chapter 4."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "9a48714f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAiO0lEQVR4nO3deXhV1b3/8feXQJgnmSEEkEFAmWMQh2utQxFRVETFCRyK2mvr7a232lq11fYWtXUeKaIiKloVoRRRQBQn5nkIJMwJQ5gCJEBCku/vjxy9+WmQ4Zxkn5PzeT1Pnpy998pZ3/Nw8slm73XWMndHREQqvypBFyAiIhVDgS8iEicU+CIicUKBLyISJxT4IiJxomrQBfyYxo0be9u2bYMuQ0QkZixYsGCnuzcp61hUB37btm2ZP39+0GWIiMQMM9t4pGO6pCMiEicU+CIicUKBLyISJxT4IiJxQoEvIhInwg58M2ttZjPNbKWZrTCzu8toY2b2jJllmNlSM+sdbr8iInJ8IjEssxD4jbsvNLO6wAIzm+buK0u1uRjoGPrqC7wY+i4iIhUk7MB3963A1tDj/Wa2CmgFlA78QcBYL5mLebaZNTCzFqGfFRGJO8XFzs7cfLbvyyfnYAF7Dx7+7quKGXec2z7ifUb0g1dm1hboBcz53qFWwOZS25mhfT8IfDMbAYwASE5OjmR5IiIVqqjY2bz7AGu27yc9O5e1O3LZknOQLTmH2Lr3IIeLyl6PpEnd6tEd+GZWB3gf+C9333eiz+Puo4BRACkpKVqdRURiQlGxs3ZHLos27WHx5hyWZu4lIzuX/MLi79o0r1eDpIY16dm6AQO6taBVgxo0rVeDBjWr0aBWIvVrVqN+zWrUqFY+42kiEvhmVo2SsH/T3T8oo0kW0LrUdlJon4hITCoqdpZl7eWrjJ18vXYnSzbvJTe/EID6NavRPak+N/VrQ8emdenYrA4dmtahbo1qgdYcduCbmQGvAKvc/YkjNJsE3GVm4ym5WbtX1+9FJNZsyTnIjLRsvkzfwTdrd7HvUEnAd2lRjyt6taJn6wb0Sm5Au8a1KYnG6BKJM/yzgBuBZWa2OLTv90AygLu/BEwBBgAZwAHg5gj0KyJSrtyd9OxcPlmxjU9Wbmdp5l4AkhrW5JLuLTizfWPObN+IRnWqB1zpsYnEKJ0vgR/9UxYanfOf4fYlIlIRNuzM48PFWUxcvIX1O/MA6JXcgHv7d+aiU5vRvkmdgCs8MVE9PbKISEXZk1fA5KVbmLAoi4WbcjCDM9o14taz23FR12Y0rVcj6BLDpsAXkbjl7szbsIc352zko2XbKCgq5pRmdbnv4s5c1qMlLRvUDLrEiFLgi0jc2XfoMBMWZvHmnI2s2Z5L3epVua5vMkNSkujaol5U3nCNBAW+iMSNzD0HGPPlBt6Zt4m8giK6J9XnscHdGdijBbUSK38cVv5XKCJxb3nWXkbNWse/l23FgEt7tOTms9rSPalB0KVVKAW+iFRac9bt4plP0/kqYxd1qlfl1rPbMfzMtpXu2vyxUuCLSKUze90unp6ezjfrdtG4TnV+P6Az16YmUy/gT7oGTYEvIpXG3PW7eWLaamav202TutV5YGBXrktNpmZiQtClRQUFvojEvNXb9vPY1DRmpGXTtG51HhzYlev6JlOjmoK+NAW+iMSsLTkHeWLaGt5fmEmd6lW5t39nhp/ZVmf0R6DAF5GYk5dfyHMzM3jly/XgcNvZ7fjFTzrQsHZi0KVFNQW+iMSM4mLnw8VZjPwojez9+VzRqxW/uagTSQ1rBV1aTFDgi0hMWLI5hz/+awWLNuXQI6k+L93Yh97JDYMuK6Yo8EUkquUcKODRqWm8PXczjetU57GrunNV7ySqVKmc0x+UJwW+iEQld2fi4i08MnklOQcPc9vZ7bj7go6BrxoVyxT4IhJ11u/M44EPl/Nlxk56tm7AG1d0o2vLekGXFfMitabtGGAgkO3up5Vx/CfARGB9aNcH7v5wJPoWkcqjsKiYl2et4+kZ6VRPqMIjg07lur5tSNDlm4iI1Bn+a8BzwNgfafOFuw+MUH8iUsms3raf/3lvCUsz9zKgW3P+eOmplWLRkWgSkcB391lm1jYSzyUi8eVwUTEvfraWZz9Np16Najx/XW8u6d4i6LIqpYq8ht/PzJYAW4B73H1FWY3MbAQwAiA5ObkCyxORipa2bR+/eXcJK7bs49IeLfnjpV1jZkHwWFRRgb8QaOPuuWY2APgQ6FhWQ3cfBYwCSElJ8QqqT0QqUHGxM+ar9Tw2dTX1alblpRv60P+05kGXVelVSOC7+75Sj6eY2Qtm1tjdd1ZE/yISPbbkHOQ37y7hm3W7uLBrM0Ze2U1n9RWkQgLfzJoD293dzSwVqALsqoi+RSR6TFycxR8+XE5RsfPo4G5cndK60q4fG40iNSzzbeAnQGMzywQeAqoBuPtLwFXAnWZWCBwErnV3Xa4RiRO5+YU8+OFyPliURe/kBjx5TU/aNKoddFlxJ1KjdIYe5fhzlAzbFJE4szxrL798exEbd+Vx9/kd+eVPO1A1oUrQZcUlfdJWRMqFu/PqVxsY+VEaJ9VO5K2fn8EZJzcKuqy4psAXkYjLOVDAPf9cwvRV2ZzfuSmPD+nBSZqrPnAKfBGJqKWZOdw5biHZ+w/xwMCu3HJWW92YjRIKfBGJCHdn3JxNPPKvlTSuk8i7t/ejl+arjyoKfBEJ24GCQn7/wTI+XLyFczs14alremq5wSikwBeRsKzfmcftb8wnPTuX/76wE3ed10GLk0QpBb6InLCZadn8avwiqlYxxt6SyjkdmwRdkvwIBb6IHDd35/mZGfx92hq6NK/Hyzf2ofVJWkg82inwReS45OYXcs+7S5i6YhuDerZk5JXdqZmYEHRZcgwU+CJyzDbuyuO21+ezbmcef7ikC7ee3U5DLmOIAl9Ejsk3a3dx55sLABh7SypndWgccEVyvBT4InJUb83ZxIMTl9O2cW1eGZaiic9ilAJfRI6osKiYP/97Fa99vYHzTmnC00N7Ua9GtaDLkhOkwBeRMu09eJi73lrIF+k7+fk57bjv4i4kaHx9TFPgi8gPbN59gFtem8f6nXk8Nrg7V5/eOuiSJAIiMim1mY0xs2wzW36E42Zmz5hZhpktNbPekehXRCJv0aY9XPHCV2zfd4ixt6Yq7CuRSK1C8BrQ/0eOX0zJouUdgRHAixHqV0QiaMqyrVw7aja1EqvywS/O4sz2GolTmUQk8N19FrD7R5oMAsZ6idlAAzNrEYm+RSR87s6Ln63lF28u5LRW9ZnwizPp0LRO0GVJhFXUNfxWwOZS25mhfVsrqH8ROYKiYuehScsZN3sTA7u34G9DelCjmj45WxlF3U1bMxtByWUfkpOTA65GpHI7WFDEL99exPRV27n93JO592edNdNlJVZRgZ8FlL7zkxTa9wPuPgoYBZCSkuLlX5pIfNqZm8+tr89naWYODw86lZv6tQ26JClnFbV0/CTgptBonTOAve6uyzkiAdmwM4/BL35N2tZ9vHRDH4V9nIjIGb6ZvQ38BGhsZpnAQ0A1AHd/CZgCDAAygAPAzZHoV0SO39LMHG5+dR7F7rz18zPo00bLEMaLiAS+uw89ynEH/jMSfYnIifsifQe3v7GAhrUSGXtrKu2baCROPIm6m7YiUj4mLs7inn8uoX2TOrx+SyrN6tUIuiSpYAp8kTjwypfreWTySlLbncQ/bkqhfk1NgBaPFPgilZi789jHq3nxs7X0P7U5T13bU2Ps45gCX6SSKip2/vDhMt6eu5nr+ibzyKDTNNtlnFPgi1RC+YVF/PqdxUxZto27zuvAby7qpKUIRYEvUtnk5Rdyx7gFfJG+kz9c0oXbzjk56JIkSijwRSqRnAMFDH91Hkszc3j8qu4MSdHUxvJ/FPgilUT2vkPc+Mpc1u/K48Ub+vCzU5sHXZJEGQW+SCWQuecAN4yeQ/b+fF4bfjpndtA89vJDCnyRGLd2Ry43jJ5DXn4h427rS+9kTZUgZVPgi8SwFVv2ctMrczGD8SP60bVlvaBLkiimwBeJUQs37WH4mLnUrl6Vcbf11bw4clQKfJEYNHvdLm59bR6N61bnzdv6ktSwVtAlSQxQ4IvEmM/X7GDE2Pkkn1SLN2/rS1NNgibHSIEvEkM+WbGNu95aRIemdXjj1lQa1akedEkSQxT4IjHiX0u28Ot3FnNqq/qMvTmV+rU046Ucn4gscWhm/c1stZllmNl9ZRwfbmY7zGxx6Ou2SPQrEi/eX5DJ3eMX0Tu5IeNuVdjLiQn7DN/MEoDngQuBTGCemU1y95Xfa/qOu98Vbn8i8Wb83E38bsIyzmrfmH/clELNRE1vLCcmEmf4qUCGu69z9wJgPDAoAs8rEvfe+GYD932wjHM7NWH0MIW9hCcSgd8K2FxqOzO07/sGm9lSM3vPzI44o5OZjTCz+WY2f8eOHREoTyQ2vfLleh6YuIILujTj5Rv7aOESCVtEruEfg38Bbd29OzANeP1IDd19lLunuHtKkyZNKqg8kejy0udreWTySi4+rTkvXN+b6lUV9hK+SAR+FlD6jD0ptO877r7L3fNDm6OBPhHoV6RSeu7TdEZ+lMalPVry7NBeJFatqPMyqewi8U6aB3Q0s3ZmlghcC0wq3cDMWpTavAxYFYF+RSqdp6av4W+frOHKXq146pqeVE1Q2EvkhD1Kx90Lzewu4GMgARjj7ivM7GFgvrtPAn5lZpcBhcBuYHi4/YpUJu7Ok9PW8MynGVzVJ4lHB3fX+rMScebuQddwRCkpKT5//vygyxApV+7O4x+v5oXP1nJNSmv+emU3qijs5QSZ2QJ3TynrmD5pKxIgd2fk1DRe/nwdQ1OT+cvlpynspdwo8EUC4u7875RV/OOL9dxwRjIPX6awl/KlwBcJgLvzl3+vYvSX67mpXxv+dNmpmCnspXxpCIBIBXN3/hwK+2EKe6lAOsMXqUDuziOTVzHmq/UMP7MtD13aVWEvFUaBL1JB3J2HJ6/k1a82cPNZbXlwoMJeKpYCX6QClA77W85qxwMDuyjspcLpGr5IOVPYS7RQ4IuUI4W9RBMFvkg5+fYG7bfX7BX2EjQFvkg5+Hbo5Ziv1usGrUQNBb5IhH0b9q98WTL0UmEv0UKBLxJB336C9tuw1zh7iSYKfJEIcXf++lHad5+gVdhLtFHgi0SAuzPyozRGzVrHjWe04Y+aLkGikAJfJEzuzqNTV/PyrHUls14OUthLdIpI4JtZfzNbbWYZZnZfGcerm9k7oeNzzKxtJPoVCZq789jHq3np87Vc17dkimOFvUSrsAPfzBKA54GLga7AUDPr+r1mtwJ73L0D8CTwaLj9igTN3fnbJ6t58bO1DE1N5s+DNJ+9RLdInOGnAhnuvs7dC4DxwKDvtRkEvB56/B5wvuk0SGKYu/PEtDU8P3Mt157eWitVSUyIROC3AjaX2s4M7SuzjbsXAnuBRmU9mZmNMLP5ZjZ/x44dEShPJPKenJ7Os59mcE1Ka/73Cq1BK7Eh6m7auvsod09x95QmTZoEXY7IDzw1fQ3PzEjn6pQkLTguMSUSgZ8FtC61nRTaV2YbM6sK1Ad2RaBvkQr19PR0npqezlV9khh5ZXeFvcSUSAT+PKCjmbUzs0TgWmDS99pMAoaFHl8FfOruHoG+RSrMszPSeXL6Ggb3TuLRwQp7iT1hL4Di7oVmdhfwMZAAjHH3FWb2MDDf3ScBrwBvmFkGsJuSPwoiMeO5T9P5+7Q1XNmrFY9d1Z0Ehb3EoIiseOXuU4Ap39v3YKnHh4AhkehLpKI9PzODv32yhit6teLxIT0U9hKzou6mrUg0efGztTz+8WoG9WzJ3xT2EuMU+CJH8OJna3l0ahqX9WjJ3xX2Ugko8EXK8NLn/xf2T1zdg6oJ+lWR2Kd3scj3vPT5WkZ+lMalCnupZPROFimldNg/qbCXSkbvZpGQl0NhP7B7C4W9VEoRGZYpEuu+vUE7sHsLnrqmp8JeKiUFvsS9Fz7L4LGpq3UZRyo9vbMlrn0b9pcp7CUO6Axf4tbzMzO++1DV34co7KXy0ztc4pLCXuKRzvAl7jw9vWTWy8t7tuTvV/fUJ2glbijwJW64O09NT+fpGekM7p2kWS8l7ijwJS58uwbts59mMKRPEiMHK+wl/ijwpdJzdx77eDUvflay4LjWoJV4pcCXSs3d+etHaYyatY7r+ibz50GnKewlboU1NMHMTjKzaWaWHvre8Ajtisxscejr+8sfipQLd+dP/1rJqFnruPGMNgp7iXvhjkW7D5jh7h2BGaHtshx0956hr8vC7FPkqIqLnQcmLue1rzdwy1nteHjQqQp7iXvhBv4g4PXQ49eBy8N8PpGwFRc7v/tgGeNmb+L2c0/mgYFdMFPYi4Qb+M3cfWvo8Tag2RHa1TCz+WY228wu/7EnNLMRobbzd+zYEWZ5Em+Kip173lvCO/M388ufduC+/p0V9iIhR71pa2bTgeZlHLq/9Ia7u5n5EZ6mjbtnmdnJwKdmtszd15bV0N1HAaMAUlJSjvR8Ij9wuKiYX7+zmMlLt/LrCzpx9wUdgy5JJKocNfDd/YIjHTOz7WbWwt23mlkLIPsIz5EV+r7OzD4DegFlBr7IicgvLOKXby3ik5Xbubd/Z+78SfugSxKJOuFe0pkEDAs9HgZM/H4DM2toZtVDjxsDZwErw+xX5DuHDhdx+xsL+GTldv54aVeFvcgRhBv4I4ELzSwduCC0jZmlmNnoUJsuwHwzWwLMBEa6uwJfIiIvv5CbX53H52t28NcruzH8rHZBlyQStcL64JW77wLOL2P/fOC20OOvgW7h9CNSln2HDnPLq/NYuGkPT1zdgyt6JQVdkkhU0ydtJSbtzitg2Ji5rNq6j+eu682Abi2CLkkk6inwJeZs33eIG0bPYdPuA/zjphTO69w06JJEYoICX2LK5t0HuH70HHbl5vPazan0a98o6JJEYoYCX2JGRnYuN4yew4GCQsbd1pdeyWVO3SQiR6DAl5iwPGsvw8bMxQzeub0fXVrUC7okkZijhTwl6s1dv5uho2ZTvWoVhb1IGHSGL1Ht07Tt3DluIa0a1mTcrX1p2aBm0CWJxCwFvkStiYuz+M27S+jcoi6v35xKozrVgy5JJKYp8CUqvfHNBh6ctILUticxelgKdWtUC7okkZinwJeo4u48NT2dp2ekc0GXpjx3XW9qVEsIuiyRSkGBL1GjqNh5aNJyxs3exFV9kvjrld2olqBxBSKRosCXqJBfWMSv31nMlGXbuOPc9tzb/xQtXCISYQp8Cdz+Q4cZMXYB36zbxR8u6cJt55wcdEkilZICXwKVve8QN782j9Xb9vPUNT25vFeroEsSqbQU+BKYjOz9DBszjz0HCnhl+Omc26lJ0CWJVGph3REzsyFmtsLMis0s5Ufa9Tez1WaWYWb3hdOnVA7zNuxm8IvfkF9YzLu391PYi1SAcIdALAeuBGYdqYGZJQDPAxcDXYGhZtY1zH4lhn20bCvXj55Do9qJTPjFmZzWqn7QJYnEhXBXvFoFHG00RSqQ4e7rQm3HA4PQurZx6dWv1vPw5JX0at2A0cNO56TaiUGXJBI3KuIafitgc6ntTKDvkRqb2QhgBEBycnL5ViYVpqjYeWTySl77egMXdW3GM0N76QNVIhXsqIFvZtOB5mUcut/dJ0a6IHcfBYwCSElJ8Ug/v1S8vPxC7h6/iOmrsrn17Hb8fkAXEqpojL1IRTtq4Lv7BWH2kQW0LrWdFNoncWD7vkPc+vo8Vm7Zx8ODTuWmfm2DLkkkblXEJZ15QEcza0dJ0F8LXFcB/UrA0rbt45ZX55Fz8DCjh6Xw087Ngi5JJK6FOyzzCjPLBPoB/zazj0P7W5rZFAB3LwTuAj4GVgHvuvuK8MqWaDdj1XYGv/A1Re78845+CnuRKBDuKJ0JwIQy9m8BBpTangJMCacviQ3uzqhZ6xg5NY3TWtbnHzel0Lx+jaDLEhH0SVuJoPzCIu6fsJz3FmRySbcW/G1ID2omaiSOSLRQ4EtE7MzN585xC5i3YQ93n9+Ru8/vSBWNxBGJKgp8CdvyrL3c/sYCdubm8+zQXlzao2XQJYlIGRT4EpaJi7O49/2lNKyVyD/v6Ef3pAZBlyQiR6DAlxNSWFTMo1PT+McX60ltdxIvXN+bxlpkXCSqKfDluOUcKOCXby/ii/Sd3NSvDQ8M7KqlCEVigAJfjsuyzL3cMW4BO/bn89jg7lx9euuj/5CIRAUFvhwTd2f8vM08NHEFTepW5907+tGzdYOgyxKR46DAl6M6WFDEAxNLxtef07ExT1/bS9Mai8QgBb78qA0787jzzYWkbdvH3ed35Ffnd9RMlyIxSoEvRzRxcRb3T1hO1QRjzPDTOe+UpkGXJCJhUODLDxwsKOJP/1rB+HmbSWnTkKeH9qJVg5pBlyUiYVLgy/9nzfb93PXWQtKzc/nP89rz6ws6UVVDLkUqBQW+ACWjcN6au4lHJq+kTvWqjL0llXM6Ngm6LBGJIAW+sDM3n3vfW8qMtGzO6diYv1/dg6Z1NaWxSGWjwI9zn6Zt57fvLWXfoUIeHNiV4We21SyXIpVUWIFvZkOAPwJdgFR3n3+EdhuA/UARUOjuKeH0K+E7UFDI/05ZxbjZm+jcvC5v3nYGpzSvG3RZIlKOwj3DXw5cCbx8DG3Pc/edYfYnETBn3S5++/5SNu46wM/Pacc9PzuF6lW1UIlIZRfuEoerAMx0CSAWHCgo5LGpq3n9mw0kNazJ2z8/g37tGwVdlohUkIq6hu/AJ2bmwMvuPupIDc1sBDACIDk5uYLKq/zmrt/N/7y3hI27DjCsXxt+278ztavrFo5IPDnqb7yZTQeal3HofnefeIz9nO3uWWbWFJhmZmnuPqushqE/BqMAUlJS/BifX45g36HDPD51NePmbNRZvUicO2rgu/sF4Xbi7lmh79lmNgFIBcoMfIkMd2fq8m08NGkFO3PzGX5mW+656BSd1YvEsXL/7Tez2kAVd98fenwR8HB59xvPtuQc5MGJy5m+KpuuLeoxeliKlh4UkbCHZV4BPAs0Af5tZovd/Wdm1hIY7e4DgGbAhNCN3arAW+4+Ncy6pQz5hUWM+XIDz36ajjvcP6ALN5/VVlMjiAgQ/iidCcCEMvZvAQaEHq8DeoTTjxzdzLRsHp68kvU787iwazMeHNiV1ifVCrosEYkiuqAb4zbszOORySuZkZbNyY1r8/otqZzbSXPgiMgPKfBj1J68Ap6bmcHYbzaQmFCF3w/ozPAz25FYVZdvRKRsCvwYc+hwEa9/vYHnZmaQl1/IkD6t+e+LOtGsniY7E5Efp8CPEUXFzoeLsnhi2hqycg7y085Nubd/Z81/IyLHTIEf5YqKnclLt/D0jHTW7cijW6v6PD6kO2e2bxx0aSISYxT4Uaq42Jm8bCvPzEgnIzuXzs3r8tINvbmoa3NNXywiJ0SBH2XyC4uYuHgLo2atIyM7l07N6vDC9b3pf6qCXkTCo8CPEvsOHebtOZsY89V6tu/Lp0uLejw7tBeXdGuhoBeRiFDgB2zjrjzGzd7I+Lmb2Z9fyFkdGvH4VT04p2NjTTstIhGlwA9AUbHz2epsxn6zkc/X7CChinHxac25/T/a0y2pftDliUglpcCvQJl7DjBhYRbj520mK+cgzepV578u6MjQ1GSNoxeRcqfAL2e5+YV8tGwr7y/MZPa63QD0O7kRf7ikCxd0bUY1TWwmIhVEgV8OcvML+Wx1Nh8t38aMVds5dLiYdo1r85sLO3F5r1aa1ExEAqHAj5BdufnMXL2Dqcu3MSt9BwWFxTSuk8jg3klc2TuJ3skNdBNWRAKlwD9B+YVFLNiwhy8ydvJF+g6WZ+0DoGX9GlzfN5mLT2tBnzYNSdCQShGJEgr8Y7Qnr4BFm/ewcGMOCzbuYfHmHA4eLqJqFaN3m4bcc1En/qNTE7q1qq8zeRGJSuGuePU4cClQAKwFbnb3nDLa9QeeBhIoWQlrZDj9lqfiYidzz0FWb9/P6m37SNu2n5Vb97FuRx4ACVWMri3qcc3prTm7Q2POaN+IOlonVkRiQLhJNQ34nbsXmtmjwO+Ae0s3MLME4HngQiATmGdmk9x9ZZh9H7eCwmJyDhSw+0ABe/IOs+dAAVl7DpK55wCZew6Suecgm/cc4EBB0Xc/0/qkmnRuXo/BvZPo06Yh3ZPqUytRAS8isSfcJQ4/KbU5G7iqjGapQEZoqUPMbDwwCCi3wL/kmS/Iyy/kcJGTX1jM4aJiCgqLOXi4qMz2datXJemkWiQ3qsWZHRrRqVldTmlel07N6ursXUQqjUim2S3AO2XsbwVsLrWdCfQ90pOY2QhgBEBycvIJFdKpWV2Kip3EqlWollCFxAQjsWoV6tWoRoPaiZxUK5GGtarRsHYiLRvUpH7NaifUj4hILDlq4JvZdKB5GYfud/eJoTb3A4XAm+EW5O6jgFEAKSkpfiLP8eQ1PcMtQ0Sk0jlq4Lv7BT923MyGAwOB8929rIDOAlqX2k4K7RMRkQoU1uf6Q6Nvfgtc5u4HjtBsHtDRzNqZWSJwLTApnH5FROT4hTuRy3NAXWCamS02s5cAzKylmU0BcPdC4C7gY2AV8K67rwizXxEROU7hjtLpcIT9W4ABpbanAFPC6UtERMKjqRpFROKEAl9EJE4o8EVE4oQCX0QkTljZQ+ejg5ntADaWw1M3BnaWw/NWJL2G6KDXEB30Gv5PG3dvUtaBqA788mJm8909Jeg6wqHXEB30GqKDXsOx0SUdEZE4ocAXEYkT8Rr4o4IuIAL0GqKDXkN00Gs4BnF5DV9EJB7F6xm+iEjcUeCLiMSJuA18M3vEzJaGZvn8xMxaBl3T8TKzx80sLfQ6JphZg6BrOl5mNsTMVphZsZnF1LA6M+tvZqvNLMPM7gu6nuNlZmPMLNvMlgddy4kys9ZmNtPMVobeR3cHXdPxMrMaZjbXzJaEXsOfyq2veL2Gb2b13H1f6PGvgK7ufkfAZR0XM7sI+LTUIvK4+71H+bGoYmZdgGLgZeAed58fcEnHxMwSgDXAhZQs2zkPGOru5bZWc6SZ2X8AucBYdz8t6HpOhJm1AFq4+0IzqwssAC6PsX8HA2q7e66ZVQO+BO5299mR7ituz/C/DfuQ2kDM/eVz909C6w1AySLySUHWcyLcfZW7rw66jhOQCmS4+zp3LwDGA4MCrum4uPssYHfQdYTD3be6+8LQ4/2UrLnRKtiqjo+XyA1tVgt9lUsexW3gA5jZX8xsM3A98GDQ9YTpFuCjoIuII62AzaW2M4mxoKlszKwt0AuYE3Apx83MEsxsMZANTHP3cnkNlTrwzWy6mS0v42sQgLvf7+6tKVl8/a5gqy3b0V5DqE3EFpEvD8fyGkTCYWZ1gPeB//re/95jgrsXuXtPSv6Xnmpm5XKJLawVr6Ld0RZgL+VNSlbkeqgcyzkhEVhEPnDH8e8QS7KA1qW2k0L7pIKFrnu/D7zp7h8EXU843D3HzGYC/YGI30yv1Gf4P8bMOpbaHASkBVXLiTrGReSlfMwDOppZOzNLBK4FJgVcU9wJ3fB8BVjl7k8EXc+JMLMm346wM7OalAwEKJc8iudROu8Dp1AyQmQjcIe7x9QZmpllANWBXaFds2NwpNEVwLNAEyAHWOzuPwu0qGNkZgOAp4AEYIy7/yXYio6Pmb0N/ISSaXm3Aw+5+yuBFnWczOxs4AtgGSW/ywC/D62jHRPMrDvwOiXvoyrAu+7+cLn0Fa+BLyISb+L2ko6ISLxR4IuIxAkFvohInFDgi4jECQW+iEicUOCLiMQJBb6ISJz4fzJKpj7SjGFEAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAuSElEQVR4nO3dd3hUZdrH8e+dQhI6IaEnJJDQFGkx2KULuyLYwXUXXZW1oO5aXlnrrq77Wnata8O14FqwL6wiSBVWQAhKhxRCCy0hkFASSLvfPzLxjZgQwkzmTLk/1zUXM6fM/EYnc89znnOeR1QVY4wxwSvE6QDGGGOcZYXAGGOCnBUCY4wJclYIjDEmyFkhMMaYIBfmdIBTERMTowkJCU7HMMYYv7Jy5cp9qhp7/HK/LAQJCQmkpaU5HcMYY/yKiGyrabkdGjLGmCBnhcAYY4KcFQJjjAlyVgiMMSbIWSEwxpgg55FCICJvikiuiKyrZb2IyAsikiUia0Skf7V1E0Qk03Wb4Ik8xhhjTp6nWgRvAyNPsH4UkOy6TQReARCRaOARYCCQCjwiIq08lMkYY8xJ8Mh1BKq6SEQSTrDJGOAdrRzzepmItBSR9sAgYI6q7gcQkTlUFpQPPJHLGOObKiqUvMPH2LG/iJwDxewsKCYsRGjVpBGtGjeiVeNwops0onPrJoSGiNNxA563LijrCOyo9jjHtay25T8jIhOpbE0QHx/fMCmNMQ0mc+8hvlq3h6837CFjz2FKyivq3KdFVDjnJcVwQbcYzk+OpUPLKC8kDT5+c2Wxqk4BpgCkpKTYbDrG+IHNeYf59w87+WrdHrJyDyMCA+Jbcf15CcS1akynVlF0atWYji2jKFflwJESDhSVcKColNyDR/luy34WZ+bx5drdACS3aco1A+MZd2Y8UY1CHX53gcNbhWAnEFftcSfXsp1UHh6qvnyhlzIZYxrI9vwinpubweerdiLAwMTWTDi7Mxed1o42zSNr3a9pRBhx0Y1/fHxlShyqSsbewyzOzOOrdXv483828OL8LK4/J4HfnJ1Ai8bhXnhHgU08NVWlq4/gC1U9vYZ1vwQmAb+gsmP4BVVNdXUWrwSqziL6HhhQ1WdQm5SUFLWxhozxPXsKj/Li/Ew+XLGD0BDhunMSuPH8LsQ2i/DYa6zYup+XF2SxID2PphFh/OqseG4bnETzSCsIdRGRlaqacvxyj7QIROQDKn/Zx4hIDpVnAoUDqOqrwEwqi0AWUARc71q3X0QeA1a4nurRuoqAMcb3lJZX8OL8LF77ZjMVqoxPjWfSkCTanuDX/6k6MyGat65PZf2uQl5ZuJnXF2UzY9Uunrj8DC7s9rOBNc1J8FiLwJusRWCM79i67wh3friK1TsKGNO3A/eM6P6TwzsN7YftB7jn49VszjvCuDPjuP+XPa11UIsGbREYY4KPqvLJyhz+NGM9oSHCS9f055dntPd6jn7xrfjyjvN5dm4Gry/K5puMPGsd1JMNMWGMqbfColImvf8D936yht6dWjDr9xc4UgSqRIaH8sdRPfn0lnNo3CiUCW8u58V5mfjjEQ8nWIvAGFMvO/YXMeHN5WzfX8T/jOzO7y7o6jMXfVW1Dv742Vr+PieDLflH+N/LehMRZqeanogVAmPMSduw6yAT3lrOsdJy3rtxIAO7tHY60s9EhofyzFV9SIxpwjNzMsg5UMxr1w6gVZNGTkfzWXZoyBhzUpZs3sfVry0lLET45JZzfLIIVBER7hiazPPj+rJqewGXvbKELfuOOB3LZ1khMMbU6Ys1u7juzRW0axHJp7ecQ7e2zZyOdFLG9O3I+zcNpLC4lEtf/pYNuw46HcknWSEwxpzQB8u3c/sHP9AnrgWf3HyO3433k5IQzee3nkNUeCjXvvEdGXsPOR3J51ghMMbUasbqXdz/+Vou7BbLv24Y6LfDOXRu3YT3bzqLsBDhmte/IzvvsNORfIoVAmNMjRZsyuWuD1dxZkI0r147gMhw/z7zJjGmCe/fNBBV5ZrXv2N7fpHTkXyGFQJjzM8s37KfW95bSY/2zfjnhBS/LwJVkto0490bB3K0rJzxry9jZ0Gx05F8ghUCY8xPrNtZyA1vr6BDyyimXp8acMM19GzfnHdvGMjBo6X86vVlHDhS4nQkx1khMMb8aOu+I1z31nKaR4Xz7g0Dad3Uc6OG+pLTO7bg7etT2VV4lJvfXUlJWd2T5AQyKwTGGAAOHS3lxnfSqFD41w2pfnd2UH0N6NyKp684g++27OeRGeuCejgKu7LYGENFhfKHD1exZd8R3r1hIF1imzodySvG9O1I5t7D/GNBFkltmnHDeYlOR3KEtQiMMTw3N4O5G3N5+OJenN3Vd68Ybgh3De/GRae15fEvN7AgPdfpOI6wQmBMkJu1bjcvzM/iqpRO/Obszk7H8bqQEOHZq/vSo11zbn//h6C84MwjhUBERopIuohkicjkGtY/KyKrXLcMESmotq682roZnshjjDk5m/Yc5K6PVtM3riWPjT0dEd8YRdTbGjcK+/E02ZveSePg0VKnI3mV24VAREKBl4BRQC9gvIj0qr6Nqv5BVfuqal/gReCzaquLq9ap6iXu5jHGnJyCohImvrOSphFhvPbrAUE/VHOHllG8em1/cg4Uc/9na4Oq89gTLYJUIEtVs1W1BJgGjDnB9uOBDzzwusaYU6Sq3PfpGnYXFvPKtQMaZG5hf5SSEM1dw7vxxZrdTFuxw+k4XuOJQtARqP5fLMe17GdEpDOQCMyvtjhSRNJEZJmIjK3tRURkomu7tLy8PA/ENiZ4fbB8B7PX7+Xei7ozoHMrp+P4lFsu7Mr5yTH8acZ60vcER3+BtzuLxwGfqGp5tWWdXZMpXwM8JyJda9pRVaeoaoqqpsTG2lykxpyqrNxDPPrFes5PjuHG87o4HcfnhIQIz1zVl2aR4Ux6/3uKS8rr3snPeaIQ7ATiqj3u5FpWk3Ecd1hIVXe6/s0GFgL9PJDJGFODY2Xl3P7BKho3CuPvV/YhxEemmPQ1sc0ieO7qvmTlHeZPM9Y7HafBeaIQrACSRSRRRBpR+WX/s7N/RKQH0ApYWm1ZKxGJcN2PAc4FNnggkzGmBk/NSmfj7oM8fcUZtLF+gRM6LzmGWwd15cO0HUxfVdtv28DgdiFQ1TJgEjAb2Ah8pKrrReRREal+FtA4YJr+tCu+J5AmIquBBcATqmqFwJgGsDA9lzf+u4UJZ3dmaM+2TsfxC38Y1o2Uzq148PN17C4M3JFKxR9PkUpJSdG0tDSnYxjjN/YdPsbI5xbRukkE0yedGzDDSnvD9vwiLnpuEamJ0bx9/Zl+fa2FiKx09cn+hF1ZbEwQeGTGeg4Wl/H8+L5WBOopvnVjJo/qwTcZeXycluN0nAZhhcCYADdr3R6+XLObO4Ym0aNdc6fj+KVfn9WZs7pE89gXG9gVgJPZWCEwJoAVFJXw0PR19GrfnN9dWOOZ2eYkhIQIT13eh3JVJgfgVcdWCIwJYI99sZH9R0p46oozCA+1P3d3VB0iWpSRx4cBdtWxfTKMCVAL03P59PscbrmwK6d3bOF0nIBw7cDKQ0R/+XJjQM13bIXAmAB06Ggp93+2lqQ2Tbl9aJLTcQJGSIjw9BV9qFANqIHprBAYE4CenLWJ3QeP8tQVZwT9qKKeFhfdmHsv6s43GXl8tW6P03E8wgqBMQEmbet+3l22nd+em0j/eBtQriH8+qzOnNahOX/+z3oOBcDcBVYIjAkgZeUVPPjvdXRoEcndI7o5HSdghYWG8Pilvck9dIxn5mQ4HcdtVgiMCSBTl25j055DPDy6F40bhTkdJ6D1jWvJtQM7M3XJVtbtLHQ6jlusEBgTIPYePMqzczK4sFssF53Wzuk4QeGei7oT3SSCBz5fS3mF/3YcWyEwJkA8/uVGSsor+PMlp/n1eDj+pEVUOA9d3JPVOYW8v3y703FOmRUCYwLAkqx9zFi9i5sv7EpCTBOn4wSVS/p04LykGJ6atYncQ0edjnNKrBAY4+dKyip4aPo64qKjuHWQDSPhbSLCo2NO41hpBU/M3OR0nFNihcAYP/fGf7ewOe8If77kNBtZ1CFdYpty0wWJfPbDTr7ffsDpOPVmhcAYP7a7sJgX5mUyrGdbhvSwyWacdOugJNo0i+DP/9lAhZ91HHukEIjISBFJF5EsEZlcw/rrRCRPRFa5bjdWWzdBRDJdtwmeyGNMsHh6VjrlFcrDF/dyOkrQaxIRxuRRPVi9o4DPf/CvqS3dLgQiEgq8BIwCegHjRaSmT+WHqtrXdfuna99o4BFgIJAKPCIidimkMSdh1Y4CPvthJzecn0h868ZOxzHA2L4d6RPXkidnbeLIsTKn45w0T7QIUoEsVc1W1RJgGjDmJPe9CJijqvtV9QAwBxjpgUzGBDRV5bEvNhDTtJF1EPuQkBDhkdG9yD10jJcXZjkd56R5ohB0BKoPzp3jWna8y0VkjYh8IiJx9dwXEZkoImkikpaXl+eB2Mb4ry/W7GbltgPcM6I7zSLDnY5jqukf34pL+3Xk9cVb2J5f5HSck+KtzuL/AAmqegaVv/qn1vcJVHWKqqaoakpsbKzHAxrjL46WlvPEV5vo2b45V6bE1b2D8br7RvYgVIS/ztzodJST4olCsBOo/mns5Fr2I1XNV9Vjrof/BAac7L7GmJ96479b2FlQzEMX9yQ0xK4g9kXtWkRy2+CuzFq/hyVZ+5yOUydPFIIVQLKIJIpII2AcMKP6BiLSvtrDS4CqMjkbGCEirVydxCNcy4wxNcg9eJSXFmQxoldbzuka43QccwI3nt+Fji2jeHzmRp8/ndTtQqCqZcAkKr/ANwIfqep6EXlURC5xbXaHiKwXkdXAHcB1rn33A49RWUxWAI+6lhljavC3r9MpLa/g/l/0dDqKqUNkeCj3XtSd9bsOMmP1LqfjnJD441RrKSkpmpaW5nQMY7xq4+6D/OKFxdxwbiIP2nUDfqGiQhn9j/9SUFTKvLsvdPzKbxFZqaopxy+3K4uN8RNPztpEs4gwJg2xOYj9RUiIcP8verKzoJh/Ld3mdJxaWSEwxg8s2byPhel53DY4iZaNGzkdx9TDuUkxXNgtlhfnZ1JQVOJ0nBpZITDGx6kqT361ifYtIplwToLTccwpmDyqB4eOlfHyws1OR6mRFQJjfNzMtXtYnVPIXcO7OX6M2Zyanu2bc0X/Trz97VZ27Pe9i8ysEBjjw0rLK3h69ia6t23GZf07OR3HuOGuEd0QwScnu7dCYIwPm7Z8O1vzi7hvVHe7eMzPtW8RxQ3nJfL5Dzt9brJ7KwTG+KjDx8p4fl4mAxOjGdy9jdNxjAfcPKgrrRqH89TsdKej/IQVAmN81D8XZ7PvcAmTR/WwyegDRPPIcG4Z1JVFGXksy853Os6PrBAY44PyDx/j9UXZjDq9Hf3ibYqOQPKbsxNo2zyCp2en4ysX9FohMMYHvbJwM8Wl5dw9orvTUYyHRYaHcsfQZFZuO8CC9Fyn4wBWCIzxObsLi3ln2TYu69+JpDZNnY5jGsBVKXF0bt2Yp2dn+MSAdFYIjPExL87PQlW5c2iy01FMAwkPDeGu4d3YuPsgX6zd7XQcKwTG+JJt+Uf4aMUOxqfGExdt8xAHstFndKBHu2Y84xpR1klWCIzxIc/NzSQsVJg02AaWC3QhIcI9I7qzNb+Ij9NynM3i6KsbY36UsfcQ/161kwnnJNCmeaTTcYwXDO3Zhv7xLXlhXiZHS8sdy+GRQiAiI0UkXUSyRGRyDevvEpENrsnr54lI52rrykVkles24/h9jQkWf/86naaNwrj5gq5ORzFeIiLce1EP9hw8yrvLnBum2u1CICKhwEvAKKAXMF5Ejp814wcgxTV5/SfAU9XWFatqX9ftEowJQqt3FDB7/V5uPL8LrZrYMNPB5OyurTk3qTWvfrOZopIyRzJ4okWQCmSparaqlgDTgDHVN1DVBapaNeTeMionqTfGuPx9TgatGofz2/MSnI5iHHDX8O7sO1zC1CXOtAo8UQg6AjuqPc5xLavNDcBX1R5HikiaiCwTkbG17SQiE13bpeXl5bkV2BhfkrZ1P4sy8rj5wq40iwx3Oo5xwIDOrRjUPZbXFm3m0NFSr7++VzuLReRaIAV4utrizq45NK8BnhORGg+QquoUVU1R1ZTY2FgvpDXGO56dm0FM0wh+c3aC01GMg/4wrBsFRaW8/e1Wr7+2JwrBTiCu2uNOrmU/ISLDgAeAS1T1WNVyVd3p+jcbWAj080AmY/zCsux8vs3K55ZBXYlqZJPOBLM+cS0Z1rMtry/OprDYu60CTxSCFUCyiCSKSCNgHPCTs39EpB/wGpVFILfa8lYiEuG6HwOcC2zwQCZjfJ6q8sycDNo0i+BXA+OdjmN8wF3Du3HwaBlvLM726uu6XQhUtQyYBMwGNgIfqep6EXlURKrOAnoaaAp8fNxpoj2BNBFZDSwAnlBVKwQmKCzdnM/yLfu5dVBXm4LSANCrQ3N+0bsdb367lQNHvDfRfZgnnkRVZwIzj1v2cLX7w2rZbwnQ2xMZjPEnVa2Bds0jGZdqrQHz/34/rBtfrdvDlMXZ3Deyh1de064sNsYBizP3kbbtALcNSbLWgPmJbm2bMfqMDrz97Vb2HT5W9w4eYIXAGC+rag10bBnFVSl2SY35uTuHJXOsrJwpi7zTV2CFwBgvW5iex6odBUwakkREmLUGzM91jW3KmL4deWepd1oFVgiM8SJV5dm5GXRqFcUVA6w1YGp3+5AkSsoqvNIqsEJgjBctSM9lTU4hdwxJJjzU/vxM7brENmWsq1WQd6hhWwX2STTGS1SV5+ZmEhcdxaX9TzQKizGVJv3YKtjcoK9jhcAYL6lqDdw+2FoD5uR0iW3K2H4d+deybQ3aKrBPozFeYK0Bc6puH5JMSVkFr33TcK0CKwTGeIG1BsypSoxpwth+HXn3u23kHjraIK9hn0hjGpi1Boy77hiSTGm5MuWbhjmDyAqBMQ3MWgPGXQkxTRjbt+FaBR4Za8gYUzNrDRhPuX1IEkUlZZSUVXj8ua0QGNOAqloDT11+hrUGjFsSYprwyrUDGuS57ZNpTAOx1oDxF1YIjGkgVa2BSYOTrDVgfJp9Oo1pAKrK83Mz6dQqisv625hCxrd5pBCIyEgRSReRLBGZXMP6CBH50LX+OxFJqLbuj67l6SJykSfyGOO0hRl5rM4p5DZrDRg/4PYnVERCgZeAUUAvYLyI9DpusxuAA6qaBDwLPOnatxeVcxyfBowEXnY9nzF+q6o10LFlFJdba8D4AU/8VEkFslQ1W1VLgGnAmOO2GQNMdd3/BBgqIuJaPk1Vj6nqFiDL9XwN4vVF2fxtdnpDPb0xACzK3MeqHQXcNjiJRmHWGjC+zxOf0o7AjmqPc1zLatzGNdl9IdD6JPcFQEQmikiaiKTl5eWdUtDsfUeYsiib3YXFp7S/MXWpbA1Uzj5m8w0Yf+E3P1dUdYqqpqhqSmxs7Ck9x22Du1KhyisLG3ZIVxO8/pu1j++3F3DLoK7WGjB+wxOf1J1AXLXHnVzLatxGRMKAFkD+Se7rMZ1aNebKlE5MW76DPYUNM3iTCV5VfQPtW0Rypc1FbPyIJwrBCiBZRBJFpBGVnb8zjttmBjDBdf8KYL6qqmv5ONdZRYlAMrDcA5lqdeugJFerIKshX8YEoSWb80nbdoBbB3W1uYiNX3G7ELiO+U8CZgMbgY9Udb2IPCoil7g2ewNoLSJZwF3AZNe+64GPgA3ALOA2VS13N9OJxEU35ooBnfhgxQ72HrRWgfGMqtZAu+aRXHVmXN07GONDPHIQU1Vnqmo3Ve2qqo+7lj2sqjNc94+q6pWqmqSqqaqaXW3fx137dVfVrzyRpy63DU6iosL6CoznLN2cz/Kt+7ltsLUGjP8Jyt6suOjGXNa/I+8v326tAuO2qjGFrDVg/FVQFgKASYOTKa9QXm3A6d9McKhqDdxqrQHjp4K2EMS3bsxl/Try/nfbybVWgTlFP2kNpFhrwPinoC0EAJOGJFFWobzaQNO/mcC3NLuyNXDLoK5EhltrwPinoC4EnVs34dJ+HXnvu23WKjD1VtUaaNs8gqutb8D4saAuBFA5/VtZhfKK9RWYelqanc/yLfu5dVCStQaMXwv6QtC5dRMu69eR976zM4jMybPWgAkkQV8IAG4fUnkGkV1XYE5WVWvglgutb8D4PysEVJ5BdLnrugIbg8jURVV5bk5la2BcarzTcYxxmxUCl0mDk6mw6wrMSVhSdd2A9Q2YAGGFwKWyVdDJWgXmhFSVZ+Zk0L5FpPUNmIBhhaCaSUMqxyB62UYmNbVYlLmPldsOcNtgaw2YwGGFoJqqkUmnLd9hs5iZn1FVnp1TOfuYXUVsAokVguPcNrhyvoKXFlirwPzUwvQ8Vu0oYNIQm4vYBBb7NB8nLroxV50Zx4crdpBzoMjpOMZHqCrPzs2gUyubi9gEHisENZg0OAlBeHGetQpMpXkbc1mTU8gdQ5IJD7U/GxNY3PpEi0i0iMwRkUzXv61q2KaviCwVkfUiskZErq627m0R2SIiq1y3vu7k8ZQOLaO4ZmA8n3yfw9Z9R5yOYxxW1RqIj27Mpf07Oh3HGI9z96fNZGCeqiYD81yPj1cE/EZVTwNGAs+JSMtq6+9V1b6u2yo383jMrYO6EhYivDAv0+koxmFfb9jL+l0HuWOotQZMYHL3Uz0GmOq6PxUYe/wGqpqhqpmu+7uAXCDWzddtcG2aRzLhnAQ+X7WTrNxDTscxDimvUJ75OoMuMU0Y27eD03GMaRDuFoK2qrrbdX8P0PZEG4tIKtAIqH757uOuQ0bPikjECfadKCJpIpKWl5fnZuyT87sLuhAVHsqzc61VEKy+WLOL9L2H+P3wboRZa8AEqDo/2SIyV0TW1XAbU307VVVAT/A87YF/AderaoVr8R+BHsCZQDRwX237q+oUVU1R1ZTYWO80KFo3jeD6cxP4cs1uNu4+6JXXNL6jrLyCZ+dk0KNdMy7u3d7pOMY0mDoLgaoOU9XTa7hNB/a6vuCrvuhza3oOEWkOfAk8oKrLqj33bq10DHgLSPXEm/Kkied3pVlkGM/MyXA6ivGyT7/PYWt+EXeP6E5IiDgdx5gG425bdwYwwXV/AjD9+A1EpBHwOfCOqn5y3LqqIiJU9i+sczOPx7VoHM6N53Vhzoa9rMkpcDqO8ZJjZeW8MC+LPnEtGdazjdNxjGlQ7haCJ4DhIpIJDHM9RkRSROSfrm2uAi4ArqvhNNH3RGQtsBaIAf7iZp4G8dvzEmjZOJy/fW2tgmDxwXfb2VlQzD0julH5O8WYwBXmzs6qmg8MrWF5GnCj6/67wLu17D/Endf3lmaR4dw6qCt/nbmJpZvzObtra6cjmQZUXFLOPxZsZmBiNOclxTgdx5gGZ6dBnKTfnJ1Au+aRPDV7E5X94iZQTV26lX2Hj3HvRd2tNWCCghWCkxQZHsqdw5L5YXsBczfW2CduAsCho6W8+s1mBnWPJSUh2uk4xniFFYJ6uHJAJxJjmvD07E2UV1irIBC9viibgqJS7h7e3ekoxniNFYJ6CAsN4a7h3cjYe5jpq3Y6Hcd4WO6ho7y+eAsXn9Ge3p1aOB3HGK+xQlBPv+zdntM6NOeZORmUlFXUvYPxGy/My6S0vIJ7RlhrwAQXKwT1FBIi3HtRd3IOFPPB8u1OxzEesmXfEaYt38H41HgSYpo4HccYr7JCcAou7BZLamI0L87PoqikzOk4xgP+9nU64aEh3D40yekoxnidFYJTICLcN7I7+w4f443FW5yOY9y0JqeAL9fs5qbzE2nTLNLpOMZ4nRWCUzSgczQjerXl1W82k3fomNNxzClSVZ74ahPRTRpx0wVdnI5jjCOsELhh8qgeHCur4Pl5NvSEv1qcuY8lm/OZNDiJZpHhTscxxhFWCNzQJbYp1wyM54PlO8jKPex0HFNPFRXKk7M20alVFL86K97pOMY4xgqBm+4cmkxUeChPfLXJ6Simnqav3sn6XQe5e0Q3IsJCnY5jjGOsELipddMIbhnUlbkb9/Jddr7TccxJKi4p56lZ6fTu2IIxfWxCehPcrBB4wA3nJdK+RSR/nbmRCht6wi+8vjib3YVHeejiXjbpjAl6Vgg8IDI8lLtHdGd1TiFfrN1d9w7GUXsPHuWVhZsZdXo7UhNtYDlj3CoEIhItInNEJNP1b6tatiuvNinNjGrLE0XkOxHJEpEPXbOZ+aVL+3WkV/vmPDVrE8fKyp2OY07g71+nU1ZRweRRPZyOYoxPcLdFMBmYp6rJwDzX45oUq2pf1+2SasufBJ5V1STgAHCDm3kcExoi3P+LnuQcKOatb7c6HcfUYv2uQj5emcN15yTQubUNJWEMuF8IxgBTXfenUjnv8ElxzVM8BKiax7he+/ui85JjGNazDS/OyyT34FGn45jjqCp/+WIjLaPCmTQk2ek4xvgMdwtBW1WtOii+B2hby3aRIpImIstEZKxrWWugQFWrBuvJAfz+9I0Hf9mL0nLlyVnpTkcxx5m7MZel2fn8YXg3WkTZxWPGVKlzzmIRmQu0q2HVA9UfqKqKSG2nzHRW1Z0i0gWY75qwvrA+QUVkIjARID7edy/+SYhpwm/PS+TVbzZz7Vnx9IuvsdvEeFlJWQX/O3MjXWObMD7Vdz8/xjihzhaBqg5T1dNruE0H9opIewDXvzXO4aiqO13/ZgMLgX5APtBSRKqKUSeg1tleVHWKqqaoakpsbGw93qL3TRqSRJtmEfxpxno7ndRHvPntFrL3HeHBX/YiPNROljOmOnf/ImYAE1z3JwDTj99ARFqJSITrfgxwLrBBK2eAXwBccaL9/VHTiDAmj+rB6pxCPv0+x+k4QW93YTEvzMtkWM+2DO7Rxuk4xvgcdwvBE8BwEckEhrkeIyIpIvJP1zY9gTQRWU3lF/8TqrrBte4+4C4RyaKyz+ANN/P4jLF9O9IvviVPzkrn0NFSp+MEtce/3Eh5hfLI6F5ORzHGJ9XZR3AiqpoPDK1heRpwo+v+EqB3LftnA6nuZPBVISHCn0afxpiXvuUf87P44y96Oh0pKH2btY8v1uzm98OSiYtu7HQcY3ySHSxtQH3iWnLlgE68+e0WNufZ6KTeVlJWwSMz1hMf3ZibL+zqdBxjfJYVggb2PyN7EBUeyoOfr6OyW8R4y9tLtpCVe5hHRvciMtxGFzWmNlYIGlhsswjuG9WDpdn5fPZ9rSdFGQ/bU3iU5+dmMrRHG4b2rO3yFmMMWCHwivFnxtM/viV/+XID+4+UOB0nKDw+cyOlFcojo09zOooxPs8KgReEhAh/vaw3h46W8b8zNzodJ+AtTM/lP6t3ccuFXYlvbR3ExtTFCoGX9GjXnBvP78LHK3NYZhPYNJjDx8p44PN1JLVpyq2DrYPYmJNhhcCL7hyaTKdWUTzw+VobqrqB/G12OrsKi3ny8t42/aQxJ8kKgRdFNQrlsbGnsznvCK99k+10nICzctsBpi7dym/O6syAzjbhjDEnywqBlw3u3oZfntGefyzIIivXri3wlGNl5dz36RraN4/k3pE24Ywx9WGFwAGPjO5F40ah3P3xasrKK5yOExBeWrCZrNzDPH5Zb5pGuHXBvDFBxwqBA9o0i+SxMaezekcBry2yQ0Tu2rTnIK8szGJs3w4M7m6DyhlTX1YIHDK6Twd+eUZ7npubwYZdB52O47dKyyu475M1NIsM52G7ZsCYU2KFwEGPjTmdFlGNuOujVZSU2SGiU/Hi/CxW5xTy6JjTiG7SyOk4xvglKwQOim7SiP+9rDeb9hzixfmZTsfxOyu37ecf8zO5rF9HLj6jg9NxjPFbVggcNrxXW64Y0ImXF25m1Y4Cp+P4jUNHS/n9h6vo0DKKP4+xQ0LGuMMKgQ94eHQv2jaL4O6PVlFUUuZ0HL/wpxkb2HmgmOeu7kuzSJuI3hh3uFUIRCRaROaISKbr35/N1C4ig0VkVbXbUREZ61r3tohsqbaurzt5/FXzyHCevrIP2fuO8Mj09U7H8XlfrtnNp9/ncNvgJFIS7MIxY9zlbotgMjBPVZOBea7HP6GqC1S1r6r2BYYARcDX1Ta5t2q9qq5yM4/fOjcphtsHJ/Hxyhw+WWnzHNdmd2Ex93++lj5xLbljaLLTcYwJCO4WgjHAVNf9qcDYOra/AvhKVYvcfN2AdOewbpzVJZoH/72WjL2HnI7jc8rKK/jDh6soLa/guav7Eh5qRzaN8QR3/5Laqupu1/09QF0zgIwDPjhu2eMiskZEnhWRiNp2FJGJIpImIml5eXluRPZdoSHCC+P60TQinNve+976C47z1Ox0lmXv59Exp5MY08TpOMYEjDoLgYjMFZF1NdzGVN9OK+dhrHUuRhFpT+Uk9rOrLf4j0AM4E4gG7qttf1WdoqopqpoSGxtbV2y/1aZ5JM+P60tW3mEe+rf1F1T5z+pdTFmUza/P6swVAzo5HceYgFLnoCyqOqy2dSKyV0Taq+pu1xd97gme6irgc1UtrfbcVa2JYyLyFnDPSeYOaOcmxXDn0GSem5vJwC7RXJUS53QkR6XvOcT/fLKGAZ1b8dDFvZyOY0zAcffQ0Axgguv+BGD6CbYdz3GHhVzFAxERKvsX1rmZJ2DcPiSZc5Na89C/1/HD9gNOx3FMYXEpv/tXGk0jw3j5V/1pFGb9AsZ4mrt/VU8Aw0UkExjmeoyIpIjIP6s2EpEEIA745rj93xORtcBaIAb4i5t5AkZVf0Hb5pHc9E4aO/YHX/96RYXy+2k/kHOgmFd+1Z+2zSOdjmRMQJLKQ/v+JSUlRdPS0pyO4RVZuYe57OVvadcikk9uOYfmQXTx1N+/TufF+Vk8OuY0fnN2gtNxjPF7IrJSVVOOX27tbB+X1KYpr147gOy8I9z23veUBsn8Be8u28aL87O4KqUTvz6rs9NxjAloVgj8wDlJMfz10t4sztzHIzPW44+tuPr4au1uHpq+jiE92vD4pb2p7EIyxjQUm8rJT1x1Zhxb8o/wysLNJLRuzMQLujodqUEs3ZzPndNW0S+uJS9d098uGjPGC6wQ+JF7R3Rne34Rf525iahGYQF3yGT9rkImvpNGfOvGvHndmUQ1CnU6kjFBwQqBHwkJEZ65ug/Hysp56N+VZ9oGSjHYnl/EdW+toGlkGO/8NpWWjW2SGWO8xdrdfiYiLJSXftWfYT3b8NC/1/Husm1OR3JbVu5hrp6ylJKyCt75bSodWkY5HcmYoGKFwA9VFYOhPdrwoJ8XgzU5BVz12lJKy5UPbjqL5LbNnI5kTNCxQuCnIsJCefna/y8G7yzd6nSkeluyeR/jpyyjcaNQPrn5bHp1aO50JGOCkhUCP1ZVDIb1bMvD09fzpxnrKfOT6wy+Xr+H695aQYeWUXxy8zkk2GiixjjGCoGfiwgL5bVfD+DG8xJ5e8lWrntrBQVFJU7HqpWq8q+lW7nlve/p2b45H/3ubNq1sKEjjHGSFYIAEBoiPHhxL56+4gyWb9nP2Je+JSvX9ya2OXi0lEkf/MBD09dzfnIM7984kFZN7OwgY5xmhSCAXJkSxwcTB3L4WDljX1rCrHV7nI70o7U5hYx+8b/MWreH+0b24M0JZ9Ikws5eNsYXWCEIMAM6RzNj0rkkxjTh5ndXcut7K8k9eNSxPKrK299u4fJXllBSVsGHE8/ilkFdCQmxYSOM8RX2kywAdWgZxWe3nsOURdk8Py+TxZn7+OOonow7M86rX8CrdhTw1y83snzrfob2aMPfruxjh4KM8UE2DHWA27LvCPd/tpal2fmkJkTz8OhenN6xRYO+5o79RTw9O50Zq3cR07QRdw3vzvjUOBs8zhiH1TYMtRWCIKCqfLwyh8e/3EhhcSlnd2nNb89LZGiPNh5tIewuLObtb7fy1rdbCQmBm87vwu8u7EpT6wswxic0SCEQkSuBPwE9gVRVrfHbWURGAs8DocA/VbVqJrNEYBrQGlgJ/FpV6zz30QrBqSksKmXaiu1MXbKVXYVHSWjdmOvPTWR0nw5En+Ihm4NHS5m1dg+f/7CTZVvyAbi8fyfuHtGN9i1sqAhjfElDFYKeQAXwGnBPTYVAREKBDGA4kAOsAMar6gYR+Qj4TFWnicirwGpVfaWu17VC4J7S8gpmrdvDG//dwqodBQAkt2lKamI0qYnRnJkQTdvmkYQe11ooK69g2/4iNuceJivvMGt2FDI/PZeSsgoSY5owtm9HxvbrQOfWdnGYMb6otkLgVptdVTe6nvxEm6UCWaqa7dp2GjBGRDYCQ4BrXNtNpbJ1UWchMO4JDw1hdJ8OjO7TgbU5hSzOymP5lv1MX7WL977b/uN2EWEhRDUKpXF4KOFhIewqKKa0/P9/OHRoEck1qfGM7deRPp1aWB+AMX7KGwdvOwI7qj3OAQZSeTioQFXLqi3vWNuTiMhEYCJAfHx8wyQNQr07taB3pxbcOqjyF/+mPYdYue0AB4pKKC4tp7iknKKSckrKKhh1enuS2jQlqU1TusY2oVkQzZ9sTCCrsxCIyFygXQ2rHlDV6Z6PVDNVnQJMgcpDQ9563WASFhrC6R1bNPhZRcYY31JnIVDVYW6+xk4grtrjTq5l+UBLEQlztQqqlhtjjPEib1xZvAJIFpFEEWkEjANmaGUv9QLgCtd2EwCvtTCMMcZUcqsQiMilIpIDnA18KSKzXcs7iMhMANev/UnAbGAj8JGqrnc9xX3AXSKSRWWfwRvu5DHGGFN/dkGZMcYEidpOH7VB54wxJshZITDGmCBnhcAYY4KcFQJjjAlyftlZLCJ5wLYGeOoYYF8DPK832XvwDfYefEMgvAfw3PvorKqxxy/0y0LQUEQkraYedX9i78E32HvwDYHwHqDh34cdGjLGmCBnhcAYY4KcFYKfmuJ0AA+w9+Ab7D34hkB4D9DA78P6CIwxJshZi8AYY4KcFQJjjAlyVgiOIyKPicgaEVklIl+LSAenM9WXiDwtIptc7+NzEWnpdKb6EpErRWS9iFSIiF+d/iciI0UkXUSyRGSy03nqS0TeFJFcEVnndJZTJSJxIrJARDa4Pkd3Op2pvkQkUkSWi8hq13v4c4O9lvUR/JSINFfVg677dwC9VPVmh2PVi4iMAOarapmIPAmgqvc5HKteRKQnUAG8Btyjqn4x3KyIhAIZwHAqp19dAYxX1Q2OBqsHEbkAOAy8o6qnO53nVIhIe6C9qn4vIs2AlcBYP/v/IEATVT0sIuHAf4E7VXWZp1/LWgTHqSoCLk0Av6uUqvp1tbmgl1E5+5tfUdWNqprudI5TkApkqWq2qpYA04AxDmeqF1VdBOx3Ooc7VHW3qn7vun+IyrlQap0T3RdppcOuh+GuW4N8H1khqIGIPC4iO4BfAQ87ncdNvwW+cjpEEOkI7Kj2OAc/+wIKNCKSAPQDvnM4Sr2JSKiIrAJygTmq2iDvISgLgYjMFZF1NdzGAKjqA6oaB7xH5exqPqeu9+Da5gGgjMr34XNO5j0Y4w4RaQp8Cvz+uNa+X1DVclXtS2WrPlVEGuRQXZ2T1wciVR12kpu+B8wEHmnAOKekrvcgItcBFwND1Uc7gurx/8Gf7ATiqj3u5FpmvMx1XP1T4D1V/czpPO5Q1QIRWQCMBDzeiR+ULYITEZHkag/HAJucynKqRGQk8D/AJapa5HSeILMCSBaRRBFpBIwDZjicKei4OlrfADaq6jNO5zkVIhJbdcafiERReQJCg3wf2VlDxxGRT4HuVJ6xsg24WVX96hediGQBEUC+a9EyPzzz6VLgRSAWKABWqepFjoY6SSLyC+A5IBR4U1UfdzZR/YjIB8AgKoc+3gs8oqpvOBqqnkTkPGAxsJbKv2WA+1V1pnOp6kdEzgCmUvk5CgE+UtVHG+S1rBAYY0xws0NDxhgT5KwQGGNMkLNCYIwxQc4KgTHGBDkrBMYYE+SsEBhjTJCzQmCMMUHu/wALyDaomNe2dgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "p = np.arange(-np.pi, np.pi, 0.1)\n",
    "phat_boson = 2.0*np.sin(p/2.0)\n",
    "phat_fermion = np.sin(p)\n",
    "\n",
    "plt.plot(p,phat_boson)\n",
    "plt.show()\n",
    "\n",
    "plt.plot(p,phat_fermion)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20194c53",
   "metadata": {},
   "source": [
    "We see that a naive fermion on a lattice has additional small momentum regions for $p\\approx \\pm \\pi$.\n",
    "\n",
    "One can show that the Euclidean gamma matrices satisfy $\\{ \\gamma_\\mu, \\gamma_\\nu \\} = 2 \\delta_{\\mu\\nu}$.\n",
    "Then the equations of motion yield\n",
    "$$\n",
    "\\begin{align}\n",
    "0 &= \\left(i\\gamma_\\mu \\sin(p_\\mu) - m \\right)  \\left(i\\gamma_\\mu \\sin(p_\\mu) + m \\right) \\psi(x) \\\\\n",
    "&=-\\left(\\sum_\\mu \\sin(p_\\mu)^2 + m^2 \\right) \\psi(x) \\,.\n",
    "\\end{align}\n",
    "$$\n",
    "The problem of fermion doublers arises both for the massive and massless case, however,\n",
    "in the massless case the above equations have simple real-valued solutions $p=0$ but also\n",
    "for any of the $p_\\mu=\\pi$.  In four dimensions there are therefore\n",
    "$$\n",
    "2^4 = 16\n",
    "$$\n",
    "solutions instead of one.  In the construction of the Hilbert space analog chapter 4, this would\n",
    "correspond to 15 **fermion doublers** being created in addition to the desired solution for $p=0$.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0745a6fd",
   "metadata": {},
   "source": [
    "## Wilson fermions\n",
    "\n",
    "The most straightforward solution to avoid the fermion doubler problem is to add a higher-dimensional operator that removes the spurious solutions.  The **Wilson fermion** solution is given by modifying the naive fermion action to\n",
    "$$\n",
    "\\begin{align}\n",
    "S_f = \\sum_x \\bar\\psi(x) \\left(\\sum_\\mu \\gamma_\\mu D_\\mu + m - \\frac12 \\sum_\\mu D_\\mu^2 \\right) \\psi(x) \\,.\n",
    "\\end{align}\n",
    "$$\n",
    "We define the corresponding **Dirac operator**\n",
    "$$\n",
    "\\begin{align}\n",
    "D(U) = \\sum_\\mu \\gamma_\\mu D_\\mu + m - \\frac12 \\sum_\\mu \\overleftarrow{D}_\\mu \\overrightarrow{D}_\\mu\n",
    "\\end{align}\n",
    "$$\n",
    "with\n",
    "$$\n",
    "\\begin{align}\n",
    "\\overleftarrow{D}_\\mu \\overrightarrow{D}_\\mu = C_\\mu^+ + C_\\mu^- - 2 \\,.\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "The equations of motion of this theory yield\n",
    "$$\n",
    "\\begin{align}\n",
    "0 &= \\sum_\\mu \\sin(p_\\mu)^2 + \\left(m + 4 - \\sum_\\mu\\cos(p_\\mu) \\right)^2 \\,.\n",
    "\\end{align}\n",
    "$$\n",
    "For small momenta $p$ the additional term does not contribute but it removes\n",
    "the spurious fermion doubler solutions at $p_\\mu=\\pi$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0d248a2",
   "metadata": {},
   "source": [
    "## Fermion correlation functions\n",
    "\n",
    "In a Berezin integral, it is straightforward to show that\n",
    "$$\n",
    "\\begin{align}\n",
    "\\int d\\psi f(\\psi) = \\int d\\psi f(\\psi + \\alpha)\n",
    "\\end{align}\n",
    "$$\n",
    "for a vector $\\alpha$ since the coefficient with exactly one factor $\\psi_i$\n",
    "each is identical in $f(\\psi)$ and $f(\\psi + \\alpha)$.\n",
    "Based on this observation, we consider additional source terms\n",
    "$$\n",
    "\\begin{align}\n",
    "\\int d\\bar\\psi d\\psi e^{-\\bar\\psi D \\psi - \\bar\\eta \\psi - \\bar\\psi \\eta} \n",
    "&=\\int d\\bar\\psi d\\psi e^{-\\bar\\psi D (\\psi - D^{-1}\\eta) - \\bar\\eta (\\psi - D^{-1}\\eta) - \\bar\\psi \\eta} \\\\\n",
    "&=\\int d\\bar\\psi d\\psi e^{-\\bar\\psi D \\psi + \\bar\\psi \\eta - \\bar\\eta \\psi + \\bar\\eta D^{-1}\\eta - \\bar\\psi \\eta} \\\\\n",
    "&=\\int d\\bar\\psi d\\psi e^{-\\bar\\psi D \\psi - \\bar\\eta \\psi + \\bar\\eta D^{-1}\\eta } \\\\\n",
    "&=\\int d\\bar\\psi d\\psi e^{-(\\bar\\psi - \\bar\\eta D^{-1}) D \\psi - \\bar\\eta \\psi + \\bar\\eta D^{-1}\\eta } \\\\\n",
    "&=\\int d\\bar\\psi d\\psi e^{-\\bar\\psi D \\psi + \\bar\\eta D^{-1}\\eta } \\\\\n",
    "&=\\det(D) e^{\\bar\\eta D^{-1}\\eta} \\,.\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Consider the fermion action with source terms\n",
    "$$\n",
    "\\begin{align}\n",
    "S_f(U) = \\sum_x \\bar\\psi(x)_i D_{i,x;j,y}(U) \\psi(x)_j \n",
    "\\end{align}\n",
    "$$\n",
    "and general correlation functions\n",
    "$$\n",
    "\\begin{align}\n",
    "\\langle O(U,\\psi(x)_i \\bar\\psi(y)_j) \\rangle &= \\frac1Z \\int d[U] d\\bar\\psi d\\psi O(U,\\psi(x)_i \\bar\\psi(y)_j) e^{-S_f(U)-S_g(U)} \\,, \\\\\n",
    "Z &= \\int d[U] d\\bar\\psi d\\psi e^{-S_f(U) -S_g(U)} = \\int d[U] \\det(D(U)) e^{-S_g(U)}\\,,\n",
    "\\end{align}\n",
    "$$\n",
    "where $S_g$ is the pure gauge action.  Using the additional source terms $\\bar\\eta$ and $\\eta$,\n",
    "we can write this as\n",
    "$$\n",
    "\\begin{align}\n",
    "\\langle O(U,\\psi(x)_i \\bar\\psi(y)_j ) \\rangle &= \\frac1Z \\int d[U] d\\bar\\psi d\\psi \\left(O(U,\\psi(x)_i\\partial_{\\eta(y)_j}) e^{-S_f(U)-S_g(U) - \\bar\\eta\\psi - \\bar\\psi \\eta}\\right)_{\\eta=\\bar\\eta=0} \\\\\n",
    "&= \\frac1Z \\int d[U] d\\bar\\psi d\\psi \\left(O(U,-\\partial_{\\eta(y)_j}\\psi(x)_i) e^{-S_f(U)-S_g(U) - \\bar\\eta\\psi - \\bar\\psi \\eta}\\right)_{\\eta=\\bar\\eta=0} \\\\\n",
    "&=\\frac1Z \\int d[U] d\\bar\\psi d\\psi \\left(O(U,\\partial_{\\eta(y)_j}\\partial_{\\bar\\eta(x)_i}) e^{-S_f(U)-S_g(U) - \\bar\\eta\\psi - \\bar\\psi \\eta}\\right)_{\\eta=\\bar\\eta=0} \\\\\n",
    "&=\\frac1Z \\int d[U] \\det(D(U)) e^{-S_g(U)} \\left( O(U,-\\partial_{\\bar\\eta(x)_i}\\partial_{\\eta(y)_j}) e^{\\bar\\eta D^{-1}(U) \\eta} \\right)_{\\eta=\\bar\\eta=0}\\,,\n",
    "\\end{align}\n",
    "$$\n",
    "where the derivative with respect to Grassmann variables\n",
    "are defined identical to the Berezin integrals\n",
    "$$\n",
    "\\begin{align}\n",
    "\\partial_{\\eta_a} f(\\eta_a) \\equiv \\int d\\eta_a f(\\eta_a) \\,.\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e2362d5",
   "metadata": {},
   "source": [
    "We first note that only correlation functions with matching powers of $\\psi$ and $\\bar\\psi$ can\n",
    "lead to a non-vanishing result.  All others vanish due to the $\\psi\\to-\\psi$ and $\\bar\\psi\\to-\\bar\\psi$\n",
    "symmetry.  The simplest non-vanishing correlation function is the two-point correlator (or propagator)\n",
    "$$\n",
    "\\begin{align}\n",
    "\\langle \\psi(x)_i \\bar\\psi(y)_j \\rangle &= \\langle \\left(-\\partial_{\\bar\\eta(x)_i}\\partial_{\\eta(y)_j} e^{\\bar\\eta D^{-1}(U) \\eta}\\right)_{\\eta=\\bar\\eta=0} \\rangle \\\\\n",
    "&= \\langle \\left(\\partial_{\\bar\\eta(x)_i} (\\bar\\eta D^{-1}(U))_{j,y} e^{\\bar\\eta D^{-1}(U) \\eta}\\right)_{\\eta=\\bar\\eta=0} \\rangle \\\\\n",
    "&= \\langle D^{-1}(U)_{i,x;j,y} \\rangle \\,.\n",
    "\\end{align}\n",
    "$$\n",
    "The calculation goes through in an analogous manner for more fermionic fields, however, one needs\n",
    "to pay attention to permutations of fields.  In general, we find **Wick's theorem**\n",
    "$$\n",
    "\\begin{align}\n",
    "\\langle \\psi(x_1)_{i_1} \\bar\\psi(y_1)_{j_1} \\cdots \\psi(x_n)_{i_n} \\bar\\psi(y_n)_{j_n} \\rangle\n",
    "&=\\sum_P (-1)^P \\langle D^{-1}(U)_{i_1,x_1;j_{P_1},y_{P_1}} \\cdots D^{-1}(U)_{i_n,x_n;j_{P_n},y_{P_n}} \\rangle  \\,,\n",
    "\\end{align}\n",
    "$$\n",
    "where the sum is over all permutations $(P_1,P_2,\\ldots,P_n)$ of $(1,2,\\ldots, n)$ with $(-1)^P$ denoting the sign of the permutations (1 if even numbers of swaps are needed to recover $(1,2,\\ldots, n)$, -1 else).  For concreteness,\n",
    "consider the next non-trivial example\n",
    "$$\n",
    "\\begin{align}\n",
    "\\langle \\psi(x_1)_{i_1} \\bar\\psi(y_1)_{j_1} \\psi(x_2)_{i_2} \\bar\\psi(y_2)_{j_2} \\rangle\n",
    "&= \\langle D^{-1}(U)_{i_1,x_1;j_{1},y_{1}} D^{-1}(U)_{i_2,x_2;j_{2},y_{2}} \\rangle - \\langle D^{-1}(U)_{i_1,x_1;j_{2},y_{2}} D^{-1}(U)_{i_2,x_2;j_{1},y_{1}} \\rangle \\,.\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f473160",
   "metadata": {},
   "source": [
    "We note that such correlation functions can be mapped to expectation values of **time-ordered products** of operators $\\hat{\\Psi}$ and $\\hat{\\bar\\Psi}$ acting on the Hilbert space by using the path-integral representation of $e^{-\\delta \\tau H}$.  This is similar to all studied cases so far.  In the Schrödinger picture, we therefore have, e.g.,\n",
    "$$\n",
    "\\begin{align}\n",
    "\\langle \\bar\\psi(t) \\Gamma_t \\psi(t) \\bar\\psi(t=0) \\Gamma_0 \\psi(t=0) \\rangle\n",
    "&=\\frac1Z \\textrm{Tr}\\left[ e^{-(T-t) H} \\hat{\\bar\\Psi} \\Gamma_t \\hat{\\Psi} e^{- t H} \\hat{\\bar\\Psi} \\Gamma_0 \\hat{\\Psi}\\right] \\,.\n",
    "\\end{align}\n",
    "$$\n",
    "We can then again consider inserting complete sets of states and studying the $t$-dependence of the correlation function to learn about the fermionic spectrum and matrix elements of the theory."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb6c541c",
   "metadata": {},
   "source": [
    "## Anti-periodic boundary conditions in time\n",
    "\n",
    "When implementing the trace in the Hilbert space in the above formula, one surprisingly has to chose **anti-periodic** boundary conditions in time for the fermionic fields.  To see this and also get more comfortable with the mapping from the Hilbert space to our extended algebra, we review the case of a single fermionic degree of freedom in detail.\n",
    "\n",
    "First, consider the states\n",
    "$$\n",
    "\\begin{align}\n",
    "\\vert q \\rangle &=  \\left(1- i \\hat{p} q \\right) \\vert 0 \\rangle = \\vert 0 \\rangle + i q \\vert 1 \\rangle \\\\\n",
    "\\langle q' \\vert  &= \\langle 0 \\vert \\hat{q} \\left(1- i  q' \\hat{p}\\right) \\\\\n",
    "&=\\langle 0 \\vert (\\hat{q} + i  q' \\hat{q}\\hat{p}) \\\\\n",
    "&=\\langle 0 \\vert (\\hat{q} - q') \\\\\n",
    "&= \\langle 1 \\vert i - \\langle 0 \\vert q' \\,,\n",
    "\\end{align}\n",
    "$$\n",
    "where we identify\n",
    "$$\n",
    "\\begin{align}\n",
    "\\vert 1 \\rangle &= \\hat{p} \\vert 0 \\rangle \\,, \\\\\n",
    "\\langle 1 \\vert &= \\langle 0 \\vert (-i\\hat{q})\n",
    "\\end{align}\n",
    "$$\n",
    "for which we have showed above that\n",
    "$$\n",
    "\\begin{align}\n",
    "\\langle 0  \\vert 0 \\rangle &= \\langle 1 \\vert 1 \\rangle = 1 \\,, \\\\\n",
    "\\langle 0 \\vert 1 \\rangle &= \\langle 1 \\vert 0 \\rangle = 0 \\,.\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "It is, however, crucial to remember that we also imposed anticommutation relations between $\\hat{q}$ and $q'$ such that the $\\vert 1 \\rangle$ and $\\langle 1 \\vert$ states anticommute with $q$ and $q'$.  It is difficult to keep track of this and therefore we often wait until the end before we express $\\langle 0 \\vert (-i\\hat{q}) = \\langle 1 \\vert$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c0054906",
   "metadata": {},
   "source": [
    "Let us remind ourselves of the completeness relation in this simple example.\n",
    "We find\n",
    "$$\n",
    "\\begin{align}\n",
    "\\vert q \\rangle \\langle q \\vert &= ( \\vert 0 \\rangle + i q \\vert 1 \\rangle)( \\langle 1 \\vert i - \\langle 0 \\vert q ) \\\\\n",
    "&=i \\vert 0 \\rangle \\langle 1 \\vert - q \\vert 0 \\rangle \\langle 0 \\vert -  q \\vert 1 \\rangle \\langle 1 \\vert \\\\\n",
    "\\int (-dq) \\vert q \\rangle \\langle q \\vert &= \\vert 0 \\rangle \\langle 0 \\vert+ \\vert 1 \\rangle \\langle 1 \\vert\n",
    "\\end{align}\n",
    "$$\n",
    "as we desired.\n",
    "\n",
    "Finally, consider the trace of a matrix.  We find\n",
    "$$\n",
    "\\begin{align}\n",
    "\\int dq \\langle -q \\vert A \\vert q \\rangle &= \\int dq \\langle 0 \\vert (\\hat{q} + q) A (\\vert 0 \\rangle + i q \\vert 1 \\rangle) \\\\\n",
    "&= \\int dq \\langle 0 \\vert q A \\vert 0 \\rangle + \\int dq \\langle 0 \\vert \\hat{q} A i q \\vert 1 \\rangle \\\\\n",
    "&= \\int dq q \\langle 0 \\vert A \\vert 0 \\rangle +\\int dq q\\langle 0 \\vert (- i\\hat{q}) A \\vert 1 \\rangle \\\\\n",
    "&= \\langle 0 \\vert A \\vert 0 \\rangle +\\langle 1 \\vert A \\vert 1 \\rangle = \\textrm{Tr} A \\,.\n",
    "\\end{align}\n",
    "$$\n",
    "So we find indeed that antiperiodic boundary conditions have to be used to implement the trace over a fermionic Hilbertspace operator."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b14847ec",
   "metadata": {},
   "source": [
    "## Quenched theory\n",
    "\n",
    "We will soon study an algorithm to efficiently generate gauge configurations for theories with both fermions and gauge degrees of freedom.  Before doing this, however, we can also study the **quenched** approximation in which we neglect the fermion determinant in the path integral and instead measure fermionic correlation functions in gauge configurations generated using only the gluonic action.\n",
    "\n",
    "## Pion correlator\n",
    "\n",
    "In the next chapter, we study the symmetries of the action in more details and will identify good quantum numbers of our Hamiltonian.  To conclude this chapter, however, we will now study for the first time the lightest particle made out of a quark and an antiquark that can be created using the operator\n",
    "$$\n",
    "\\begin{align}\n",
    "\\hat{O}_\\pi(t) = i\\sum_\\vec{x} \\hat{\\bar{u}}(\\vec{x},t) \\gamma_5 \\hat{d}(\\vec{x},t)\n",
    "\\end{align}\n",
    "$$\n",
    "with $\\gamma_5 =  \\gamma_1 \\gamma_2 \\gamma_3 \\gamma_4$\n",
    "where $u$ and $d$ two different quenched fermion degrees of freedom.\n",
    "Note that from the definition of the conjugate momentum we require\n",
    "$\\hat{\\bar{u}} = \\hat{u}^\\dagger \\gamma_4$ such that\n",
    "$$\n",
    "\\begin{align}\n",
    "\\hat{O}_\\pi(t)^\\dagger = \\sum_\\vec{x}(i\\hat{u}(\\vec{x},t)^\\dagger \\gamma_4 \\gamma_5 \\hat{d}(\\vec{x},t))^\\dagger = -\\sum_\\vec{x}i\\hat{d}^\\dagger(\\vec{x},t) \\gamma_5 \\gamma_4 \\hat{u}(\\vec{x},t)\n",
    "=\\sum_\\vec{x}i\\hat{\\bar{d}}(\\vec{x},t) \\gamma_5 \\hat{u}(\\vec{x},t)\n",
    "\\end{align}\n",
    "$$\n",
    "since $\\gamma_\\mu^\\dagger = \\gamma_\\mu$.\n",
    "\n",
    "Applying the Wick theorem, we find\n",
    "$$\n",
    "\\begin{align}\n",
    "C(t) &= \\langle \\hat{O}_\\pi(t) \\hat{O}^\\dagger_\\pi(0) \\rangle \\\\\n",
    "&= \\sum_{\\vec{x},\\vec{y}}\\langle \\textrm{Tr}[ D^{-1}(U)_{\\vec{x},t; \\vec{y},0} \\gamma_5 D^{-1}(U)_{\\vec{y},0; \\vec{x},t} \\gamma_5] \\rangle \\,.\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Next, one uses that the Dirac operator is $\\gamma_5$-Hermitian, i.e.,\n",
    "$$\n",
    "\\begin{align}\n",
    "(\\gamma_5 D(U))^\\dagger = \\gamma_5 D(U)\n",
    "\\end{align}\n",
    "$$\n",
    "which combined with $\\gamma_5^\\dagger=\\gamma_5$ and $\\gamma_5^2=\\mathbb{1}$ yields\n",
    "$$\n",
    "\\begin{align}\n",
    "D(U)^\\dagger = \\gamma_5 D(U) \\gamma_5\n",
    "\\end{align}\n",
    "$$\n",
    "and therefore\n",
    "$$\n",
    "\\begin{align}\n",
    "C(t)\n",
    "&= \\sum_{\\vec{x},\\vec{y}}\\langle \\textrm{Tr}[ D^{-1}(U)_{\\vec{x},t; \\vec{y},0}  (D^{-1}(U)_{\\vec{x},t; \\vec{y},0})^\\dagger] \\rangle \\,.\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Finally, we can use translational symmetry to write\n",
    "$$\n",
    "\\begin{align}\n",
    "C(t)\n",
    "&= V\\sum_{\\vec{x}}\\langle \\textrm{Tr}[ D^{-1}(U)_{\\vec{x},t; \\vec{0},0}  (D^{-1}(U)_{\\vec{x},t; \\vec{0},0})^\\dagger] \\rangle \\,.\n",
    "\\end{align}\n",
    "$$\n",
    "with spatial volume $V$.  If we therefore found the vector\n",
    "$$\n",
    "\\begin{align}\n",
    "\\alpha_{\\vec{x},t} = D^{-1}(U)_{\\vec{x},t; \\vec{0},0}\n",
    "\\end{align}\n",
    "$$\n",
    "for a given gauge-field configuration, we could create the correlator $C(t)$ in a straightforward way.\n",
    "We conclude this chapter by demonstrating code that finds this vector for a given gauge field configuration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "e9454285",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GPT :  203973.351320 s : Reading 8c32_5.7/su3.200\n",
      "GPT :  203973.353464 s : Switching view to [1,1,1,1]/Read\n",
      "GPT :  203973.369661 s : Read 0.00219727 GB at 0.13545 GB/s (0.165431 GB/s for distribution, 0.748903 GB/s for reading + checksum, 5.95349 GB/s for checksum, 1 views per node)\n",
      "GPT :  203973.373617 s : Read 0.00219727 GB at 0.772312 GB/s (2.66976 GB/s for distribution, 1.08936 GB/s for reading + checksum, 5.82922 GB/s for checksum, 1 views per node)\n",
      "GPT :  203973.376494 s : Read 0.00219727 GB at 1.12103 GB/s (3.78792 GB/s for distribution, 1.59806 GB/s for reading + checksum, 6.46283 GB/s for checksum, 1 views per node)\n",
      "GPT :  203973.379951 s : Read 0.00219727 GB at 0.797025 GB/s (1.93939 GB/s for distribution, 1.35729 GB/s for reading + checksum, 6.8065 GB/s for checksum, 1 views per node)\n",
      "GPT :  203973.380776 s : Completed reading 8c32_5.7/su3.200 in 0.029685 s\n"
     ]
    }
   ],
   "source": [
    "U = g.load(\"8c32_5.7/su3.200\")\n",
    "grid = U[0].grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 158,
   "id": "4d85d44c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GPT :  207178.063015 s : Initializing gpt.random(13,vectorized_ranlux24_389_64) took 0.000413179 s\n",
      "(678.8217289003439+1316.0853819848078j)\n",
      "(678.8217289003376+1316.0853819848253j)\n",
      "(2760.203862794375+1.8189894035458565e-12j)\n"
     ]
    }
   ],
   "source": [
    "Uapt = g.copy(U)\n",
    "# anti-periodic fermion in time can be achieved modifying last link in time\n",
    "Uapt[3][:,:,:,grid.gdimensions[3]-1] *= -1\n",
    "\n",
    "def DiracWilson(U, src, m):\n",
    "    dst = g((4.0 + m) * src)\n",
    "    for mu in range(4):\n",
    "        fmu = g(U[mu] * g.cshift(src, mu, 1))\n",
    "        bmu = g(g.cshift(g.adj(U[mu]) * src, mu, -1))\n",
    "        dst += 0.5*(g.gamma[mu] * fmu - fmu)\n",
    "        dst -= 0.5*(g.gamma[mu] * bmu + bmu)\n",
    "    return dst\n",
    "\n",
    "rng = g.random(\"13\")\n",
    "a, b = rng.cnormal([g.vspincolor(grid),g.vspincolor(grid)])\n",
    "\n",
    "# test gamma5 hermiticity\n",
    "print(g.inner_product(b, g.gamma[5] * DiracWilson(Uapt, a, -0.5)))\n",
    "print(g.inner_product(a, g.gamma[5] * DiracWilson(Uapt, b, -0.5)).conjugate())\n",
    "\n",
    "# norm real\n",
    "print(g.inner_product(a, g.gamma[5] * DiracWilson(Uapt, a, -0.5)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 159,
   "id": "3c6c91f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "g.default.set_verbose(\"fgmres_convergence\", False)\n",
    "inv = g.algorithms.inverter.fgmres(eps = 1e-8, maxiter = 400)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 160,
   "id": "e8591111",
   "metadata": {},
   "outputs": [],
   "source": [
    "def G5DiracWilson(dst, src):\n",
    "    dst @= g.gamma[5] * DiracWilson(Uapt, src, -0.5)\n",
    "    return dst\n",
    "\n",
    "invG5D = inv(G5DiracWilson)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 161,
   "id": "00409fec",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "lattice(ot_matrix_spin_color(4,3),double)"
      ]
     },
     "execution_count": 161,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "src = g.mspincolor(grid)\n",
    "g.create.point(src, [0, 0, 0, 0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 162,
   "id": "15d1777b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GPT :  207288.383994 s : fgmres: NOT converged in 400 iterations;  computed squared residual 1.251610e-12 / 1.200000e-15;  true squared residual 1.251610e-12 / 1.200000e-15\n"
     ]
    }
   ],
   "source": [
    "dst = g(invG5D * g.gamma[5] * src)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7b01cee",
   "metadata": {},
   "source": [
    "This is a very inefficient way to invert this matrix but it is easy to implement.  We will learn better methods soon."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 166,
   "id": "60bdfe2b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Error of inversion: 1.25160956143534e-12\n"
     ]
    }
   ],
   "source": [
    "print(\"Error of inversion:\", g.norm2(DiracWilson(Uapt, dst, -0.5) - src))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 164,
   "id": "ad0b9891",
   "metadata": {},
   "outputs": [],
   "source": [
    "Ct = g.slice(g.trace(g.adj(dst)*dst),3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6e86492",
   "metadata": {},
   "source": [
    "Let us now study the effective mass of this particle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 165,
   "id": "4b979d0b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAUyElEQVR4nO3df5Bd5X3f8fenQrZl7LHkoLigH0BTBuwYG9k72CmZGscBZIcCdTNTUcfBqT2a8RjHSVtaaGbCDJ6paenkNw7W2Cp260BSDERNbcsa45Q2Da5WQPgZGRX/QGtSKRZgJ2gMkr/9Y4/sy7Kre7V7d+/uo/drZmfveZ5z7v2eRXzuuc957jmpKiRJ7fo7oy5AkjS/DHpJapxBL0mNM+glqXEGvSQ1zqCXpMb1Dfok65J8JckjSR5O8pFp1kmS30myJ8kDSd7U03dFkse6nyuGvQOSpKNLv3n0SU4GTq6qe5O8EtgFXFZVj/Ss8y7gw8C7gLcAv11Vb0nyamAcGAOq2/bNVfXUvOyNJOlF+h7RV9WTVXVv9/h7wKPAmimrXQp8pibdA6zs3iAuAnZU1YEu3HcAG4e6B5KkozrhWFZOchqwAfjqlK41wBM9y3u7tpnap3vuzcBmgBNPPPHNZ5111rGUJknHtV27dv11Va2erm/goE/yCuBzwK9U1XeHVdwRVbUF2AIwNjZW4+Pjw34JSWpWkm/O1DfQrJsky5kM+c9W1e3TrDIBrOtZXtu1zdQuSVogg8y6CfAp4NGq+o0ZVtsG/GI3++atwDNV9SSwHbgwyaokq4ALuzZJ0gIZZOjmPOC9wINJ7u/a/i2wHqCqbgI+z+SMmz3As8AvdX0HknwU2Nltd11VHRha9ZKkvvoGfVX9LyB91ingQzP0bQW2zqo6SdKc+c1YSWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpccd045Hj3Z33TXDD9t18++mDnLJyBVdddCaXbZj2PiqStGgY9AO6874Jrrn9QQ4+fxiAiacPcs3tDwIY9pIWNYduBnTD9t0/DPkjDj5/mBu27x5RRZI0GIN+QN9++uAxtUvSYmHQD+iUlSuOqV2SFguDfkBXXXQmK5Yve0HbiuXLuOqiM0dUkSQNxpOxAzpywtVZN5KWmr5Bn2QrcDGwr6peP03/VcB7ep7vtcDq7n6x3wC+BxwGDlXV2LAKH4XLNqwx2CUtOYMM3dwMbJyps6puqKpzquoc4Brgf0y5Afjbu/4lHfKStFT1Dfqquhs40G+9zuXALXOqSJI0VEM7GZvk5Uwe+X+up7mALyXZlWTzsF5LkjS4YZ6M/UfAn00ZtvnpqppI8uPAjiR/2X1CeJHujWAzwPr164dYliQd34Y5vXITU4Ztqmqi+70PuAM4d6aNq2pLVY1V1djq1auHWJYkHd+GEvRJXgW8DfjjnrYTk7zyyGPgQuChYbyeJGlwg0yvvAU4HzgpyV7gWmA5QFXd1K32j4EvVdXf9mz6GuCOJEde5w+q6ovDK12SNIi+QV9Vlw+wzs1MTsPsbXsceONsC5MkDYeXQJCkxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMZ5K8FF6M77JrxloaShMegXmTvvm+Ca2x/k4POHAZh4+iDX3P4ggGEvaVYcullkbti++4chf8TB5w9zw/bdI6pI0lJn0C8y33764DG1S1I/Dt0sMqesXMHENKF+ysoVQ30dzwNIxw+P6BeZqy46kxXLl72gbcXyZVx10ZlDe40j5wEmnj5I8aPzAHfeNzG015C0eBj0i8xlG9bwsXefzZqVKwiwZuUKPvbus4d6tO15AOn44tDNInTZhjXzOozieQDp+OIR/XFopvH+YZ8HkLQ4GPTHoYU4DyBp8egb9Em2JtmX5KEZ+s9P8kyS+7ufX+/p25hkd5I9Sa4eZuGavYU4DyBp8RhkjP5m4PeAzxxlnf9ZVRf3NiRZBtwIXADsBXYm2VZVj8yyVg3RfJ8HAKdwSotF3yP6qrobODCL5z4X2FNVj1fVc8CtwKWzeB4tQU7hlBaPYY3R/1SSv0jyhSQ/2bWtAZ7oWWdv1zatJJuTjCcZ379//5DK0qg4hVNaPIYR9PcCp1bVG4HfBe6czZNU1ZaqGquqsdWrVw+hLI2SUzilxWPOQV9V362qv+kefx5YnuQkYAJY17Pq2q5NxwGncEqLx5yDPsnfTZLu8bndc34H2AmckeT0JC8BNgHb5vp6WhoWcgrnnfdNcN71d3H61f+d866/y/MA0hR9Z90kuQU4HzgpyV7gWmA5QFXdBPw88MEkh4CDwKaqKuBQkiuB7cAyYGtVPTwve6FF58jsmvmedeP1+6X+MpnJi8vY2FiNj4+PugwtAeddf9e0V/tcs3IFf3b1z4ygImk0kuyqqrHp+vxmrJY0T/pK/Rn0WtI86Sv1Z9BrSfO6PVJ/XqZYS9pCnfSVljKDXkveQly3R1rKHLqRpMYZ9JLUOINekhpn0EtS4zwZKw1gIW6i4o1aNF8MeqmPhbiejtfs0Xxy6EbqYyFuotLajVq8ouji4hG91MdCXE+npWv2+Olk8fGIXupjIa6n09I1e1r7dNICg17qYyGup9PSjVpa+nQCbQxDOXQj9bEQ19Np6UYtp6xcMe09Aob96WShZkK1MAzljUek48hC3KhlajjC5KeTj7377HmbpTQfrwELd2ObYbxpeeMRScDCDKtctmENH3v32axZuYIwGYrDDuCFOg+wEH+vI29aE08fpPjRp4ZhDhENcs/YrcDFwL6qev00/e8B/g0Q4HvAB6vqL7q+b3Rth4FDM73bSFoYCzWsMt9XFF2o8wAL8fc62pvWsP6GgxzR3wxsPEr/14G3VdXZwEeBLVP6315V5xjy0ui1cqOWhZqltBB/r4V40+ob9FV1N3DgKP3/u6qe6hbvAdYOqTZJQ7YQwyoLYaHesBbi77UQb1oDnYxNchrwJ9MN3UxZ718BZ1XVB7rlrwNPAQV8oqqmHu33brsZ2Aywfv36N3/zm98cdB8kHYdauTbQsE4sH+1k7NCCPsnbgY8DP11V3+na1lTVRJIfB3YAH+4+IRyVs24kHU/me9bNUObRJ3kD8EngnUdCHqCqJrrf+5LcAZwL9A16STqezPfJ6zlPr0yyHrgdeG9Vfa2n/cQkrzzyGLgQeGiurydJOjaDTK+8BTgfOCnJXuBaYDlAVd0E/DrwY8DHk8CPplG+BrijazsB+IOq+uI87IMk6Sj6Bn1VXd6n/wPAB6Zpfxx44+xLkyQNg9+MlaTGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUuIGCPsnWJPuSTHtz70z6nSR7kjyQ5E09fVckeaz7uWJYhUuSBjPoEf3NwMaj9L8TOKP72Qz8PkCSVzN5M/G3AOcC1yZZNdtiJUnHbqCgr6q7gQNHWeVS4DM16R5gZZKTgYuAHVV1oKqeAnZw9DcMSdKQDWuMfg3wRM/y3q5tpvYXSbI5yXiS8f379w+pLEnSojkZW1VbqmqsqsZWr1496nIkqRnDCvoJYF3P8tqubaZ2SdICGVbQbwN+sZt981bgmap6EtgOXJhkVXcS9sKuTZK0QE4YZKUktwDnAycl2cvkTJrlAFV1E/B54F3AHuBZ4Je6vgNJPgrs7J7quqo62kldSdKQDRT0VXV5n/4CPjRD31Zg67GXJkkahkVzMlaSND8MeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWrcQEGfZGOS3Un2JLl6mv7fTHJ/9/O1JE/39B3u6ds2xNolSQPoe8/YJMuAG4ELgL3AziTbquqRI+tU1a/2rP9hYEPPUxysqnOGVrEk6ZgMckR/LrCnqh6vqueAW4FLj7L+5cAtwyhOkjR3gwT9GuCJnuW9XduLJDkVOB24q6f5ZUnGk9yT5LKZXiTJ5m698f379w9QliRpEMM+GbsJuK2qDve0nVpVY8A/A34ryU9Mt2FVbamqsaoaW7169ZDLkqTj1yBBPwGs61le27VNZxNThm2qaqL7/Tjwp7xw/F6SNM8GCfqdwBlJTk/yEibD/EWzZ5KcBawC/rynbVWSl3aPTwLOAx6Zuq0kaf70nXVTVYeSXAlsB5YBW6vq4STXAeNVdST0NwG3VlX1bP5a4BNJfsDkm8r1vbN1JEnzLy/M5cVhbGysxsfHR12GJC0ZSXZ150NfxG/GSlLjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklq3EBBn2Rjkt1J9iS5epr+9yXZn+T+7ucDPX1XJHms+7limMVLkvrre3PwJMuAG4ELgL3AziTbprnJ9x9W1ZVTtn01cC0wBhSwq9v2qaFUL0nqa5Aj+nOBPVX1eFU9B9wKXDrg818E7KiqA1247wA2zq5USdJsDBL0a4Anepb3dm1T/ZMkDyS5Lcm6Y9yWJJuTjCcZ379//wBlSZIGMayTsf8NOK2q3sDkUfunj/UJqmpLVY1V1djq1auHVJYkaZCgnwDW9Syv7dp+qKq+U1Xf7xY/Cbx50G0lSfNrkKDfCZyR5PQkLwE2Adt6V0hycs/iJcCj3ePtwIVJViVZBVzYtUmSFkjfWTdVdSjJlUwG9DJga1U9nOQ6YLyqtgG/nOQS4BBwAHhft+2BJB9l8s0C4LqqOjAP+yFJmkGqatQ1vMjY2FiNj4+PugxJWjKS7Kqqsen6/GasJDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGDRT0STYm2Z1kT5Krp+n/F0keSfJAki8nObWn73CS+7ufbVO3lSTNr743B0+yDLgRuADYC+xMsq2qHulZ7T5grKqeTfJB4D8A/7TrO1hV5wy3bEnSoAY5oj8X2FNVj1fVc8CtwKW9K1TVV6rq2W7xHmDtcMuUJM3WIEG/BniiZ3lv1zaT9wNf6Fl+WZLxJPckuezYS5QkzUXfoZtjkeQXgDHgbT3Np1bVRJK/B9yV5MGq+r/TbLsZ2Aywfv36YZYlSce1QY7oJ4B1Pctru7YXSPKzwK8Bl1TV94+0V9VE9/tx4E+BDdO9SFVtqaqxqhpbvXr1wDsgSTq6QYJ+J3BGktOTvATYBLxg9kySDcAnmAz5fT3tq5K8tHt8EnAe0HsSV5I0z/oO3VTVoSRXAtuBZcDWqno4yXXAeFVtA24AXgH81yQA36qqS4DXAp9I8gMm31SunzJbR5I0z1JVo67hRcbGxmp8fHzUZUjSkpFkV1WNTdfnN2MlqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxg0U9Ek2JtmdZE+Sq6fpf2mSP+z6v5rktJ6+a7r23UkuGmLtkqQB9A36JMuAG4F3Aq8DLk/yuimrvR94qqr+PvCbwL/vtn0dsAn4SWAj8PHu+SRJC2SQI/pzgT1V9XhVPQfcClw6ZZ1LgU93j28D3pEkXfutVfX9qvo6sKd7PknSAjlhgHXWAE/0LO8F3jLTOlV1KMkzwI917fdM2XbNdC+SZDOwuVv8myS7B6htOicBfz3LbRebVvallf0A92UxamU/YG77cupMHYME/YKoqi3Alrk+T5LxqhobQkkj18q+tLIf4L4sRq3sB8zfvgwydDMBrOtZXtu1TbtOkhOAVwHfGXBbSdI8GiTodwJnJDk9yUuYPLm6bco624Arusc/D9xVVdW1b+pm5ZwOnAH8n+GULkkaRN+hm27M/UpgO7AM2FpVDye5Dhivqm3Ap4D/nGQPcIDJNwO69f4IeAQ4BHyoqg7P074cMefhn0WklX1pZT/AfVmMWtkPmKd9yeSBtySpVX4zVpIaZ9BLUuOaCfp+l2lYKpKsS/KVJI8keTjJR0Zd01wlWZbkviR/Mupa5iLJyiS3JfnLJI8m+alR1zQbSX61+7f1UJJbkrxs1DUNKsnWJPuSPNTT9uokO5I81v1eNcoaBzXDvtzQ/ft6IMkdSVYO47WaCPoBL9OwVBwC/mVVvQ54K/ChJbwvR3wEeHTURQzBbwNfrKqzgDeyBPcpyRrgl4Gxqno9kxMsNo22qmNyM5OXU+l1NfDlqjoD+HK3vBTczIv3ZQfw+qp6A/A14JphvFATQc9gl2lYEqrqyaq6t3v8PSbDZNpvEy8FSdYCPwd8ctS1zEWSVwH/kMkZZlTVc1X19EiLmr0TgBXdd15eDnx7xPUMrKruZnJmX6/eS7B8GrhsIWuaren2paq+VFWHusV7mPzu0Zy1EvTTXaZhyYbjEd1VQDcAXx1xKXPxW8C/Bn4w4jrm6nRgP/CfumGoTyY5cdRFHauqmgD+I/At4Engmar60mirmrPXVNWT3eO/Al4zymKG6J8DXxjGE7US9M1J8grgc8CvVNV3R13PbCS5GNhXVbtGXcsQnAC8Cfj9qtoA/C1LZ4jgh7rx60uZfOM6BTgxyS+Mtqrh6b6oueTnjCf5NSaHcT87jOdrJeibutRCkuVMhvxnq+r2UdczB+cBlyT5BpPDaT+T5L+MtqRZ2wvsraojn65uYzL4l5qfBb5eVfur6nngduAfjLimufp/SU4G6H7vG3E9c5LkfcDFwHtqSF90aiXoB7lMw5LQXd75U8CjVfUbo65nLqrqmqpaW1WnMfnf5K6qWpJHj1X1V8ATSc7smt7B5De+l5pvAW9N8vLu39o7WIInlafovQTLFcAfj7CWOUmykcmhzkuq6tlhPW8TQd+dvDhymYZHgT+qqodHW9WsnQe8l8mj3/u7n3eNuigB8GHgs0keAM4B/t1oyzl23SeS24B7gQeZzIAlcwmBJLcAfw6cmWRvkvcD1wMXJHmMyU8s14+yxkHNsC+/B7wS2NH9v3/TUF7LSyBIUtuaOKKXJM3MoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mN+/+RnsNWFWGbDAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.ylim(0,2)\n",
    "plt.plot(range(13),[np.log(Ct[t] / Ct[t+1]).real for t in range(13)], marker=\"o\", ls=\"\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd8c0480",
   "metadata": {},
   "source": [
    "This is the effective mass of a pion propagating on a single quenched gauge configuration.  We note that the particle is quite heavy.  Next, we will study better methods for the inversion, which will allow us to study lighter masses as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 181,
   "id": "a8dba566",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test D: 1.7159776039272344e-25\n",
      "GPT :  213012.700159 s : cg: converged in 40 iterations\n",
      "GPT :  213012.930778 s : cg: converged in 40 iterations\n",
      "GPT :  213013.152907 s : cg: converged in 40 iterations\n",
      "GPT :  213013.362946 s : cg: converged in 40 iterations\n",
      "GPT :  213013.594801 s : cg: converged in 40 iterations\n",
      "GPT :  213013.855515 s : cg: converged in 40 iterations\n",
      "GPT :  213014.080708 s : cg: converged in 40 iterations\n",
      "GPT :  213014.305992 s : cg: converged in 40 iterations\n",
      "GPT :  213014.499045 s : cg: converged in 38 iterations\n",
      "GPT :  213014.728914 s : cg: converged in 39 iterations\n",
      "GPT :  213014.973023 s : cg: converged in 39 iterations\n",
      "GPT :  213015.204403 s : cg: converged in 39 iterations\n",
      "Test D^{-1}: 1.3723199796550752e-11\n"
     ]
    }
   ],
   "source": [
    "w = g.qcd.fermion.wilson_clover(U, mass=-0.5, csw_r=0.0, csw_t=0.0, nu=1.0, xi_0=1.0,\n",
    "                                isAnisotropic=False,\n",
    "                                boundary_phases=[1,1,1,-1])\n",
    "\n",
    "print(\"Test D:\",g.norm2(w(b) - DiracWilson(Uapt, b, -0.5)))\n",
    "\n",
    "inv = g.algorithms.inverter\n",
    "pc = g.qcd.fermion.preconditioner\n",
    "g.default.set_verbose(\"cg_convergence\", False)\n",
    "cg = inv.cg({\"eps\": 1e-6, \"maxiter\": 1000})\n",
    "invD = w.propagator(inv.preconditioned(pc.eo1_ne(), cg))\n",
    "\n",
    "dst2 = g(invD * src)\n",
    "\n",
    "print(\"Test D^{-1}:\",g.norm2(dst-dst2) / g.norm2(dst2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 184,
   "id": "2de31a0d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GPT :  213203.588247 s : cg: converged in 40 iterations\n",
      "GPT :  213203.806645 s : cg: converged in 40 iterations\n",
      "GPT :  213204.011701 s : cg: converged in 40 iterations\n",
      "GPT :  213204.211629 s : cg: converged in 40 iterations\n",
      "GPT :  213204.415629 s : cg: converged in 40 iterations\n",
      "GPT :  213204.635100 s : cg: converged in 40 iterations\n",
      "GPT :  213204.840919 s : cg: converged in 40 iterations\n",
      "GPT :  213205.039177 s : cg: converged in 40 iterations\n",
      "GPT :  213205.233666 s : cg: converged in 38 iterations\n",
      "GPT :  213205.428045 s : cg: converged in 39 iterations\n",
      "GPT :  213205.623858 s : cg: converged in 39 iterations\n",
      "GPT :  213205.833268 s : cg: converged in 39 iterations\n",
      "GPT :  213206.770324 s : cg: converged in 48 iterations\n",
      "GPT :  213207.058377 s : cg: converged in 48 iterations\n",
      "GPT :  213207.309152 s : cg: converged in 48 iterations\n",
      "GPT :  213207.542193 s : cg: converged in 48 iterations\n",
      "GPT :  213207.819426 s : cg: converged in 48 iterations\n",
      "GPT :  213208.113768 s : cg: converged in 48 iterations\n",
      "GPT :  213208.385462 s : cg: converged in 48 iterations\n",
      "GPT :  213208.657899 s : cg: converged in 49 iterations\n",
      "GPT :  213208.907246 s : cg: converged in 46 iterations\n",
      "GPT :  213209.160413 s : cg: converged in 48 iterations\n",
      "GPT :  213209.413827 s : cg: converged in 47 iterations\n",
      "GPT :  213209.669169 s : cg: converged in 48 iterations\n",
      "GPT :  213210.652197 s : cg: converged in 61 iterations\n",
      "GPT :  213210.991799 s : cg: converged in 61 iterations\n",
      "GPT :  213211.347376 s : cg: converged in 62 iterations\n",
      "GPT :  213211.676221 s : cg: converged in 62 iterations\n",
      "GPT :  213212.008118 s : cg: converged in 62 iterations\n",
      "GPT :  213212.328634 s : cg: converged in 61 iterations\n",
      "GPT :  213212.644664 s : cg: converged in 62 iterations\n",
      "GPT :  213212.981035 s : cg: converged in 62 iterations\n",
      "GPT :  213213.366373 s : cg: converged in 59 iterations\n",
      "GPT :  213213.702463 s : cg: converged in 61 iterations\n",
      "GPT :  213214.035591 s : cg: converged in 60 iterations\n",
      "GPT :  213214.363354 s : cg: converged in 62 iterations\n",
      "GPT :  213215.353759 s : cg: converged in 85 iterations\n",
      "GPT :  213215.830689 s : cg: converged in 85 iterations\n",
      "GPT :  213216.308263 s : cg: converged in 86 iterations\n",
      "GPT :  213216.808890 s : cg: converged in 86 iterations\n",
      "GPT :  213217.374935 s : cg: converged in 85 iterations\n",
      "GPT :  213217.925377 s : cg: converged in 84 iterations\n",
      "GPT :  213218.441550 s : cg: converged in 85 iterations\n",
      "GPT :  213218.980990 s : cg: converged in 85 iterations\n",
      "GPT :  213219.409078 s : cg: converged in 81 iterations\n",
      "GPT :  213219.918039 s : cg: converged in 84 iterations\n",
      "GPT :  213220.468912 s : cg: converged in 83 iterations\n",
      "GPT :  213221.040040 s : cg: converged in 85 iterations\n",
      "GPT :  213222.792975 s : cg: converged in 135 iterations\n",
      "GPT :  213223.667487 s : cg: converged in 136 iterations\n",
      "GPT :  213224.576976 s : cg: converged in 137 iterations\n",
      "GPT :  213225.565169 s : cg: converged in 138 iterations\n",
      "GPT :  213226.570092 s : cg: converged in 137 iterations\n",
      "GPT :  213227.423486 s : cg: converged in 135 iterations\n",
      "GPT :  213228.310716 s : cg: converged in 139 iterations\n",
      "GPT :  213229.255515 s : cg: converged in 140 iterations\n",
      "GPT :  213230.114226 s : cg: converged in 129 iterations\n",
      "GPT :  213230.983376 s : cg: converged in 137 iterations\n",
      "GPT :  213231.839474 s : cg: converged in 133 iterations\n",
      "GPT :  213232.730879 s : cg: converged in 139 iterations\n",
      "GPT :  213235.227432 s : cg: converged in 331 iterations\n",
      "GPT :  213237.343689 s : cg: converged in 332 iterations\n",
      "GPT :  213239.387726 s : cg: converged in 338 iterations\n",
      "GPT :  213241.258662 s : cg: converged in 328 iterations\n",
      "GPT :  213243.151086 s : cg: converged in 319 iterations\n",
      "GPT :  213244.864590 s : cg: converged in 324 iterations\n",
      "GPT :  213246.766639 s : cg: converged in 349 iterations\n",
      "GPT :  213248.462825 s : cg: converged in 342 iterations\n",
      "GPT :  213250.160066 s : cg: converged in 319 iterations\n",
      "GPT :  213251.871504 s : cg: converged in 336 iterations\n",
      "GPT :  213253.860384 s : cg: converged in 333 iterations\n",
      "GPT :  213255.942031 s : cg: converged in 328 iterations\n"
     ]
    }
   ],
   "source": [
    "def effective_mass_pion_single_config(m):\n",
    "    w = g.qcd.fermion.wilson_clover(U, mass=m, csw_r=0.0, csw_t=0.0, nu=1.0, xi_0=1.0,\n",
    "                                isAnisotropic=False,\n",
    "                                boundary_phases=[1,1,1,-1])\n",
    "    invD = w.propagator(inv.preconditioned(pc.eo1_ne(), cg))\n",
    "    dst = g(invD * src)\n",
    "    Ct = g.slice(g.trace(g.adj(dst)*dst),3)\n",
    "    return [np.log(Ct[t] / Ct[t+1]).real for t in range(13)]\n",
    "\n",
    "emp = {}\n",
    "for m in [-0.5, -0.6, -0.7, -0.8, -0.9, -1.0]:\n",
    "    emp[m] = effective_mass_pion_single_config(m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 193,
   "id": "f132b3e9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEGCAYAAACdJRn3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAASYUlEQVR4nO3df5BdZ13H8feHpoUZAVGyKjYt6QwByWAVXSpYmUQpM2mVVhSZRlFRIH9gGRxAreKQbflHZET8UcCAtYBIqcBAhDAFS1McsditQEsSi6GlNLXapSDCMIKxX/+4N3C7bLL3yZ6zdy/7fs2cyX3Oee7Z75OduZ89P+5zUlVIktTiQZMuQJI0fQwPSVIzw0OS1MzwkCQ1MzwkSc02TLqAk7Fx48bavHnzpMuQpKly8803f76qZrrY11SGx+bNm5mfn590GZI0VZLc2dW+PG0lSWpmeEiSmhkekqRmhockqZnhIUlqZnhIkpoZHpKkZoaHJKmZ4SFJamZ4SJKa9RoeSa5Mcm+STy3T70lJjiZ5Vp/1SJK60feRx1XAjhN1SHIK8Crggz3XIknqSK/hUVUfAb6wTLcXAe8C7u2zFklSdyZ6zSPJ6cAzgdeP0XdXkvkk8wsLC/0XJ0k6rklfMH8t8DtVdf9yHatqT1XNVtXszEwn09FLkk7SpJ/nMQtcnQRgI3BBkqNV9Z6JViVJOqGJhkdVnXXsdZKrgPcZHJK09vUaHkneDmwHNiY5AuwGTgWoqjf0+bMlSf3pNTyqamdD3+f2WIokqUOTvmAuSZpChockqZnhIUlqZnhIkpoZHpK0WrZvh82bJ11FJwwPSVpNd9456Qo6YXhIkpoZHpLUt7k5SOCGGwbtZLDMzU2yqhUxPCSpb3NzUAXbtg3aVYPF8JAkrSeGhyStpkc/etIVdGLSU7JL0vqxf/+kK+iMRx6SpGaGhySpmeEhSWpmeEiSmhkekqRmhockqZnhIUlqZnhIkpr1Gh5Jrkxyb5JPHWf7LyW5JcmtST6a5If6rEeS1I2+jzyuAnacYPsdwLaq+kHglcCenuuRJHWg1+lJquojSTafYPtHR5o3Apv6rEeS1I21dM3jecAHjrcxya4k80nmFxYWVrEsSdJiayI8kvwkg/D4neP1qao9VTVbVbMzMzOrV5ykbz/btw8WnbSJz6qb5GzgTcD5VXXfpOuRtE589rOTrmCqTfTII8mZwLuBX66qT0+yFknrzJ13TrqCqdbrkUeStwPbgY1JjgC7gVMBquoNwCuARwKvSwJwtKpm+6xJkrRyfd9ttXOZ7c8Hnt9nDZL0DXNzcNll32wP/miF3bun+nnik7AmLphL0qqYm4Mq2LZt0K4aLAZHM8NDktTM8JC0/uzfPzhVpZNmeEhanzxVtSKGhySpmeEhSWpmeEiSmhkekqRmhockqZnhIUlqZnhIkpoZHpKkZoaHJKmZ4SFJamZ4SJKaGR6SpGaGhySpmeEhabKc3XYqGR6SJmf79gc+FlZTw/CQJDXrNTySXJnk3iSfOs72JPnTJIeT3JLkR/qsR9IaMTcHCdxww6CdDBZPYU2Nvo88rgJ2nGD7+cCW4bILeH3P9UhaC+bmoAq2bRu0qwaL4TE1eg2PqvoI8IUTdLkIeEsN3Ag8Ismj+qxJ0hqyf/+kK9BJmvQ1j9OBu0baR4brvkWSXUnmk8wvLCysSnGSVsHu3ZOuQCdh0uExtqraU1WzVTU7MzMz6XIkdcVTVVNp0uFxN3DGSHvTcJ0kaQ2bdHjsBX5leNfVk4EvVdU9E65JkrSMDX3uPMnbge3AxiRHgN3AqQBV9QZgH3ABcBj4KvBrfdYjSepGr+FRVTuX2V7Ab/RZgySpe5M+bSVJmkKGhySpmeEhSWpmeEiSmhkekqRmhockqZnhIUlqZnhIkpoZHpKkZoaHJKmZ4SFJamZ4SPKZGmpmeEjr3fbtcNllk65CU8bwkCQ1Mzyk9WpuDhK44YZBOxksnsLSGDJ4pMZ0mZ2drfn5+UmXIX37SGAKPwvUJsnNVTXbxb488pAkNTM8JMHu3ZOuQFPG8JDkdQ41MzwkSc16D48kO5LcluRwkkuX2H5mkuuTfDzJLUku6LsmSdLK9BoeSU4BrgDOB7YCO5NsXdTt94FrquqJwMXA6/qsSZK0chuW65Dkmqp6dpJbgdF7+QJUVZ19grefAxyuqtuH+7oauAg4ONKngIcPX38n8O8N9UuSJmDZ8ABePPz3Z05i/6cDd420jwA/tqjPHPDBJC8CvgM4b6kdJdkF7AI488wzT6IUSVJXlj1tVVX3DP+9E/ga8EPA2cDXhutWaidwVVVtAi4A3prkW+qqqj1VNVtVszMzMx38WEnSyRr7mkeS5wP/DPwc8CzgxiS/vszb7gbOGGlvGq4b9TzgGoCq+ifgIcDGceuSJK2+cU5bHfNbwBOr6j6AJI8EPgpceYL33ARsSXIWg9C4GPjFRX0+BzwNuCrJ4xmEx0JDXZKkVdZyt9V9wJdH2l8erjuuqjoKXAJcCxxicFfVgSSXJ7lw2O2lwAuSfBJ4O/DcmsYJtyRpHWk58jgMfCzJexncIXURcEuSlwBU1WuWelNV7QP2LVr3ipHXB4FzG+uWJE1QS3h8Zrgc897hvw/rrhxJ0jQYOzyq6oSPGkvyZ1X1opWXJEla67r8hrmnniRpnXBiRElSM8NDktSsy/BIh/uSJK1hY18wTzILvBx49PB9iydG/JPuy5MkrUUtRx5vA/4K+HngGQwmSnzGsY1VdVWnlUnrkU/005RoCY+FqtpbVXdU1Z3Hlt4qk9ajy054R7y0ZrR8SXB3kjcB1zGYXReAqnp351VJkta0liOPXwN+GNjB4HTVsVNXklZibg6SwQLffO0pLK1hGXcOwiS3VdXjeq5nLLOzszU/Pz/pMqTuJeC8oOpJkpuraraLfbUceXx0ieePS5LWoZZrHk8GPpHkDgbXPMZ5hrmkFrt3T7oCaSwt4bGjtyokDXidQ1OiZVZdb8uVJAHObSVJOgmGhySpmeEhSWpmeEiSmvUeHkl2JLktyeEklx6nz7OTHExyIMnf9F2TJGllWm7VbZbkFOAK4OnAEeCmJHur6uBIny3A7wLnVtUXk3xPnzVJklau7yOPc4DDVXV7VX0duBq4aFGfFwBXVNUXAarq3p5rkiStUN/hcTpw10j7yHDdqMcCj03yj0luTLLklxGT7Eoyn2R+YWGhp3IlSeNYCxfMNwBbgO3ATuCNSR6xuFNV7amq2aqanZmZWd0KJUkP0Hd43A2cMdLeNFw36giwt6r+t6ruAD7NIEwkSWtU3+FxE7AlyVlJTgMuBvYu6vMeBkcdJNnI4DTW7T3XJUlagV7Do6qOApcA1wKHgGuq6kCSy5NcOOx2LXBfkoPA9cBvVdV9fdYlSVqZsR8GtZb4MChJajeph0FJkgQYHpKkk2B4SJKaGR6SpGaGhySpmeEhSWpmeEiSmhke0mJzc5OuQFrzDA9p1PbtcNllk65CWvMMD0lSM8NDgsGpqgRuuGHQTgaLp7CkJRkeEgxCogq2bRu0qwaL4SEtyfCQJDUzPKRR+/fD7t2TrkJa8wwPaTFPVUnLMjwkSc0MD0lSM8NDktTM8JAkNTM8JEnNeg+PJDuS3JbkcJJLT9Dv55NUkk4ezi5J6k+v4ZHkFOAK4HxgK7AzydYl+j0MeDHwsT7rkSR1o+8jj3OAw1V1e1V9HbgauGiJfq8EXgX8T8/1SJI60Hd4nA7cNdI+Mlz3DUl+BDijqt5/oh0l2ZVkPsn8wsJC95VKksY20QvmSR4EvAZ46XJ9q2pPVc1W1ezMzEz/xUmSjqvv8LgbOGOkvWm47piHAU8A9if5LPBkYK8XzSVpbes7PG4CtiQ5K8lpwMXA3mMbq+pLVbWxqjZX1WbgRuDCqprvuS5J0gr0Gh5VdRS4BLgWOARcU1UHklye5MI+f7YkqT8b+v4BVbUP2Ldo3SuO03d73/VIklbOb5hLkpoZHpKkZoaHJKmZ4SFJamZ4SJKaGR6SpGaGhySpmeEhSWpmeEiSmhkekqRmhockqZnhIUlqZnhIkpoZHpKkZoaHJKmZ4SFJamZ4SJKaGR6SpGaGhySpmeEhSWrWe3gk2ZHktiSHk1y6xPaXJDmY5JYk1yV5dN81SZJWptfwSHIKcAVwPrAV2Jlk66JuHwdmq+ps4J3AH/ZZkyRp5fo+8jgHOFxVt1fV14GrgYtGO1TV9VX11WHzRmBTzzVJklao7/A4HbhrpH1kuO54ngd8YKkNSXYlmU8yv7Cw0GGJkqRWa+aCeZLnALPAq5faXlV7qmq2qmZnZmZWtzhJ0gNs6Hn/dwNnjLQ3Ddc9QJLzgJcD26rqaz3XJElaob6PPG4CtiQ5K8lpwMXA3tEOSZ4I/AVwYVXd23M9kqQO9BoeVXUUuAS4FjgEXFNVB5JcnuTCYbdXAw8F/jbJJ5LsPc7uJElrRN+nraiqfcC+ReteMfL6vL5rkCR1a81cMJckTQ/DQ5LUzPCQJDUzPCRJzQwPSVIzw0OS1MzwkCQ1MzwkSc0MD0lSM8NDktTM8JAkNTM8JEnNDA9JUjPDQ5LUzPCQJDUzPCRJzQwPSVIzw0OS1MzwkCQ1MzwkSc16D48kO5LcluRwkkuX2P7gJO8Ybv9Yks191yRJWplewyPJKcAVwPnAVmBnkq2Luj0P+GJVPQb4Y+BVfdYkSVq5vo88zgEOV9XtVfV14GrgokV9LgLePHz9TuBpSdJzXZKkFdjQ8/5PB+4aaR8Bfux4farqaJIvAY8EPj/aKckuYNew+ZUkt/VS8fI2sqi2b3PrbbzgmNeL9Tjmx3W1o77DozNVtQfYM+k6ksxX1eyk61gt62284JjXi/U65q721fdpq7uBM0bam4brluyTZAPwncB9PdclSVqBvsPjJmBLkrOSnAZcDOxd1Gcv8KvD188CPlxV1XNdkqQV6PW01fAaxiXAtcApwJVVdSDJ5cB8Ve0F/hJ4a5LDwBcYBMxaNvFTZ6tsvY0XHPN64ZhXIP6RL0lq5TfMJUnNDA9JUjPDYxlJfiHJgST3JznhbX1JTkny8STvW636ujbOeJOckeT6JAeHfV+82nV2adzf8XJT7UyTJN+d5ENJ/m3473cdp98fDv9vDiX502n+Am/DmM9M8sHhmA9O85RJ44552PfhSY4k+fNx9m14LO9TwM8BHxmj74uBQ/2W07txxnsUeGlVbQWeDPzGEtPOTJNlxzzmVDvT5FLguqraAlw3bD9Akh8HzgXOBp4APAnYtppFdmzZMQ+9BXh1VT2ewSwZ965SfX0Yd8wAr2S8zznA8FhWVR2qqmW/zZ5kE/DTwJv6r6o/44y3qu6pqn8Zvv4yg8A8fTXq68OYv+NxptqZJqPTAr0Z+Nkl+hTwEOA04MHAqcB/rkZxPVl2zMM/CDZU1YcAquorVfXVVauwe+P8nknyo8D3Ah8cd8eGR3deC/w2cP+E61hVw0P6JwIfm3ApfVtqqp2pDUzge6vqnuHr/2DwwfEAVfVPwPXAPcPl2qqa5iPrZccMPBb4ryTvHp6CfvXwqHNaLTvmJA8C/gh4WcuOp2Z6kj4l+Xvg+5bY9PKqeu8Y7/8Z4N6qujnJ9o7L69xKxzuyn4cC7wJ+s6r+u6v6+tDVmKfJicY82qiqSvIt9+wneQzweAYzQwB8KMlTq+ofOi+2IysdM4PPxKcy+IPoc8A7gOcy+D7amtTBmF8I7KuqIy2XtAwPoKrOW+EuzgUuTHIBg8P8hyf566p6zsqr614H4yXJqQyC421V9e6VV9WvDsY8zlQ7a8qJxpzkP5M8qqruSfIolj6v/0zgxqr6yvA9HwCeAqzZ8OhgzEeAT1TV7cP3vIfBdb01Gx4djPkpwFOTvBB4KHBakq9U1QlvCvG0VQeq6neralNVbWbwDfkPr9Xg6MLwjpu/BA5V1WsmXc8qGWeqnWkyOi3QrwJLHX19DtiWZMPwj4VtTPcNIeOM+SbgEUlmhu2fAg6uQm19WXbMVfVLVXXm8PPrZcBblgsOMDyWleSZSY4wSOf3J7l2uP77k+ybbHXdG3O85wK/DPxUkk8MlwsmVPKKjTPmqjoKHJtq5xBwTVUdmFTNHfgD4OlJ/g04b9gmyWySYzd9vBP4DHAr8Engk1X1d5MotiPLjrmq/o/BB+h1SW4FArxxQvV2YZzf80lxehJJUjOPPCRJzQwPSVIzw0OS1MzwkCQ1MzwkSc0MD0lSM8NDktTM8JBWKMnmJP+a5Kokn07ytiTnJfnH4XMUzpl0jVLXDA+pG49hMDPpDwyXXwR+gsG3lX9vgnVJvTA8pG7cUVW3VtX9wAEGD+ApBlN7bJ5oZVIPDA+pG18beX3/SPt+nL1a34YMD0lSM8NDktTMWXUlSc088pAkNTM8JEnNDA9JUjPDQ5LUzPCQJDUzPCRJzQwPSVKz/wcPBSod7VbIjgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "quark_masses = []\n",
    "pion_masses = []\n",
    "pion_masses_err = []\n",
    "\n",
    "for m in emp:\n",
    "    quark_masses.append(m)\n",
    "    pion_masses.append(0.5*(emp[m][7] + emp[m][8]))\n",
    "    pion_masses_err.append(0.5*(emp[m][7] - emp[m][8]))\n",
    "    \n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "plt.ylim(0,1.5)\n",
    "plt.xlim(-1.5,-0.4)\n",
    "plt.xlabel(\"m\")\n",
    "plt.ylabel(\"m_pi\")\n",
    "ax.errorbar(quark_masses, pion_masses, pion_masses_err, marker='+', ls='', c='red', label=\"m_pi(m)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a618f9ac",
   "metadata": {},
   "source": [
    "We note that for quenched QCD we may have approximately massless Wilson fermions at around $m \\approx -1.2$.\n",
    "The mass $m_c$ at which the quark effectively becomes massless is also called the **critical mass**.  The fact that we can use the pion as a proxy is due to the pion being a Goldstone boson of chiral symmetry in this limit.  We will return to the topic of symmetries and in particular also chiral symmetry in later chapters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "080f9e27",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.9.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
