{
 "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 8: strong coupling expansion\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d5caf31",
   "metadata": {},
   "source": [
    "In this chapter we perform a strong-coupling expansion of the path integral and the expectation value\n",
    "of a Wilson loop\n",
    "$$\n",
    "\\begin{align}\n",
    "W_C = \\textrm{Tr}\\left[\\prod_{\\mu_i,x_i,m_i \\in C} U^{m_i}_{\\mu_i}(x_i)\\right] \\,,\n",
    "\\end{align}\n",
    "$$\n",
    "where $m_i \\in \\{ \\dagger, 1 \\}$.\n",
    "The path $C$ could, e.g., describe a plaquette, a rectangle, or more complicated closed paths.\n",
    "\n",
    "We will perform the expansion using the Wilson gauge action\n",
    "$$\n",
    "\\begin{align}\n",
    "S_W &= \\beta \\sum_x \\sum_{\\mu<\\nu} ( 1 - P_{\\mu\\nu}(x)) \n",
    "\\end{align}\n",
    "$$\n",
    "with scalar plaquette\n",
    "$$\n",
    "\\begin{align}\n",
    "P_{\\mu\\nu}(x) &= \\frac1{N_c} \\textrm{Re}\\: \\textrm{Tr}\\: U_{\\mu\\nu}(x) \\,, \\\\\n",
    "U_{\\mu\\nu}(x) &= U_\\mu(x) U_\\nu(x+a \\hat{\\mu}) U_\\mu^\\dagger(x + a \\hat{\\nu}) U_\\nu^\\dagger(x) \\,,\n",
    "\\end{align}\n",
    "$$\n",
    "and\n",
    "$$\n",
    "\\begin{align}\n",
    "\\beta &= \\frac{2 N_c}{g^2} \\,.\n",
    "\\end{align}\n",
    "$$\n",
    "The expansion in a large coupling parameter $g$, therefore corresponds to an expansion in a small\n",
    "number $\\beta$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "81cad55a",
   "metadata": {},
   "source": [
    "The expectation value of a Wilson loop is then given by\n",
    "$$\n",
    "\\begin{align}\n",
    "\\langle W_C \\rangle &= \\frac{\n",
    "\\int d[U]  \\textrm{Tr}\\left[\\prod_{\\mu_i,x_i \\in C} U_{\\mu_i}(x_i)\\right] \n",
    "\\exp\\left[ -\\beta \\sum_x \\sum_{\\mu<\\nu} ( 1 - P_{\\mu\\nu}(x)) \\right]\n",
    "}{\n",
    "\\int d[U]\n",
    "\\exp\\left[ -\\beta \\sum_x \\sum_{\\mu<\\nu} ( 1 - P_{\\mu\\nu}(x)) \\right]\n",
    "} \\\\\n",
    "&= \\frac{\n",
    "\\int d[U]  \\textrm{Tr}\\left[\\prod_{\\mu_i,x_i \\in C} U_{\\mu_i}(x_i)\\right] \n",
    "\\exp\\left[ \\beta \\sum_x \\sum_{\\mu<\\nu} P_{\\mu\\nu}(x) \\right]\n",
    "}{\n",
    "\\int d[U]\n",
    "\\exp\\left[ \\beta \\sum_x \\sum_{\\mu<\\nu} P_{\\mu\\nu}(x) \\right]\n",
    "} \\\\\n",
    "&= \n",
    "\\int d[U]  \\textrm{Tr}\\left[\\prod_{\\mu_i,x_i \\in C} U_{\\mu_i}(x_i)\\right] \n",
    "\\exp\\left[ \\beta \\sum_x \\sum_{\\mu<\\nu} P_{\\mu\\nu}(x) \\right] ( 1 + O(\\beta)) \\\\\n",
    "&= \n",
    "\\int d[U]  \\textrm{Tr}\\left[\\prod_{\\mu_i,x_i \\in C} U_{\\mu_i}(x_i)\\right] \n",
    "\\exp\\left[ \\frac{\\beta}{N_c} \\sum_x \\sum_{\\mu<\\nu} \\textrm{Re}\\: \\textrm{Tr}\\: U_{\\mu\\nu}(x) \\right] ( 1 + O(\\beta)) \\\\\n",
    "&= \n",
    "\\int d[U]  \\textrm{Tr}\\left[\\prod_{\\mu_i,x_i \\in C} U_{\\mu_i}(x_i)\\right] \n",
    "\\exp\\left[ \\frac{\\beta}{2N_c} \\sum_x \\sum_{\\mu<\\nu} (\\textrm{Tr}\\: U_{\\mu\\nu}(x) + \\textrm{Tr}\\: U^\\dagger_{\\mu\\nu}(x)) \\right] ( 1 + O(\\beta)) \\\\\n",
    "&= \n",
    "\\int d[U]  \\textrm{Tr}\\left[\\prod_{\\mu_i,x_i \\in C} U_{\\mu_i}(x_i)\\right] \n",
    "\\exp\\left[ \\frac{\\beta}{2N_c} \\sum_x \\sum_{\\mu \\neq \\nu} \\textrm{Tr}\\: U_{\\mu\\nu}(x) \\right] ( 1 + O(\\beta)) \\\\\n",
    "&= \n",
    "\\int d[U]  \\textrm{Tr}\\left[\\prod_{\\mu_i,x_i \\in C} U_{\\mu_i}(x_i)\\right] \n",
    "\\sum_{n=0}^\\infty \\frac1{n!} \\left(\\frac{\\beta}{2N_c}\\right)^n (\\sum_P \\textrm{Tr}\\: U_P)^n ( 1 + O(\\beta)) \\,,\n",
    "\\end{align}\n",
    "$$\n",
    "where in the last line we sum over all possible plaquettes $P$ (every spacetime-point and every $x$)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "645aa455",
   "metadata": {},
   "source": [
    "Let us now first consider the leading order $n=0$ term for which we need to perform the integral\n",
    "$$\n",
    "\\begin{align}\n",
    "\\int d[U]  \\textrm{Tr}\\left[\\prod_{\\mu_i,x_i \\in C} U_{\\mu_i}(x_i)\\right] \\,.\n",
    "\\end{align}\n",
    "$$\n",
    "Note that for any closed Wilson loop, each link $U_\\mu(x)$ or its adjoint only appears once.  In the case of the plaquette, e.g.,\n",
    "$$\n",
    "\\begin{align}\n",
    "U_{\\mu\\nu}(x) &= U_\\mu(x) U_\\nu(x+a \\hat{\\mu}) U_\\mu^\\dagger(x + a \\hat{\\nu}) U_\\nu^\\dagger(x)\n",
    "\\end{align}\n",
    "$$\n",
    "such that if we integrate over $U_\\mu(x)$ first, we are performing an integral of the type\n",
    "$$\n",
    "\\begin{align}\n",
    "M_{ab} = \\int d[U_\\mu(x)] (U_\\mu(x))_{ab} \\,.\n",
    "\\end{align}\n",
    "$$\n",
    "One can show that for any $V \\in $ SU(3), we have\n",
    "$$\n",
    "\\begin{align}\n",
    "(VM)_{ab} = \\int d[U_\\mu(x)] (V U_\\mu(x))_{ab} = \\int d[U_\\mu(x)] (U_\\mu(x))_{ab} = M_{ab}\n",
    "\\end{align}\n",
    "$$\n",
    "due to the invariance of the Haar measure.  In matrix form, we have\n",
    "$$\n",
    "\\begin{align}\n",
    "VM = M\n",
    "\\end{align}\n",
    "$$\n",
    "for any group element $V$.  One can show that this can only be satisfied for\n",
    "$$\n",
    "\\begin{align}\n",
    "M = 0 \\,.\n",
    "\\end{align}\n",
    "$$\n",
    "The leading term $n=0$ therefore vanishes.\n",
    "\n",
    "Before moving on to the $n=1$ term, we first verify the integral equation numerically by a straightforward extension of our Haar integration routine to SU(3)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "id": "c5022108",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\\int dU U_00 (0.012514483442892249, 0.0040934094282299)\n",
      "\\int dU U_01 (0.006719000502069609, 0.004046532911347319)\n",
      "\\int dU U_02 (-0.001452925711356024, 0.004136127321453077)\n",
      "\\int dU U_10 (0.0009426477768964693, 0.0040616879424338545)\n",
      "\\int dU U_11 (-0.00044307235985169323, 0.004088854622644644)\n",
      "\\int dU U_12 (-0.002986773309435577, 0.004102720114800084)\n",
      "\\int dU U_20 (-0.005149322584176619, 0.004073162320270486)\n",
      "\\int dU U_21 (-0.00013018166498728664, 0.004102016198185097)\n",
      "\\int dU U_22 (-0.0015228645553870606, 0.004093229510335719)\n"
     ]
    }
   ],
   "source": [
    "import gpt as g\n",
    "import numpy as np\n",
    "import scipy.linalg\n",
    "\n",
    "su3_generators = g.object_type.ot_matrix_su_n_fundamental_algebra(3).generators(np.complex128)\n",
    "\n",
    "U = [np.matrix(scipy.linalg.expm(2j*sum([np.random.normal()*p.array for p in su3_generators]))) for i in range(40)]\n",
    "for u in U:\n",
    "    assert np.linalg.norm(u * u.H - np.identity(3)) < 1e-14\n",
    "\n",
    "np.random.seed(13)\n",
    "\n",
    "def haar_integrate(f, U, N, NskipToMakeIndependent):\n",
    "    u0 = U[np.random.randint(0,len(U))]\n",
    "    val = 0.0\n",
    "    val2 = 0.0\n",
    "    n = 0.0\n",
    "    for i in range(N):\n",
    "        v = f(u0).real\n",
    "        val += v\n",
    "        val2 += v**2\n",
    "        n += 1\n",
    "        for j in range(NskipToMakeIndependent):\n",
    "            u0 = U[np.random.randint(0,len(U))] * u0\n",
    "    # just use biased estimator for error\n",
    "    return (val/n, (val2/n - (val/n)**2)**0.5 / n**0.5 )\n",
    "\n",
    "for a in range(3):\n",
    "    for b in range(3):\n",
    "        print(f\"\\int dU U_{a}{b}\",haar_integrate(lambda u: u[a,b], U, 10000, 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28c56bef",
   "metadata": {},
   "source": [
    "We now consider the next term in the strong coupling expansion, i.e., the $n=1$ term.  We have to compute\n",
    "$$\n",
    "\\begin{align}\n",
    "\\int d[U]  \\textrm{Tr}\\left[\\prod_{\\mu_i,x_i \\in C} U_{\\mu_i}(x_i)\\right] \n",
    " \\frac{\\beta}{2N_c} \\sum_P \\textrm{Tr}\\: U_P ( 1 + O(\\beta)) \\,.\n",
    "\\end{align}\n",
    "$$\n",
    "Let us for concreteness look at the term, where $P$ is the counter-clockwise plaquette at the origin in the $0-1$ plane\n",
    "and that we consider a wilson loop in the $0-1$ plane with size $6 \\times 4$.  We plot this scenario below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "id": "a9433fed",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAU8AAADnCAYAAACe9dTFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAANuklEQVR4nO3de8wldX3H8feXvbKLGGVXXFgvNdTVotRLisXCYtHWKllo7IVNjFZqm4akRZMqBVPaSi0pLLghJm60F2ijDaFNWmkUKTa0UHuJaU1EsBayjaIULYLA7sJef/1j5oHjdnfPzPfMeWbOc96v5OScM8/M/L6Tc+bzfM+cy0QpBUlSO8f1XYAkzSLDU5ISDE9JSjA8JSnB8JSkBMNTkhIMT0lKGGR4RsTZEfHjfdfRhYhYFhE/HxEv7ruWLkTECRHxrog4oe9auhARp0bE1ogY5L7QVkS8NiLO67uOLkRlS0Rs6ruWIxnUE6YOzX8B7gQ+2Xc9k6hDcyuwE7gZeF/PJU2kDs0rgIeAG4GZ3kHr0Pwk8ADwF8Dzey5pInVo3gH8G/CpvuuZRB2aFwBfB/4a+P2eSzqi5W0XiIgoHX8tKSLOBrYBZwBrulz3mHEDoMvtiYhlwC8A11DtkIvWoU3psTkB+A3gCmAZ1ePzeJdjHGPsaWzPqcDvAu+iah5WAk93OcYxxp7G9rwWuBZ4I3A8EF2uf8zYnW5PvT9uAa4DNrCI+05GtNn2iHgTVVeoQQrAr9tKOfHNUg69pOncbTvPM+rrdS2XG+dngKuAk4G19bQHgLd1PM7hlgHfBw50uM5VVF3N5cAKnt2em4A/6HCcIyjLgEc6XulJwPuB91B1aquAJ6i2746OxzpcoXp8uvTDwJVUhx1WUD0H9gKbgUc7HutwB4AnO17nZqp952U8+6rtEeCsjsc5XAC7gH0drnM51au2K6n2m4V953bg1zsc50j+C0q79yVKKY0vwKXUr3K7vlA9GBdUG8Fe4CvTGGexLsBqqpe7j1LtNNf3XdOE27Me+CiwG9gPXNB3TRNuzyupjqftAQ4C6/quaYJtCeDNwJepDkE81HdNE27PCuBi4OH6uXbLIox5sG22tX3ZfilwQyllasdVRo57HCilfG5a4yyWiFgN/BJwVynla33XM6mIWE/VWd9USpl2pzZ1EfFKquDZUUo52Hc9k6j3nfOAk0opt/Rdz6QiYgXwTuDeUsqXpjzWQeC4Ntk2uPCUpMWWCc9BfVRJkmaF4SlJCYanJCUYnpKUYHhKUoLhKUkJhqckJRiekpRgeEpSguEpSQmGpyQlGJ6SlGB4SlKC4SlJCYanJCUYnpKUYHhKUiILDU9Jc60+ffPC7dc1Xc7wlDTvrhm5fW3ThQxPSXOr7jrPHpl0VtPus1F4RsSKiFgHnFjfX1dfPBGcpFl2DdVpwhccT8Pus2nn+WHgO/U1wIPAd4E3NVxekgalbv5eCOwZmbwH2NCkMWx06uGIOAf4PLBmZPIuYGMp5fFWFUvSwEREAej81MOllLuBe0YmPQVcb3BKmleNOk94pvv8O6rjA3adkpaMqXWe9UrvBnbWd+06Jc21th9V+nh9vb3rQiRplrQNz68C2HVKmnd+SF6SEgxPSUowPCUpwfCUpATDU5ISDE9JSjA8JSnB8JSkBMNTkhIMT0lKMDwlKcHwlKQEw1OSEgxPSUowPCUpwfCUpATDU5ISDE9JSjA8JSlhecv5T59KFZLUr6dp2Uy2Dc+280vSLFhJy/D0ZbskJRiekpRgeEpSguEpSQmGpyQlGJ6SlGB4SlKC4SlJCYanJCUYnpKUYHhKUoLhKUkJhqckJRiekpRgeEpSguEpSQmGpyQlGJ6SlGB4SlKC4SlprkVEUGdhfbsRw1PSvDt/5PaWpgsZnpLmVt1pXjcy6bqm3afhKWmenQ+cOnJ/Aw27zyiljJ8p4oXAmcDPAe8GLqz/dFspZX+rUiVpAOoO82vApsP+dD+wqYwJx+UNx/lN4H3AQlD+OfAc4O3A7Y2rlaRhuQv4DrB55P79TRZs2nmeBtwDrB6Z/DDwolLKgValStLAREQBKKV0+257KeUB4DPAwXrSLuByg1PSvGrUecIz3ee9wErsOiUtIVPrPOuVPgD8U33XrlPSXGv7UaU/qq8/3XUhkjRL2obn/wDYdUqad35IXpISDE9JSjA8JSnB8JSkBMNTkhIMT0lKMDwlKcHwlKQEw1OSEgxPSUowPCUpwfCUpATDU5ISDE9JSjA8JSnB8JSkBMNTkhIMT0lKMDwlCZqdRnhE2/A8ExqfmVOSZkXrYGsbnnsTAS1JQ3eo7QK+bJekBMNTkhIMT0lKMDwlKcHwlKQEw1OSEgxPSUowPCUpwfCUpATDU5ISDE9JSjA8JSnB8JSkBMNTkhIMT0lKMDwlKcHwlKQEw1OSEgxPSUowPCXNtYhYRZ2FEbG66XKGp6R5996R27/adCHDU9LcqrvOq0Ymfbhp92l4Sppn7wVGw3IlDbvPKGX8edgj4seArcBPAG8APlr/6epSyvdalSpJAxARATwMrOPZRvIQ8BiwvowJx6bheRXwIWDZyOQCnFtKuTtRtyT1LiIuAl4K/GE96XLgG6WUm8cu2zA8nwc8CKwdmfxl4PXj0lmShi4iCkApJZou0+iYZynlMeAGYF89aTfwAYNT0rxq1HnCM93nQ1QHV+06JS0ZU+s865U+BvxVfdeuU9Jca/tRpYWDqHd2XYgkzZK24bkLwK5T0rzzQ/KSlGB4SlKC4SlJCYanJCUYnpKUYHhKUoLhKUkJhqckJRiekpRgeEpSguEpSQmGpyQlGJ6SlGB4SlKC4SlJCYanJCUYnpKUYHhKUsLylvM/ZypVSFK/9gONz5wJ7cNzU8v5JWkWLKdleLZ92b6/5fySNAtan9TSY56SlGB4SlKC4SlJCYanJCUYnpKUMLjwjIgNEbEjIj7Sdy2SdDSDCc+F0AR2Ar8CXNBzSZJ0VG0/JN9YRHwQeFnDGs4CXlHfX1Zfvzoivpoc/gng/FLKY8nlJemYphaewLUdrOP0CZbdAWztoAZJ+n+m/bL9tlJKjLsAJwJXAruAp+pl72my7FHWB2DXKWlqBnHMs5TyZCnlI8ApwNXAbmBPv1VJ0tENIjwXjIToBuCtfdcjSUczzWOeaaWUJ/uuQZKOZVCdpyTNCsNTkhIMT0lKMDwlKcHwlKSEpRyeG/suQNLMaHX+IliC4RkRm+ubmyNiyW2fpG5FxDrq8IyI9U2XW4rhsq2+Ph742R7rkDQbrhi5/aGmCy2p8Ky7zlfVd1cA2+w+JR1N3XVeMjLp15p2n0stWLYBa0buvwC7T0lHdzk/eLzzOBp2n1HK+NMVR8S7gQ9QvQnzPOAeqvMcX1RK+c+R+U4GPgusBF5dT74H+HQp5ZomBWVFxEbgQarf8jwROEh1nvkvllLeMs2xJc2eiAjgEWAVsLaevBvYB5xUxoRj0++2rwJeXl9DFYwHjzDfofpvK0embQJOaDhOWinlWxFxJrAauAv4LnAR8N/THlvS7CmllIg4G1hHlRkAbwMeGRec0LzzXEnV1b2gnnQQuLWU8o4jzLud6hjCQtDuBjaWUr4/dqCOREQB7i2lvGrszJLmXp0ZjPwe8FiNjnmWUvYBv82zP1S8j6MfF7iaqgOF6mXz9sUMTklaDG3eMPozYG99+/OjxzpHlVL+F/jEwl3g+nx5kjRMjcOz7j7/uL477t2oq+vru+06JS1FbT+qdBvA0brOBXX3CXBnpihJGrq24Xmkd9iP5UDL+SVpJiy1D8lL0qIwPCUpwfCUpATDU5ISDE9JShhkeEbEhoj4y4h4c/3lfUkalEGGJ3AasAX4DPAfhqikoWn6q0pZF0fEuYnlNlAF+1rgNcDngO9FxN8A32i4jqH+Y5C0BEwzPPdS/Rzdpg7WtZIqUC8ZN+OIP+1gXEk6oqmFZylldXbZiDgH+FvgucAu4NtUP8b82Sa/sydJ0zbtl+1ZB6iC8+sYmpIGaKjh+a9UJ3K7z9CUNESDDM86MO/tuw5JOhrfkZakBMNTkhIGGZ4RsSwi3hoRq8bPPXwRcWJEnLdUPugfET8UET/adx1diYg3RMQpfdfRhaj8VESsHT/38EXE8RHx0xExuKwaXEG1N1J9MP7bEXHJEgjRS4AvAPdFxPlLIERvpPrm1+0R8Zq+i+nAF4CdEbEjIjb0XcyEXg7cDjwUER9cAiG6leoMFjsj4heHFKKDKeQIngROAraxNEL0IPAK4GZmP0SD6rnzFuCfl0CIBtWpsn+ZpRGiu4ATgd+jCtHLZjxEnwJeAvwJAwrRRudtf2bm6quW/zC1aiSpP0+UUp7bdOa24RnAhcCLE4W1sQl4D7CG6hzwB4AvUZ3S+JsdjhPA65j+uZa2AD8JLKM6530BbgU+BTze4ThrqD4fu3fcjBO6DDi1vv00sB+4iWqb9nc4zinAyVTPgWm6YeT2U8CjwMeBL3Y8zmlUv9cwTacA76fqpAvV8+1eYAdwf8djnUG1D03TucDbgdVU++kh4A6qQ0ePdDjOcuArIyevHK+UMrgLcA7Vg76H6peVfqTvmibcnt+qH/g9wHZgfd81Tbg9/1g/Po8ClwKr+65pwu3ZBeymCpcLqZuKWbxQNR4L+87fA6/vu6YJt+diqn/IT1Gd+nxj3zUtXFp1noslIp4PXAHcWEq5r+96JhURpwPvBLaXNv/ZBioi3gG8CPhEKeXpvuuZVERcRvVV4FvLEHeIFupjm78D3FJK+fe+65lURLyU6g3Xj5VSvtVzOT9gkOEpSUPX+ztWkjSLDE9JSjA8JSnB8JSkBMNTkhIMT0lKMDwlKcHwlKSE/wMWsTcjr/yE9gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def draw_U(x,mu):\n",
    "    plt.arrow(x[0],x[1],mu[0],mu[1])\n",
    "    plt.arrow(x[0],x[1],mu[0]*0.5,mu[1]*0.5, head_width=0.1, head_length=0.1, fc='k', ec='k')\n",
    "    \n",
    "def draw_P(x):\n",
    "    a = 0.2; oma = 1.0 - a; l = 1.0 - 2*a;\n",
    "    draw_U([x[0]+a,x[1]+a],[l,0.0])\n",
    "    draw_U([x[0]+oma,x[1]+oma],[-l,0.0])\n",
    "    draw_U([x[0]+oma,x[1]+a],[0.0,l])\n",
    "    draw_U([x[0]+a,x[1]+oma],[0.0,-l])\n",
    "    \n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "plt.axis(\"off\")\n",
    "ax.set_aspect('equal')\n",
    "\n",
    "nx = 6\n",
    "ny = 4\n",
    "plt.xlim(-0.1,nx + 0.1)\n",
    "plt.ylim(-0.1,ny + 0.1)\n",
    "\n",
    "for x in range(nx):\n",
    "    draw_U([x,0],[1,0])\n",
    "    draw_U([x+1,ny],[-1,0])\n",
    "for y in range(ny):\n",
    "    draw_U([0,y+1],[0,-1])\n",
    "    draw_U([nx,y],[0,1])\n",
    "    \n",
    "draw_P([0,0])\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83c452f0",
   "metadata": {},
   "source": [
    "In this case there are a lot of links that only appear once and therefore their integrals vanish.  In addition, however, there are some links that appear twice as in the integral\n",
    "$$\n",
    "\\begin{align}\n",
    "M^{(2)}_{ab,cd} = \\int d[U_\\mu(x)] (U_\\mu(x))_{ab} (U_\\mu(x))_{cd} \\,.\n",
    "\\end{align}\n",
    "$$\n",
    "One can show that such integrals have to vanish.  We check this again numerically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "id": "0a1a251e",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(13)\n",
    "for a in range(3):\n",
    "    for b in range(3):\n",
    "        for c in range(3):\n",
    "            for d in range(3):\n",
    "                v, e = haar_integrate(lambda u: u[a,b]*u[c,d], U, 4000, 2)\n",
    "                if abs(v) >= 3*e:\n",
    "                    print(a,b,c,d, v, e)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "484b959c",
   "metadata": {},
   "source": [
    "We see that the configurations with minimal power of $\\beta$ that can survive are those where each link in $W_C$ is matched by the inverse link of a given plaquette.  The corresponding configuration is plotted below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "id": "b070bc5f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAU8AAADnCAYAAACe9dTFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAusUlEQVR4nO2debgdVZW333UzD5AAAWVwYFBQBqdPBAQCCIoCgcZPpRsVQRywFRyYosgQBIUEbaRbxAlQRAaxhRZaBGVUnJoPERGUSYbIkIaEDBeS3Lu+P9Yuz7mXm3tr2OfsqvOs93nOU1VJ1d5r1d61alXd2r8tqorjOI5TjL7UBjiO4zQRD56O4zgl8ODpOI5TAg+ejuM4JfDg6TiOUwIPno7jOCXw4Ok4jlOCWgZPEdlZRHZIbUcMRGSciPxfEXlpaltiICLTReR9IjI9tS0xEJGNReQgEanltVAUEXmdiOyR2o4YiLGfiGyZ2paRqFWHCUHzNuAG4Bup7alCCJoHAQ8AlwBHJTapEiFozgUWAucDjb5AQ9D8BnAfcDGwbmKTKhGC5nXAb4CLUttThRA05wD3Av8JnJrYpBEZX/QAERGNPCxJRHYG5gPbAVNjlj1GvQIQ0x8RGQe8CzgDuyC7lqF1qG2mA58A5gLjsPZZErOOUeruhD8bAycB78OSh4nAczHrGKXuTvjzOuBMYCdgCiAxyx+j7qj+hOtxP2ABsCFdvHbKIEV8F5HdsKzQqSUC+HBbxymHPKw6+LK8exfNPLcLy1kFjxuLvYF5wIuAaeHf7gPeHrme4YwDFgOrI5Y5Cctqjgcm0PLnAuC0iPWMgI4DFkUudD3gk8AHsExtEvAs5t91kesajmLtE5NXAJ/HXjtMwPrA88CuwNOR6xrOamBp5DJ3xa6dzWg9tS0Cdoxcz3AEWAasjFjmeOyp7fPYdZNdO9cCH49Yz0j8BbTY3yVUNfcPOJLwlBv7hzXGHHOC54E7O1FPt37AZOxx92nsojkrtU0V/Vkf+DKwHFgFzEltU0V/XoW9T1sBDACzUttUwRcB3gL8P+wVxMLUNlX0ZwJwKPB46GuXdaHOgaKxrehj+5HA2arasfcqbe89VqvqNZ2qp1uIyGTgEOBmVf1zanuqIiLrY5n1Bara6Uyt44jIq7DAc66qDqS2pwrh2tkDWE9VL0ttT1VEZAJwMPAnVf1dh+saAPqKxLbaBU/HcZxuUyZ41upTJcdxnKbgwdNxHKcEtQ2eItK17z27QS/5IyLjRWRSajtiISJTs29+e4Fe6mtQX39qGTxF5DXAEhG5oq5Ds4ogIh8CnhaR+SIS+zOvFPwAWCgiH+uRIHoPcLeI7NP0IBqGAT8jIteKyGtT21MVETkAWCwiXxeRjVLb004tgyewNtAP7A/c0QNBdF3sG7aPA3/rgSC6AebTmcBjPRBE1wW2wobRNj2ITsE+79kT+FUPBNF1sE/9DgXur1UQLfgtVO7vPJG+B7GPnFP9on2HSOtD6pT+TIvoz56IDCb05fHI3+h9GelL2DZyfVR/pO/2tH1NPh63fWRpWn/YNEcfqs93niKiWMApM3pnEhawsnqyC/1B8o0C2R74s6q+ukTdL0BEJmK+LC9ZxLRh24PYyIz7sI+0x2J7YFdVvaVk/UMQkWOwrDGmP8swf8Zq75cAG+bpQ3kRkUXYSKgY/mQXxCKsv43FG4Bxkf1RLHssM3pnAjZeP2MwLB8Gnsxx/PbAElWdWaLuEQn+9LfZUoRpWJu0x4LVwP3kG621PXCaqp4who2FP1XqXOZpDh9XMpPYBQsqy7Bg+QlgcoHjFbgr3p2TiaHMSSWPPw67sy3FRlDNIXxjW8CfXSL685m87biG42/COvEK4Hrg9QWOnV+l7jWU+URFf7Jhhv3AN4FNCh4b2x8Fvlfy2C2DH8uxYbRzgekF617SAX+2LXnsoVjishS7AfwzdrMqUve8HPsVzjwLqyp1icewu/7XgW+qaldUbzrIXcCfgBOA/9LQWg3mV1hmdKyq3p7amAjcBjwEnKKqjya2pSpP05LZO0dVlyW2pyr3An/FdCEu0xqNAqtl8FTVB4CtU9sRC1W9Grg6tR2xUNW5qW2IiaruldqGWKjqU8C2qe2Ihar+CtgmtR0jUde/tjuO49QaD56O4zgl8ODpOI5TAg+ejuM4JfDg6TiOUwIPno7jOCXoSPAUkWzERqfnUekWrxm2bDo7wT9mxuwFZqQ2IBYikn0+2BN9TUReFlbrHgsKx8JOZZ6fCMu9Gi4YkXF6WH4xqRURCO2xd9js9KRaHUdEtsKG82ZTWDedg8PyVSIyM6UhkZgXlkfXVWwlTN+crb8+73HRg2fIOrOPqMcDh8euo5uEE/vmsLlj+4luKIfTGhwxt+0poamc3rY+P5kVEQhZ55eyTeCYhOZUJmSd7w6bGwN1HYxwRtv6mXkP6kTm+QlaF+dEYF7Ds88zCJkNNiNm7pNbN0I7zKMlHDGe1lNC4whZ595t/7Rdw7PPg4HsVco44KiGZ5/zMD/ApkVeULfsMyRD7X1mx7zZZ67gKSITgv7k2mF7VvjJsP2mAZ+lNX80mL5gV7JPEVmnTSdzYrBxSoXyXofNi52dJwF26Ub2KcasNn/WCdtVhtQejrVHxlTgs93IPkVkYvBlWtieJSLrVSz2dIYqCE2lS9mniEwL/vSF7VkiUvrda1vW2f4eehxdyj5FZEZbXxsX/Cmt4B6yzoMwlaeMLehS9jns2pketieOsOsZWFKUMYWcCVLezPMUTLnmlLD9CCZvtduw/TbA1EnapcEGMaHZjiIiGwNPBdsAXgH8HbiiQrGbYfJXmWzcirD9igpl5mV/7Bxn/vwQa4NPVihzS4bKgi0P2xtUKDMvC7B5uI8I248AT4nImyqUuQk2T3nGcqDqDSYvvwQW0roZPQo8XiHgzMCCZfu1M0AX+lr4w+HjmA9gN7iFwM0Vin0ZplSV+TOA9bVXVigzFyKyE0OvnU9i/n1p2H4CvJihspArgA1zZcg5ZZ12wU5Cu8DoUmDGGFJQ18SUthrDxj5MTabdxmXARyKUnWkKrtVFfzYIDdnuzwpguwhlX0VkGbUcdb41tEe7P4uBqRHKLiwnFqHO04a1zyBwS6SyFXigy/7cNkJfOzmiPwd20ZdpwJIRYsEeY9hYqA/lyjzVRHj/2PZP/cBZqrokz/HdQFUHsUecdgmufuD8NBZVQ1WfxLQlnw//NAjcpKp3prOqEtcBD7RtrwBOVdU8YtB1ZD5Ds/h+mv0HnmMYmoENAl9JZEslVHU59mVMuz9/AW6IWU+RPxgdQ+sRaYB6ntgrsUdbML3Jz6tqGTXuunA6rQt0JSaq3EjUbu9H0+pDq4Bz01lUDVVdDJxNq31uV9Vfp7OoGqp6K5DdmAepWXJUgnOwOAXW544OfTAauYNnyD6zzKGWJ7Yt+wQ7cd9JaE5lVPUJ4Fth8zcNzjozrqM1FcQXGpx1ZixoW29y1pnR7sOXk1kRgbbsE+z9bdSsE4p/qvS1sKxj1plxZVje1vCsM+O0sDxj1L0aQLjzZ69Rvjbavk1AVZ/BLkyanHVmhOwT4OE6JkclOCcsvxU764TiwfMugDqf2JB9gsn3N56QfYLNLd4L/BGgB7LOjIWpDegATZ+KBABtTUFyVyfKd2EQx3GcEnjwdBzHKYEHT8dxnBJ48HQcxymBB0/HcZwS1DJ4BlGMcWPv2Qy6NNa6a7g/9aUXr526KTFl1Cp4hoY/ABtKdVVicyojIpNF5CjgSRH5VGp7qiIiG4jI2cBiEam7MviYiMjWInIVJlCydmp7qhCunb2AO7Bx6o0mKLkdDjxGTXVaa3HHDXeW/bGT9GJMlqs/qVEVEJHJwEeAkzFJrmmYAlAjEZENgM8BH8JuuKuB9ZMaVQER2RobfbInJmk3wFBpu8YQrp09sdFOm2N97e9JjaqAiEwADsEGh0wLv5cmNWoNdDp4vl1Eyn7Zv22FYyFoj0bm2QpPEJ8WkU9XqDvmo4sCVDi/k4Ara/Q0tQqYVMGfcVj2GdGkymxawZ8NK147a1U4dk1cUeH8vquiP5PH3qU4nQyeV2ITjY3ltAAzGSqamjE4wr/lQJaBRhubq6rLReRuLNsay5/xmD/CCwNeSX94jLgjWT6GTciVp0NOwm5EI/X8su1zU7nj1sjh2FC8PMNxpzNUrDsjky4ryiDSt2Ds3QpxIbAP+c7vTEbOmsu2zXOgc8ferxC/wbRxxzq/4zB/+oh27chToNHHtQP59DzbNO9m0yHdROxivpWWbuid3dL/64AvfcC7gIcw3VPFxFSS21bSn2nAsZhG4nJM6mtOarsq+LMRpujUj0n+rQRmpbargj+vAX4a2mUQWJjapgq+CHbj+HPbtXNZF+rtjJ5nN1DV21R1Z+z9zS8Zqh/aKFR1UFUvx+62h2FqVH9Oa1V5VHW5qp4JbIi9x32aBo9/VtWFqnoE1j7fxnx5bvSj6ouq/kFV98YSkJ/RobHc3UCNq4FXA+8B7qamsUBC1M23s8hs4EZVrdXLIcdxnCpk71SLxLbaZJ6O4zhNwoOn4zhOCTx4Oo7jlMCDp+M4Tgk8eDqO45TAg6fjOE4Jio4w2rojVjiO46TlOQomk0WDZy2ERBzHcSIzkYLB0x/bHcdxSuDB03EcpwQePB3HcUrgwdNxHKcEHjwdx3FKUMvgKSJ9IrJzr0zMJSJTReRNqe2IhYhsLCJbprYjFiKynYjMSm1HDMJcRjuJyKTUtsRARCaKyJvrOAlcLYMnpkB/M/CwiLy/B4Lox4HbROT3IrJbYlti8D3gbhG5okeC6K3AIyIyvweC6Cswfx4TkY/1QBD9F+AWrL/tU6cgWtfgKcCzmPjuf9D8ICrYJGNvAH7SA0F0HNZ39gfu6IEg2ofNc/Nx4G8ND6ICLAPWA86k+UFUMIX8rYBLqFMQLShVfyQ5peqx6QA04e/UiBL94xP7osCGEf15f2p/Ik+h8HNMzDbV74nI/ixP5Ef2uyiyP0n7GrBbDhsHivbLTmZygj16/7bEsZti85hMxpxfBdwLXAD8LcfxPwTeCXy+RN0jkWXoX8FOclF2B16LZWyrMZ9uwu6ki3Mc/0NgC+JNKZtNG1x24rIPAuuE9eexG+WPgKsZexK292EZa0y2wy6Asv4c3bb+HBa8voc9/o7FRcAGJetdE1Ox6SeuKXHs+tijbjah4vPAw8D5wF9yHP9DYE6JesfiAmBRiePeCOyATUQ4EH6/Bb4PPJXj+B8CbwFuLFH36BS8gxTJPBU4ruSdahes0ZcD1wOvL3GnuyvinXNiKHNSyeOPw4JmP/BNYJMS/uwS0Z/P5G3HNRx/ExYknwXmAtMLHDu/St1rKPOJiv4sC33tYeCfgXEFj43tjwLfK3nslqFtlgG/LtpvQt1LOuDPtiWPPTT40w/8ANi8RN3zcuxXq8yzCn8GvgN8U1VvT21MBH6Kzdg4X1UbO3FaG+cB1wJfVdVlqY2JwFnAPdgsjWWeLOrEY9gN+hJVvSW1MRG4BfPny6p6f2pj2qll8FTVRcARqe2Ihar+ATgqtR2xUNWLU9sQE1U9KbUNsQg3s39NbUcsVPU+aupPXf/a7jiOU2s8eDqO45TAg6fjOE4JPHg6juOUwIOn4zhOCTodPCd2uHynHL3WLr2YBPRaG9Xyy54qdKTTicjGYXWvTpSfgH3C8h1JrYjH22BIOzWWMMZ5Zlgfl9aa6ojIjLC6fVJDIiEiO4TVToxaikLoQ31t67no1B37xLDcQUTWH3XPmhNO5hlh84xaCBJUIAhe7BQ2T05oSix2x4a9ArwrpSGR+FRYvkRENk9qSRzODMsjRWTCqHumY5+29f3yHhQ9eIZs5v3ZJjZ8r8nsi6k7EZb7JrQlBnOxdgF4b5Ozz3AjW0DLnzOanH2GrPMzYbMP+EJCcyoTss43hM3JwCEJzRmRtj6UsSBvgtSJzPPEtnLHAx9tavbZdmKnh3+aToGTWzdC1vlRWu+f+mh29rk78Mq27XVpdvb5KVpZtAAHNDz7PBOYEtanAqfVMPvcB2hPIDYkb/aZc3D9i7F3FhdiA+3nhN+EYfttTEs1KPutwsalRhMaGMXOtwa7FFNcmQO8skJ5+wX72/1ZCczpgi9T286zAp8L67MqlHnWCP6sBjbugj8bB/uvaetD+1FAhGOEMm8f5oti4h6lyyxQ97bBh8E2f3atUN7avFCKbgD4Qad9CfXPbutrq8P6NhXKe9MIsWAl8KEu+DJ+2LXz47C+0bD9BNM0GN6H/gLImPXkNGZ+cDxr3MWhYd82bL91gesw1Z2sEW4CPtKFE7ZJqPPZtqD9HHB9xQ5wEyZNpsAvw/aOXfDnXeH8LQ51Lw1t8NkKZX6EliJSJot3HbBuF/w5F1PKWtHWhwaBnSuU+dW2vpb58908HT+CP3djSj9Z8Fwa+tzUkuVNAX4ygj9zu+DL9NDXlrYF7X7gjxXK3BKTgcv8+WtY36cL/swO7ZJdOytCnz9n2H4CfGOEc/6tmMFzi3Ay26Pz34HxoxyjwDWdPlHD6ryBoSLMS4G3Ryh3WihvrS76MgV4Ztg57yeCKDJwFZFl1HLUuU1b4Mx+DwJ9EcouLCcWoc7D24JNdrP+bqSyFXigy/78gKFPJUuBQyL6c2AXfRmH6f6297UVwFZj2FioD+V656mmbHIlLSHgZcDxqro6z/Fd5FgswGQ8jMnBNQ5V7QdOwbJ9sDvn+aoaSxC5q6jqXcDPsZsbWB86RlUH13xUrbkQuyAzVtP6yqSJnID5kLEMExxuHGqygsdhPoD1uWtV9Z6Y9RT5g9EJDA2etTuxqvo7Wsr1q4GjNdxWGsp5WDYAdmc8NaEtMZhL6wJdhKnPNxJVXYW9h86C/+Wq+lA6i6qhppX547A5QD2ToyJcDjwd1lfTga9+cgfPkH1m0xLU+cQeG5ZLaWjWmdGWfYJdnI3MOjNC9vnHsNnkrDPjQlrBs8lZZ8YJYTlIDZOjIrRlnwC3x846ofinSt8My9qe2JB9AtzU8Kwz47ywPCepFfH4blg2NuvMCNnnA2H9obTWVEdbSu1/qXFyVITLw/KCThReNHj+HaAhJ/bx1AbEIGSfkG+yqybwGEAPZJ0Zi1Mb0AGeSW1ADLQ1pcrCTpTfi4IKjuM4HceDp+M4Tgk8eDqO45TAg6fjOE4JPHg6juOUoJbBU0Q2EJEXpbYjFiLyChGZnNqOGIixTVOVpYYjIhNEZKvUdsRCRGaKyEtS2xELEdlURNZKbcdI1Cp4hqB5NvAopr7TaERkaxG5CriX1sfujSQEzb2AO7AP3Rutqh+C5uFYX7tbRNZNbVMVQtA8Dfuc8Nax9q87IrKZiFwM3EfrW+daUYt5RURkA2yo24ewgD4h/BqJiGwNfBHYE5uLRqjJuS5KyDD3xHRNN8dEUp6lpTvZKIKe5CHAaZjs33RMfatWiUReRGQmcAxwFObDZJp97WyGiUD/E3bN9FHTa6fTRr1dRMqO8tm2wrEAMyocO5zsg+5nKzytflpEPl3BhpijpRYBVDi/awNX1ujJfQmwTgV/JgNP1cgfgE0r+LNhxWtnWoVj18QVFc7vuyr605GbSSeD50GYxmYeG/YAdsMytMxRxfQNy7AUkwyLgqr2i8hR5GuEtTCR3+2wu2aW0SwCbitpwr2YQlQsjsAku5aNtSOwEXAgJmY8HmujlcGmh0rULZgkXky2Bw4mnz/bAPtjAaJ9hsrrsAy0KOOAeSWOG405DFXIXxOCaVfuifWzdn/+q2Tdq4h47QQ+jN1wx2IK9jrojWE7i0/LgV+UrPsRyp+L0SmokzebDukmYjMgzsMugJXAnd3S/+uQP5sCF9ESzD0rtU0V/dkRe5e2HBM17riafgd96cPEph/CbrSrqaDQn/qH3QiOxTLwfmBhapsq+rMRJp7dH9rmsi7UWVjPU8KBuRCR2cCNqtqx55vwDufTwLOqumCM3WuPiGyKKe58X1WvT21PVURkR+CTwHHacDEMEekD3olleoer6vOJTaqEiEwDPgbMVNXPpbanKiKyEab0dJOqXtrhuiyCFohttQuejuM43aZM8GzkXxgdx3FS48HTcRynBB48HcdxSuDB03EcpwQePB3HcUrgwdNxHKfECL6iwXN7G9TgOI7TUxQObEWD5/Nxh1g7juPUgsITEvpju+M4Tgk8eDqO45TAg6fjOE4JPHg6juOUwIOn4zhOCWoZPEVkbRH5rIi8LLUtMRCRLUTkGBHJIwhbe0TkrSLyfhGp5fQIRRGRD4nI7qntiIGITBaR40Vky9S2xEBENhaRuSIyK7Utw6ll8AReA5wM3CMiF/RAEH0nNqfRYyJyUg8E0c8B3wIe7pEg+hXgahH5vYjsltiWqrwMOBW4Q0Su6IEg+lZMJP1vIjK/TkG0qJ7nkcDZeTTvRGR9oOyUrq/DJoFaC1OSHsCm5Pg34LEcxyvwqKquLln/CxCRLYANSx5+CPABbMqG57Bvyr4BfId8U0esUtVHS9Y9IiLyOmzyszKcC2wd1ldgPpyBTa8xkOP4paq6qGTdLyBMUrcL5T9Cvp7WFBb9wP3YBHG/yXn8IlVdWrLuFyAiM7BpXMrwcuDr2OR2A9j1cwMwH3gwZxmPqerKkvW/gJD8vLTk4XMwgeep2AwGCnwP64OLcxw/ADyiYwQ6ERkA+gppFReUqj+SnFL1wcmUvysiSvRPqIE/W0X055jU/kSeQuHBxP6siuxP6r72Pz3mz/tz2DhQtF92+rH9eFWVoj9gV2w+FrDMZinwWWCtnMcDxHxcycqcXNKf47EMACxLewT4F2B8AX/Wj+jPIEAZX4I9N7eVtRzL0HbNeWwnplaZWtGf5aGc1diTwSXAFgWO7cRri4tK+rIVdr0QfFkOnA6sW6CvbdEBf7Yr6c9htNpnGfAkNqHcpMT+1HM+5MB0rBN8EThHVfM83taZPixoHodNaJXn8bauKHan/j1wjKrektieqij2SPifwAmqen9ie6oyDQs4ZwMLVPWZxPZUZTIWNE8ALtSIrxSqUNfg+Rvgn4H/7oGgCfBtbKre/2p40Mz4V2AdVb01tSGROAB4qAeCJsB9wLuBX/RA0AS4AngWu3ZqETQzahk8w0m6PLUdsVD748iPU9sRC1X9U2obYqKqP09tQyzCzfmK1HbEQlWfpab+1PVTJcdxnFrjwdNxHKcEHjwdx3FK4MHTcRynBB48HcdxStDp4LlOh8vvFhOHLZtOr7RLxoTUBnSAGakNiIGIZDFmraSGdICOBM8wbhpg306Un4DDhi2bzhwY0k6NJYiSzAjr0xKbU5k2EZwdkxoSjwPD8gMpjRgNEZlEiIUiMjnvcZ3KPM8Iyy1FpKw4SC0IJ/aksHli2G4soT2yNjkzpS2ROJjW8NlPpDQkEvPCch0R2TmpJRUJWef8sPleEZmZ0JzR+GDb+ofyHhQ9eIZsJmv0Pmx4ZZM5HMgC5qSw3WROx9SdAHZqcvYZss4v0Qqec5ucfYas891hcxytwNNUDgQyCbk+TJCmVoRkaF7bP52SN/vsROZ5Bq1g0we8ranZZ9uJzS7IacC8pmafoR32ptXuU2h29nkwQ2X1xtPs7HMerRsbwHZNzT7bss6sfSYBR9Uw+/wgNnY+YyJ5s0/NJyn1RuAs4NeYiMJZ4bfesP22DP+/mpZ4xABwcZ56qvywC+cLwS4FVoX1/SqUeWjwZVVbmauBw7rgz0Zt51mxIWpnAdtWKPP7bW3S3k7R5O5GqXvnYP89bX1oAbB2hTL/1uZL1t+WYGpVnfbnvW1tk/kzl6CRW6K8DWjpbyqmfDUA3NAFXwQTuB7uz79UKPMdwf7s2lkdfp/rgj8zh107fwzrO43g9xMj9KFFedoxrzHz2ho1+w0CuwzbbwLwGUw5SIGFYX2HLpywFw87CdnvqgplbhLsPzGUdWLY3qQL/uw7gj+rgaMqlPmmYP8jobzjgKOBCV3w58sj+DMIvKFCmQe29bXMn/d12pdQ9+0j9LUVwNSS5QmWNWf+rAzre3TBl+mYCPRwf35bocx1g/2ZP5eG9Vd2wZ/tQ98afu2cMcK+7xmhDx2Uq56cxqyDaem1G3P7aNE57HNNNzpyW53fCZ2uvTNvF6HcaaG8tbroSx+mkNN+zheXvTiHlX0VBYVfI9S5QWiP9sB5S6SyCwvZRqjzrcOuiRXAKZHKVuCBLvtz2rD2WRYrcIfyDuyyP7eNcGNbfwwbC/WhXO881aStzg6BCUwr8GgNtdaIk2hNA6HATap6Z0J7SqOqg9gL9kySrx/4gqquSGdVeVT1SeCb2GMcmD+1+wNCAa4DHmjbHsSy66YynyCSHfgLNn1HUzkGC5hgfe7rqvpUzAqK/MFoAa2TW8sTq6qPAD8ImwNYCt5krsTeyYA9dnwtoS0xOA27qQHcrqq/TmlMFULicDSta+IsVV0yyiG1RlUXYwkSmE91TI5yo6Y1myVOSge++skdPEP2+cOwWecTe1JYLmxq1pnRln0C/HtTs86MkH1eHzabnHVmXIcp0EOzs86M7NOofmqYHJUg62P/HTvrhOKfKl0SlrU9sSH7BGj61BAZV4bl95JaEY/LAJqcdWaEBOKvYb2xWWdGyD4B7qlxcpQbbc10cGknyi8aPJfBPzpN3Yk2FWxKQvYJNplXL9AL06q00yvt0k5/agMi05E+56pKjuM4JfDg6TiOUwIPno7jOCXw4Ok4jlMCD56O4zglqF3wFGMvEdkjtS0xEJEJInKIiGyR2pYYiMhMEfmoiMxIbUsMRGQzEfmgiIwbe+/6IyJvFpH9UtsRAxEZJyLvEZHXpLZlJGoTPLOgCdwBXA38W1KDKhKC5geBR4FvA0ckNqkSIWiehvlzDjA7sUmVCEHzYuBPwDdo+NQkIWj+CvgFcF5qe6oQguZB2PDXizDFp9oxvsPl9+XUvnwLJmq7GS3tzHEiUrd5Tybm8Gc8pjN5KuZL5s/EuvmTs21mAJ/GFH8E0wB9FphS0p/c0xwUJac/mwInA/tjbTUe+1ZzLRF5fpTj1kTHEpCc/uyIDT3cDpia2VS3vgaMz+FPH/BOTLB7HVpaoBNq6E/Hg+fbKP8R8auxi7QsZS6ENbE6LBdVKOPj4VeWmP7cDuxH+bZZm9ZoszpwP/AKyvszmaEiH0UZHHuXQigW5Mv68yKqXTvLKxy7JqqM8jmA+sSCfyBFBguJyGzgRlWVHPuOY6gq9mjsjimab04rU7sb0+Urw6CqRh0lISITaE33MBrtmed0Wv6cg4nllmG1qkbrACIi5J9xcgbwSeAoLDPIMs8PAz8pacJKVV019m75CKrleROBl2OZ5z8xNPPcivI3xxUxR90VvHbehD21vRZrGwEeB0q/Y1fVqMEzTJeSJ0MXTKP1S5geaJZ5/hgTny5DrmtHREyXLkds+wcFNfJm0yHdROzEvQV757kSuLOb+n8d8GcCpkT/OJa5npXapor+zMREsZdhEl9zUttU0Z9Nsfdp/ZgC16zUNlX0Z0fgVizLWpjanoq+9GFzOT0U+tplXaizsJ5nxzLPsoSsaA9gQFVv7FQ93SJkrO8BblPV+1PbU5XwV/b3AJdqD4hhiMjLgd2A76nqwOh71x8R2RFYV1WvTm1LVcITxYHAvar6xw7XVTjzrF3wdBzH6TZlgmdtPlVyHMdpEh48HcdxSuDB03EcpwQePB3HcUrgwdNxHKcERUcY1W6IlOM4TgRWkW8QzD8oGjy3LLi/4zhOExhPweBZ9LE92pA6x3GcGlF4eK2/83QcxymBB0/HcZwSePB0HMcpgQdPx3GcEnjwdBzHKUEtg6eIbCIi14vIO4MsVaMRkdki8lMRKSvuXCtE5GgRuUBEXpbalhiIyEUicpKIrJ3alqqIyDoi8jMReV8QIW40IvI6EblORHZLbctw6hqYNgXeDFwA3NcDQXQHYC/gBhG5qQeC6H6Ysvc9PRJEDwCOBx7rgSC6ASZa/jXgYRF5f8OD6Gsxf34iIr+vUxAtqud5JHB2Hs07kb6jQQ+j3Pwu62KdIJuKYBWm+H0zNnvjWKwDHKGqz5So+wUEgeYLgdeXLGIzbJ6c7LytBp4CbgIW5zi+T1U/UrLuFyAi07AZSmeVLGLrtvVB7Bu5+4BfYrMAjMVjqvqFknW/ABHZGriY/FNXDKfdn9WYP3cA/5Pj2HHATar6/ZJ1vwCRvg+AHku5a2c60H4zWxV+t5JvnqaZwAkxhbtF5KuYwHkZNsGmssluAKuxa+ZG8k2bMgU4XFVXj7aTiAxg11nHxJALBE+JNqdLSW5T1Z1iFCQiE+nQJFIFeJOq/jZGQSIyD/h8jLLKElNQW0QWAevFKq8Mkf1Jfe08oqovjVVYDfw5XlXPGG2HMsGz04/Cx6uqFP0Bu9KaLW858Btg1wLHg83wGJvJJf05HpsnZzU22dglwBYF/ckzDW1eloBd8CX9uTmU8xzWPqdjUz/kOXZBRD8yBir6k014tgx4EpvcblLBY2NzUUlftgp+EGy7C9ibEBhy9rUZHfBnu5L+HAaswLLwfmzSwW0LXjvT1mBTJer6LuQ5LPj9BjhGVW9JbE9VnsNuVJcS+ZEoEcsxn74MLNBIr0cSkt0ETgAuVNU8rx7qykosWNwFHA38TLXA42X9eA6bj/4nwFxVvSuxPf+glsFTVX8nIi9R1TzvN5vAOcDlqrowtSGReBcwsQeCZsargcUND5oAqOqDIvISbAbNJgfNjEuAm1X1sdSGDKeWwROghwInqjoI9ErgRG1e7049rnYdVX0ytQ0xqWOgKUu4AdTSnyZ//uM4jpMMD56O4zgl8ODpOI5TAg+ejuM4JfDg6TiOU4KOBM8wnBGGDhNrMtnolbLDGevGy2BIOzWdyakN6AAbpTYgBiIyNaxuktSQsSl8LXQq89w3LPfvkQv02LA8JqkVEQjtcWDY3He0fZuAiMwgzOoqIhsnNqcyIrJDWN1BRMqO1a8TR4TlJ5JaMQoiMosQPEVk/bzHRQ+e4eLMhuDNAnaPXUc3CSc2E+X4cJGTW1N2BzIfFvTAze1TbesnpzIiIvPDciI2GKGxhKwz01DYQ0Q2T2nPKMxtW/9s3oM6kXnuS+uRYyLNv0Dn0krp+xh6ohtF241tYvinjWhw9hmyzs/Qap/3Njn7DFlnptw1Hjij4dnnEbQG4owHoilpxSIkR0e0/dNH8iZIUYNn28U5ve2fX0lDs89wYj9K653aJOCjDc4+d8faI2M6zb65fYqhMnR9wEmJbInBmZiEWsa6NDT7bMs6M1GOccABNcw+j2fo+84+8mafqjrmD3g/cCfwNKZ1eCfwB2CrYfttHf5/cViuxAb2X56nnio/TAnm18E2bbPzCxXK/DCmhbi4za9VwEe74M8bgv2ZP38L6wdUKPPS0B7PD2unbbrgz0eC/cvb2uYOYNMKZT7aVp4CS8NvfBf8+foIfe3n2Jj/MuW9KJSzJCwHCBq2XfBlEnDDCP78R4Uy96elvZm1zUrglC74s3mIT3e21X0npuvZvp8A/4upUGV+LwtxTsaqJ+/Y9klYxpLJom0bGnc4dwM7YSn6zcCfgCOBv+aspwqDwDYMlZ96FSYCW5aLgD9j2cC1wEHYxZpHJLcqfZj9WRu9FOt8VfQIjgL+HVND+j+YIvwA1k6dZgqwJa1XBtvSEh4uyx5Y0Mkk8t4BPKNjCN9GYmOsv2VsCzxOOQFjVPUJEXkjdp5uxoLOAdhNs9MMYuIoG7T929ZAFfWva2g9cd6MPZH+AvhjhTKLsDWtp5LpDI1fgI2bF5Gdsb/NZH3o7cAiDdF1VHJG8onAE7Si82rgR2Mco8A1nb7LDKvzBEz7L7OzH9gwQrnTQnlrddmfG2gptSvwIKbLWLXcq0Lf6aYvU4Bn2nxZBXw3UtkDCfzZZlhfWwocEqlsBR7osj+HBx8yf1Yw7Mmyoj8HdtmfH4Q+lvnzv5hG62g2FupDud55qkl1nRCCEVgGlPuvUl3kbFoZ8QBwvqr+PaE9VTmW1jlfgWmblspsUqOq/cAptBT5VwMnprOoGmq6kj+nlTkvA6JNxZGAC7E+BubTtap6T0J7qnIC1sfAXlWdqKpRZ4Mo8gejC2l1/J/W8cSq6lKgXW7/1FS2xEBVfwdkU288DfwooTkxOI9Wh75cVR9KaEsM2r+8OF6787qgI6jqKuBz2SYN/qoEQE1w/MdhcxXwrdh15A6eIfvMDKhj1plxdlje0/CsMyP7QP+LTc06M0L2eVnYbGzWmRGyzyVhs8lZZ8aFYflMHZOjEpwQlhfHzjqh+KdK/w1Q5xMbsk+A25IaEomQfQL8NKkh8bgWoAeyzoz7AZqcdWaE7BPgL0kNiYS2prvpyLVTNHiO9Bf2utL4zjyMXphSAXqvXZp0TeSl13zqiD+uquQ4jlMCD56O4zgl8ODpOI5TAg+ejuM4JfDg6TiOU4LaBU8RmSkip4nIcaltiYGIbCYiF4rIXqltiYGIvFlELheRTVPbUhURGSci7xGRi0Vk0thH1BsRmSYix4vIF1PbEgMR2VhEvi4iB6W2ZSRqEzyzoIkp5RwDHJzYpEqEoHkxJrrxPmDvxCZVIgTNXwE/A+ZgQhiNJATNg4AHsIEf7yao0TeRLGgCCzFB6EPSWlSNEDTPA+7DxtwfOMYhSaii0JOHzXJkkOOBtwBvxoJ5ZtO2InJNJ40rwSdFZOUY+6wF7IMJR/TRukF9QERe1UnjCpBNOZAnu98IU/fZEJgQ/m0l8EUR+WiJul9U4phc5PRna0xNaiothSeAy0TkuRLVblnimDzMyOGPALMx9aI+Wu2zYcVrpxNJ1QEi8oox9pmCqRq9HvMtU0Xat4axIJ+qUpvyyGxyKo9gghaa5icKHB1RoWVcOl/+8Vs3oj/vSezLYGQFnR8l9Uf67o/sz6J0vogC50b2J/W188a8Nhbyq+BJmF20gpzlrg98BVN1eS44cmfserr1w3Q4fxz8yTQrz0ptV0lfBHsyuIOWaOwSYE5q20r6MwE4FNPezCTY+oFZqW0r6c9MYF5om0wib2Fquyr4symmo9tPS1Lusi7UWzh41uKdp6o+paqfwqbEPRd7LBzr8bi2qOqfVfUATA3+J1jDrBr1oJqixs+B12Hq4H8A1qahwyxVdZWqng+8BBPqfgITyW2k6IqqLlbVE7GpfRdgQafJ186DqvpeTJz5UqxdanntSIi6+XYWmQ3cqKodnfMmzBEkqvpkJ+vpFiKyBfCoqpZ5p1YrwnxHrwbu1iKdp6aIyARgM1W9N7UtMRCRmcB0VX00tS0xEJGXY8ruyzpcj6WfBWJbLYOn4zhONykTPGvx2O44jtM0PHg6juOUoJbBM3zE/LZeGPUBICJri8geDZ4ffQgisqmIvCa1HbEQkTeJyEap7YiBGHuJyLSx964/IjJFRN4qIrWLVbUzKLATNnXpYyJyRA8E0SOA64G7RWSfHgii5wO3i8i1IvLa1MZE4HrgARE5V0Q2TG1MRV6JqfUvFJFjeiCIHoTNYPGAiLy7TkG0NoaMwFJgPWA+vRFEB4CtgEtofhAVrO/sCfyqB4KoYJ8rHUZvBNFl2OdkJ2NB9NiGB9F+7DPGb1OjIFrqr+0ds8ZxHCcdz6rqjLw7Fw2egn0o/dIShhVhS+AD2PjjQeyD7N9hU9c+HLEewcbRdvqD7/2w8cfjsA+YFbgKG0mxZJTjijIVG1MffabAYRwLbBzWn8M+Yr4A8ynmB80bYWPhO/0B+9lt6/3YNM9fA34ZuZ4tgE5ngBsBn8QyacX625+wwSd/jVzXdgSdhA4yG3gHMBm7TgeB67BXR4si1jMeG9X4VO4jUg/HWsNQqV2wRl8BXAm8OrVNFf05LjT8CmwY6vqpbaroz02hfZ7GRulMTm1TRX+WAcux4LI/Ialo4g9LPLJr5+fAG1LbVNGfQ7Ebcj+mgLVJapuyX6HMs1uIyLrAXOB8Vb07tT1VEZGtMYm9r2iRO1tNEZEDseGN52lvjJo6FrgXuErreEEUILzbPBEbD/4/qe2pShhhdARwjtZs1FQtg6fjOE7dSf4XK8dxnCbiwdNxHKcEHjwdx3FK4MHTcRynBB48HcdxSuDB03EcpwQePB3HcUrgwdNxHKcE/x9DZzn1jEJzmwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "plt.axis(\"off\")\n",
    "ax.set_aspect('equal')\n",
    "\n",
    "def draw_invP(x):\n",
    "    a = 0.2; oma = 1.0 - a; l = 1.0 - 2*a;\n",
    "    draw_U([x[0]+oma,x[1]+a],[-l,0.0])\n",
    "    draw_U([x[0]+a,x[1]+oma],[l,0.0])\n",
    "    draw_U([x[0]+oma,x[1]+oma],[0.0,-l])\n",
    "    draw_U([x[0]+a,x[1]+a],[0.0,l])\n",
    "\n",
    "nx = 6\n",
    "ny = 4\n",
    "plt.xlim(-0.1,nx + 0.1)\n",
    "plt.ylim(-0.1,ny + 0.1)\n",
    "\n",
    "for x in range(nx):\n",
    "    draw_U([x,0],[1,0])\n",
    "    draw_U([x+1,ny],[-1,0])\n",
    "for y in range(ny):\n",
    "    draw_U([0,y+1],[0,-1])\n",
    "    draw_U([nx,y],[0,1])\n",
    "    \n",
    "for x in range(nx):\n",
    "    for y in range(ny):\n",
    "        draw_invP([x,y])\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae8eef2e",
   "metadata": {},
   "source": [
    "In this configuration, we only have to consider integrals of the type\n",
    "$$\n",
    "\\begin{align}\n",
    "N(X)_{ab} = \\int d[U_\\mu(x)] (U_\\mu(x) X U_\\mu(x)^\\dagger)_{ab}\n",
    "\\end{align}\n",
    "$$\n",
    "for which we can show that for every group element $V$, we have\n",
    "$$\n",
    "\\begin{align}\n",
    "(V N(X))_{ab} = \\int d[U_\\mu(x)] (V U_\\mu(x) X U_\\mu(x)^\\dagger)_{ab} = \\int d[U_\\mu(x)] (U_\\mu(x) X U_\\mu(x)^\\dagger V)_{ab} = (N(X) V)_{ab} \\,.\n",
    "\\end{align}\n",
    "$$\n",
    "We therefore have that $N(X)$ commutes with every group element $V$ and therefore by Schur's lemma\n",
    "it must be proportional to the identity matrix, i.e.,\n",
    "$$\n",
    "\\begin{align}\n",
    "N(X)_{ab} = \\int d[U_\\mu(x)] (U_\\mu(x) X U_\\mu(x)^\\dagger)_{ab} = n(X) \\delta_{ab} \\,.\n",
    "\\end{align}\n",
    "$$\n",
    "We also have that $N(X)$ must be linear in $X$ such that we must be able to write\n",
    "$$\n",
    "n(X) = n_{ab} X_{ba} \\,.\n",
    "$$\n",
    "Since this needs to hold for arbitrary $X$, we now have the identity\n",
    "$$\n",
    "\\begin{align}\n",
    "\\int d[U_\\mu(x)] (U_\\mu(x))_{ab} (U_\\mu(x)^\\dagger)_{cd} = n_{bc} \\delta_{ad} \\,.\n",
    "\\end{align}\n",
    "$$\n",
    "We now determine $n_{bc}$ by setting $a=d$ and summing over them such that we have\n",
    "$$\n",
    "\\begin{align}\n",
    "\\delta_{bc} = n_{bc} N_c \\,.\n",
    "\\end{align}\n",
    "$$\n",
    "We therefore finally arrive at the most important integration formula for this chapter\n",
    "$$\n",
    "\\begin{align}\n",
    "\\int d[U_\\mu(x)] (U_\\mu(x))_{ab} (U_\\mu(x)^\\dagger)_{cd} = \\frac1{N_c} \\delta_{bc} \\delta_{ad} \\,.\n",
    "\\end{align}\n",
    "$$\n",
    "Let us first verify this numerically:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "id": "578077fa",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "checking all a=0, b=0\n",
      "checking all a=0, b=1\n",
      "checking all a=0, b=2\n",
      "checking all a=1, b=0\n",
      "checking all a=1, b=1\n",
      "checking all a=1, b=2\n",
      "checking all a=2, b=0\n",
      "checking all a=2, b=1\n",
      "checking all a=2, b=2\n"
     ]
    }
   ],
   "source": [
    "np.random.seed(13)\n",
    "for a in range(3):\n",
    "    for b in range(3):\n",
    "        print(f\"checking all a={a}, b={b}\")\n",
    "        for c in range(3):\n",
    "            for d in range(3):\n",
    "                v, e = haar_integrate(lambda u: u[a,b]*u.H[c,d], U, 2000, 2)\n",
    "                err = abs(v - (1.0/3.0 if b == c and a == d else 0.0))\n",
    "                assert err < 3*e"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a73108a4",
   "metadata": {},
   "source": [
    "The lowest-order $\\beta$ contribution for a $n_s \\times n_t$ Wilson loop is therefore given by\n",
    "$$\n",
    "\\begin{align}\n",
    "\\langle W_C \\rangle\n",
    "&= \n",
    "\\int d[U]  \\textrm{Tr}\\left[\\prod_{\\mu_i,x_i \\in C} U_{\\mu_i}(x_i)\\right] \n",
    "\\frac1{(n_s n_t)!} \\left(\\frac{\\beta}{2N_c}\\right)^{n_s n_t} (\\sum_P \\textrm{Tr}\\: U_P)^{n_s n_t} ( 1 + O(\\beta)) \\,,\n",
    "\\end{align}\n",
    "$$\n",
    "We now need to count how many link pairs we have in this case.  We first note that we have\n",
    "$(n_s n_t)!$ permutations of the individual plaquettes.  Let us now first integrate over the internal plaquette\n",
    "pairs.  We can write such a pair as\n",
    "$$\n",
    "\\int d[U_\\mu(x)] \\textrm{Tr}\\: [A U_\\mu(x)] \\textrm{Tr}\\: [B U_\\mu(x)^\\dagger] = \\frac1{N_c} \\textrm{Tr} \\: [B A] \\,. \n",
    "$$\n",
    "Therefore for each removed link we get a factor of $\\frac1{N_c}$ and obtain a joint loop with the integrated link removed.  This could look like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "id": "fd648c1b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAU8AAADnCAYAAACe9dTFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAuoElEQVR4nO2debgdVZW333UzD5CAEIWgCIKgART9REEwSIOtIoHGRulGRRBFbAUHpigyRFEhQT+kW8ChBUVEBlsQUATFAIJT04iIoEwiiQj5ICHDDSS56/tj7eKce7m5t4Z9zq46vd7nOU+dOk/tvdeq2rVqVZ3avy2qiuM4jlOMvtQGOI7jNBEPno7jOCXw4Ok4jlMCD56O4zgl8ODpOI5TAg+ejuM4JfDg6TiOU4JaBk8R2V1EXpfajhiIyBgR+WcReVFqW2IgIlNF5N0iMjW1LTEQkZkicrCI1PJcKIqI7Cwie6W2IwZi7Cci26W2ZThq1WFC0LwNuBH4amp7qhCC5sHAA8AlwDGJTapECJpzgcXAN4FGn6AhaH4VuA+4GNg4sUmVCEHzeuBXwEWp7alCCJpzgHuB/wI+k9ikYRlbtICIiEYeliQiuwPzgZ2AyTHrHqVdAYjpj4iMAQ4CzsBOyK5laB06NlOBjwBzgTHY8VkWs40R2u6EPzOBU4B3Y8nDeGB1zDZGaLsT/uwMnAnsBkwCJGb9o7Qd1Z9wPu4HLAA2o4vnThmkiO8isieWFTq1RAAfbus45ZCHVQe2zLt10cxzp7DcpGC50XgzMA94PjAl/HYf8JbI7QxlDLAUWBuxzglYVnMiMI6WPxcAp0dsZxh0DLAkcqXPAz4KvBfL1CYAT2H+XR+5raEodnxisi3waeyxwzisDzwNvAF4InJbQ1kLLI9c5xuwc2drWndtS4BdI7czFAFWAM9ErHMsdtf2aey8yc6d64APR2xnOP4EWux/CVXN/QGOJtzlxv5gB2OOOcHTwJ2daKdbH2Aidrv7BHbSnJXapor+bAp8EVgJrAHmpLapoj8vw56nrQLWAZuktqmCLwL8A/A/2COIxaltqujPOOAw4NHQ1y7tQpvrisa2orftRwNnq2rHnqu0PfdYq6rXdqqdbiEiE4FDgZtU9Y+p7amKiGyKZdYXqGqnM7WOIyIvwwLPuaq6LrU9VQjnzl7A81T10tT2VEVExgGHAH9Q1d90uK11QF+R2Fa74Ok4jtNtygTPWr2q5DiO0xQ8eDqO45SgtsFTRLr2vmc36CV/RGSsiExIbUcsRGRy9s5vL9BLfQ3q608tg6eIvAJYJiJX1HVoVhFE5P3AEyIyX0Riv+aVgu8Ci0XkQz0SRO8B7haRfZseRMMw4CdF5DoReWVqe6oiIgcAS0XkPBHZPLU97dQyeAIbAv3A/sAdPRBEN8beYfsw8JceCKIzMJ/OBBb1QBDdGNgeG0bb9CA6CXu9Z2/g1h4Iohthr/odBtxfqyBa8F2o3O95In0PYi85p/pEew+R1ovUKf2ZEtGfvREZSOjLo5Hf0fsi0pfw2MgNUf2RvtvT9jX5cNzjI8vT+sNWOfpQfd7zFBHFAk6Z0TsTsICVtZOd6A+SbxTILsAfVfXlJdp+DiIyHvNlZckqpgxZH8BGZtyHvaQ9GrsAb1DVm0u2PwgROQ7LGmP6swLzZ7Tj/UJgszx9KC8isgQbCRXDn+yEWIL1t9F4NTAmsj+KZY9lRu+Mw8brZwyE5cPAYznK7wIsU9XpJdoeluBPf5stRZiCHZP2WLAWuJ98o7V2AU5X1ZNGsbHwq0qdyzzN4RNKZhJ7YEFlBRYsPwJMLFBegbviXTkZH+qcULL8CdiVbTk2gmoO4R3bAv7sEdGfT+Q9juspvxDrxKuAG4BXFSg7v0rb66nz7xX9yYYZ9gNfA7YoWDa2Pwp8u2TZ7YIfK7FhtHOBqQXbXtYBf3YsWfYwLHFZjl0A/gW7WBVpe16O7QpnnoVVlbrEIuyqfx7wNVXtiupNB7kL+ANwEvBDDUerwdyKZUbHq+rtqY2JwG3AQ8BpqvpIYluq8gQtmb1zVHVFYnuqci/wZ0wX4lKt0SiwWgZPVX0AmJXajlio6jXANantiIWqzk1tQ0xUdZ/UNsRCVR8HdkxtRyxU9VZgh9R2DEdd/213HMepNR48HcdxSuDB03EcpwQePB3HcUrgwdNxHKcEHjwdx3FK0JHgKSLZiI1Oz6PSLV4xZNl0doNnZ8bsBaalNiAWIpK9PtgTfU1Etgxf6x4LCsfCTmWeHwnLfRouGJHxubD8fFIrIhCOx5vDaqcn1eo4IrI9Npw3m8K66RwSli8TkekpDYnEvLA8tq5iK2H65uz7q/KWix48Q9aZvUQ9FjgidhvdJOzY14fVXdt3dEM5gtbgiLltdwlN5XNt3+cnsyICIev8QrYKHJfQnMqErPMdYXUmUNfBCGe0fT8zb6FOZJ4foXVyjgfmNTz7PIOQ2WAzYubeuXUjHId5tIQjxtK6S2gcIet8c9tPOzU8+zwEyB6ljAGOaXj2OQ/zA2xa5AV1yz5DMtTeZ3bNm33mCp4iMi7oT24Y1jcJHxmy3RTgk7TmjwbTF+xK9ikiG7XpZI4PNk6qUN/O2LzY2X4SYI9uZJ9ibNLmz0ZhvcqQ2iOw45ExGfhkN7JPERkffJkS1jcRkedVrPZzDFYQmkyXsk8RmRL86Qvrm4hI6WevbVln+3PoMXQp+xSRaW19bUzwp7SCe8g6D8ZUnjK2oUvZ55BzZ2pYHz/MpmdgSVHGJHImSHkzz9Mw5ZrTwvpfMXmrPYdsNwNTJ2mXBhvAhGY7iojMBB4PtgFsC/wNuKJCtVtj8leZbNyqsL5thTrzsj+2jzN/LseOwUcr1Lkdg2XBVob1GRXqzMsCbB7uo8L6X4HHReS1FercApunPGMlUPUCk5dfAItpXYweAR6tEHCmYcGy/dxZRxf6Wvjj8FHMB7AL3GLgpgrVbokpVWX+rMP62ksr1JkLEdmNwefORzH/vjBkOwFewGBZyFXAZrky5JyyTntgO6FdYHQ5MG0UKahrY0pbjWJjH6Ym027jCuDICHVnmoIbdNGfGeFAtvuzCtgpQt1XEVlGLUebbwrHo92fpcDkCHUXlhOL0ObpQ47PAHBzpLoVeKDL/tw2TF87NaI/B3bRlynAsmFiwV6j2FioD+XKPNVEeH/f9lM/cJaqLstTvhuo6gB2i9MuwdUPfDONRdVQ1ccwbcmnw08DwEJVvTOdVZW4HnigbX0V8BlVzSMGXUfmMziL76fZf/Acx+AMbAD4UiJbKqGqK7E3Y9r9+RNwY8x2ivxhdBytW6R11HPHXond2oLpTX5aVcuocdeFz9E6QZ/BRJUbidrl/VhafWgNcG46i6qhqkuBs2kdn9tV9ZfpLKqGqt4CZBfmAWqWHJXgHCxOgfW5Y0MfjEbu4BmyzyxzqOWObcs+wXbcfyY0pzKq+nfg62H1Vw3OOjOupzUVxGcbnHVmLGj73uSsM6Pdhy8msyICbdkn2PPbqFknFH9V6SthWcesM+PKsLyt4VlnxulhecaIWzWAcOXPHqN8ZaRtm4CqPomdmDQ568wI2SfAw3VMjkpwTlh+PXbWCcWD510Add6xIfsEk+9vPCH7BJtbvBf4PUAPZJ0Zi1Mb0AGaPhUJANqaguSuTtTvwiCO4zgl8ODpOI5TAg+ejuM4JfDg6TiOUwIPno7jOCWoZfAMohhjRt+yGXRprHXXcH/qSy+eO3VTYsqoVfAMB/4AbCjVVYnNqYyITBSRY4DHRORjqe2piojMEJGzgaUiUndl8FERkVkichUmULJhanuqEM6dfYA7sHHqjSYouR0BLKKmOq21uOKGK8v+2E56ASbL1Z/UqAqIyETgSOBUTJJrCqYA1EhEZAbwKeD92AV3LbBpUqMqICKzsNEne2OSdusYLG3XGMK5szc22uklWF/7W1KjKiAi44BDscEhU8LnRUmNWg+dDp5vEZGyb/bvWKEsBO3RyDxV4Q7i4yLy8Qptx7x1UYAK+3cCcGWN7qbWABMq+DMGyz4jmlSZrSr4s1nFc2eDCmXXxxUV9u9BFf2ZOPomxelk8LwSm2hsNKcFmM5g0dSMgWF+y4GsAI02NldVV4rI3Vi2NZo/YzF/hOcGvJL+sIi4I1k+hE3IladDTsAuRMP1/LLHZ2G5cuvlCGwoXp7huFMZLNadkUmXFWUA6Vsw+maFuBDYl3z7dzrDZ81lj81q0Lmjb1eIX2HauKPt3zGYP31EO3fkcdDo49qBfHqebZp3s+mQbiJ2Mt9CSzf0zm7p/3XAlz7gIOAhTPdUMTGV5LaV9GcKcDymkbgSk/qak9quCv5sjik69WOSf88Am6S2q4I/rwB+HI7LALA4tU0VfBHswvHHtnPn0i602xk9z26gqrep6u7Y85tfMFg/tFGo6oCqXoZdbQ/H1Kj+mNaq8qjqSlU9E9gMe477BA0e/6yqi1X1KOz4fAPzZfXIpeqLqv5OVd+MJSA/oUNjubuBGtcALwfeCdxNTWOBhKibb2OR2cDPVbVWD4ccx3GqkD1TLRLbapN5Oo7jNAkPno7jOCXw4Ok4jlMCD56O4zgl8ODpOI5TAg+ejuM4JSg6wmhWR6xwHMdJy2oKJpNFg2cthEQcx3EiM56CwdNv2x3HcUrgwdNxHKcEHjwdx3FK4MHTcRynBB48HcdxSlDL4CkifSKye69MzCUik0XktantiIWIzBSR7VLbEQsR2UlENkltRwzCXEa7iciE1LbEQETGi8jr6zgJXC2DJ6ZAfxPwsIi8pweC6IeB20TktyKyZ2JbYvBt4G4RuaJHgugtwF9FZH4PBNFtMX8WiciHeiCI/itwM9bf9q1TEK1r8BTgKUx89z9ofhAVbJKxVwNX90AQHYP1nf2BO3ogiPZh89x8GPhLw4OoACuA5wFn0vwgKphC/vbAJdQpiBaUqj+anFL12HQAmvDzmYgS/WMT+6LAZhH9eU9qfyJPofBTTMw21efvkf1ZmciP7HNRZH+S9jVgzxw2rivaLzuZyQl26/3rEmW3wuYxmYg5vwa4F7gA+EuO8pcDbwc+XaLt4cgy9C9hO7kobwReiWVsazGfFmJX0qU5yl8ObEO8KWWzaYPLTlz2PmCj8P1p7EL5feAaRp+E7d1YxhqTnbAToKw/x7Z9X40Fr29jt7+jcREwo2S762MyNv3EtSXKbord6mYTKj4NPAx8E/hTjvKXA3NKtDsaFwBLSpR7DfA6bCLCdeHza+A7wOM5yl8O/APw8xJtj0zBK0iRzFOBE0peqfbADvpK4AbgVSWudHdFvHKOD3VOKFn+BCxo9gNfA7Yo4c8eEf35RN7juJ7yC7Eg+RQwF5haoOz8Km2vp86/V/RnRehrDwP/AowpWDa2Pwp8u2TZ7cKxWQH8smi/CW0v64A/O5Yse1jwpx/4LvCSEm3Py7FdrTLPKvwR+E/ga6p6e2pjIvBjbMbG+ara2InT2jgfuA74sqquSG1MBM4C7sFmaSxzZ1EnFmEX6EtU9ebUxkTgZsyfL6rq/amNaaeWwVNVlwBHpbYjFqr6O+CY1HbEQlUvTm1DTFT1lNQ2xCJczP4ttR2xUNX7qKk/df233XEcp9Z48HQcxymBB0/HcZwSePB0HMcpgQdPx3GcEnQ6eI7vcP1OOXrtuPRiEtBrx6iWb/ZUoSOdTkRmhq/7dKL+BOwblm9NakU8/hEGHafGEsY4Tw/fx6S1pjoiMi183SWpIZEQkdeFr50YtRSF0If62r7nolNX7JPD8nUisumIW9acsDPPCKtn1EKQoAJB8GK3sHpqQlNi8UZs2CvAQSkNicTHwvKFIvKSpJbE4cywPFpExo24ZTr2bfu+X95C0YNnyGbek61iw/eazNswdSfC8m0JbYnBXOy4ALyrydlnuJAtoOXPGU3OPkPW+Ymw2gd8NqE5lQlZ56vD6kTg0ITmDEtbH8pYkDdB6kTmeXJbvWOBDzY1+2zbsVPDT1MpsHPrRsg6P0jr+VMfzc4+3wi8tG19Y5qdfX6MVhYtwAENzz7PBCaF75OB02uYfe4LtCcQm5E3+8w5uP4F2DOLC7GB9nPCZ9yQ7WbSUg3KPmuwcanRhAZGsPNNwS7FFFfmAC+tUN9+wf52f54B5nTBl8lt+1mBT4Xvm1So86xh/FkLzOyCPzOD/de29aH9KCDCMUydtw/xRTFxj9J1Fmh7x+DDQJs/b6hQ34Y8V4puHfDdTvsS2p/d1tfWhu87VKjvtcPEgmeA93fBl7FDzp0fhO+bD9lOME2DoX3oT4CM2k5OY+YHx7ODuzQc2H8cst3GwPWY6k52EBYCR3Zhh20R2nyqLWivBm6o2AEWYtJkCvwirO/aBX8OCvtvaWh7eTgGn6xQ55G0FJEyWbzrgY274M+5mFLWqrY+NADsXqHOL7f1tcyfb+Xp+BH8uRtT+smC5/LQ5yaXrG8ScPUw/sztgi9TQ19b3ha0+4HfV6hzO0wGLvPnz+H7vl3wZ3Y4Ltm5syr0+XOGbCfAV4fZ51+PGTy3CTuzPTr/DRg7QhkFru30jhrS5o0MFmFeDrwlQr1TQn0bdNGXScCTQ/Z5PxFEkYGriCyjlqPNHdoCZ/Z5EOiLUHdhObEIbR7RFmyyi/W3ItWtwANd9ue7DL4rWQ4cGtGfA7voyxhM97e9r60Cth/FxkJ9KNczTzVlkytpCQGvAE5U1bV5yneR47EAk/EwJgfXOFS1HzgNy/bBrpzfVNVYgshdRVXvAn6KXdzA+tBxqjqw/lK15kLshMxYS+stkyZyEuZDxgpMcLhxqMkKnoD5ANbnrlPVe2K2U+QPo5MYHDxrt2NV9Te0lOvXAsdquKw0lPOxbADsyviZhLbEYC6tE3QJpj7fSFR1DfYcOgv+l6nqQ+ksqoaaVuYPwuo66pkcFeEy4InwfS0deOsnd/AM2Wc2LUGdd+zxYbmchmadGW3ZJ9jJ2cisMyNkn78Pq03OOjMupBU8m5x1ZpwUlgPUMDkqQlv2CXB77KwTir+q9LWwrO2ODdknwMKGZ50Z54flOUmtiMe3wrKxWWdGyD4fCN8fSmtNdbSl1P6nGidHRbgsLC/oROVFg+ffABqyYx9NbUAMQvYJ+Sa7agKLAHog68xYmtqADvBkagNioK0pVRZ3ov5eFFRwHMfpOB48HcdxSuDB03EcpwQePB3HcUrgwdNxHKcEtQyeIjJDRJ6f2o5YiMi2IjIxtR0xEGOHpipLDUVExonI9qntiIWITBeRF6a2IxYispWIbJDajuGoVfAMQfNs4BFMfafRiMgsEbkKuJfWy+6NJATNfYA7sBfdG62qH4LmEVhfu1tENk5tUxVC0Dwde53wltG2rzsisrWIXAzcR+td51pRi3lFRGQGNtTt/VhAHxc+jUREZgGfB/bG5qIRarKvixIyzL0xXdOXYCIpT9HSnWwUQU/yUOB0TPZvKqa+VatEIi8iMh04DjgG82EizT53tsZEoP8JO2f6qOm502mj3iIiZUf57FihLMC0CmWHkr3Q/VSFu9WPi8jHK9gQc7TUEoAK+3dD4Moa3bkvAzaq4M9E4PEa+QOwVQV/Nqt47kypUHZ9XFFh/x5U0Z+OXEw6GTwPxjQ289iwF7AnlqFljiqmb1iG5ZhkWBRUtV9EjiHfQdgAE/ndCbtqZhnNEuC2kibciylExeIoTLJrxWgbApsDB2JixmOxY/RMsOmhEm0LJokXk12AQ8jnzw7A/liAaJ+h8nosAy3KGGBeiXIjMYfBCvnrQzDtyr2xftbuzw9Ltr2GiOdO4APYBXc0JmGPg14T1rP4tBL4Wcm2/0r5fTEyBXXyZtMh3URsBsR52AnwDHBnt/T/OuTPVsBFtARzz0ptU0V/dsWepa3ERI07rqbfQV/6MLHph7AL7VoqKPSn/mAXguOxDLwfWJzapor+bI6JZ/eHY3NpF9osrOcpoWAuRGQ28HNV7dj9TXiG83HgKVVdMMrmtUdEtsIUd76jqjektqcqIrIr8FHgBG24GIaI9AFvxzK9I1T16cQmVUJEpgAfAqar6qdS21MVEdkcU3paqKrf63BbFkELxLbaBU/HcZxuUyZ4NvIfRsdxnNR48HQcxymBB0/HcZwSePB0HMcpgQdPx3GcEnjwdBzHKTGCr2jw3MUGNTiO4/QUhQNb0eD5dNwh1o7jOLWg8ISEftvuOI5TAg+ejuM4JfDg6TiOUwIPno7jOCXw4Ok4jlOCWgZPEdlQRD4pIlumtiUGIrKNiBwnInkEYWuPiLxJRN4jIrWcHqEoIvJ+EXljajtiICITReREEdkutS0xEJGZIjJXRDZJbctQahk8gVcApwL3iMgFPRBE347NabRIRE7pgSD6KeDrwMM9EkS/BFwjIr8VkT0T21KVLYHPAHeIyBU9EETfhImk/0VE5tcpiBbV8zwaODuP5p2IbAqUndJ1Z2wSqA0wJel12JQc/xdYlKO8Ao+o6tqS7T8HEdkG2Kxk8UOB92JTNqzG3in7KvCf5Js6Yo2qPlKy7WERkZ2xyc/KcC4wK3xfhflwBja9xroc5Zer6pKSbT+HMEndHpR/CfkGWlNY9AP3YxPE/Spn+SWqurxk289BRKZh07iU4cXAedjkduuw8+dGYD7wYM46FqnqMyXbfw4h+XlRyeJzMIHnydgMBgp8G+uDS3OUXwf8VUcJdCKyDugrpFVcUKr+aHJK1QcnU36uiCjRP64G/mwf0Z/jUvsTeQqFBxP7syayP6n72n/3mD/vyWHjuqL9stO37SeqqhT9AG/A5mMBy2yWA58ENshZHiDm7UpW58SS/pyIZQBgWdpfgX8FxhbwZ9OI/gwAlPEl2HNTW10rsQztDTnLdmJqlckV/VkZ6lmL3RlcAmxToGwnHltcVNKX7bHzheDLSuBzwMYF+to2HfBnp5L+HE7r+KwAHsMmlJuQ2J96zoccmIp1gs8D56hqntvbOtOHBc0TsAmt8tze1hXFrtS/BY5T1ZsT21MVxW4J/ws4SVXvT2xPVaZgAedsYIGqPpnYnqpMxILmScCFGvGRQhXqGjx/BfwL8KMeCJoA38Cm6v1hw4Nmxr8BG6nqLakNicQBwEM9EDQB7gPeAfysB4ImwBXAU9i5U4ugmVHL4Bl20mWp7YiF2p8jP0htRyxU9Q+pbYiJqv40tQ2xCBfnK1LbEQtVfYqa+lPXV5Ucx3FqjQdPx3GcEnjwdBzHKYEHT8dxnBJ48HQcxymBB898jB+ybDobpTYgMuNSG9ABpqU2IAYiksWYDZIa0gE8eObj8CHLpjMHnh3f3miCKMm08H1KYnMq0yaCs2tSQ+JxYFi+N6URIyEiEwixUEQm5i3nwXMUwo49JayeHNYbi4hsT0uw5cyUtkTiEFrDZz+S0pBIzAvLjURk96SWVCRknfPD6rtEZHpCc0bifW3f35+3kAfP0TkCyALmhLDeZD6HqTsB7Nbk7DNknV+gFTznNjn7DFnnO8LqGFqBp6kcCGQScn2YIE2tCMnQvLafTsubfXrwHIG2HZudkFOAeU3NPkPW+WZax30Szc4+D2GwrN5Ymp19zqN1YQPYqanZZ1vWmR2fCcAxNcw+34eNnc8YT97sU/NJSr0GOAv4JSaicFb4PG+EMgqckqf+2J/Q9l0R6jkMU9pZE+pcE9YP74IPm7ftZ8WGqJ0F7Fihzu9ggh7rQp1rwzKa3N0Ibe8e7L+nrQ8tADasUOdf2nzJxEqWYWpVnfbnXW3HJvNnLkEjt0R9M2jpbyqmfLUOuLELvggmcD3Un3+tUOdbg/3ZubM2fD7VBX+mDzl3fh++7zaM338fpg8tyXMc8xozr+2gZp8BYI/1bL9l2ObWTu+o9bQfK3hugakgnRzqPDmsb9EFH9425KBmHfCYCnW+Ntj/11DfCcCxwLgu+PPFYfwZAF5doc4Dgw/a5s+7u9THbh/ii2LyiZNL1idY1pz580z4vlcXfJmKiUAP9efXFercONif+fO98P2lXfBnl9C3hp47Zwyz7TuH6UMH52onpzEbYVp67cbcvr7oDFxIK1N7QTc685D2owTPtvqmhDo36KIPfZhCTvs+X1r25BxS91UUFH6N0OaMEFzaA+fNkeouLGQboc03DTknVgGnRapbgQe67M/pQ47PiliBO9R3YJf9uW2YC9umo9hYqA/leuapJm11NnY1BNMKPFZDq+0MeegtwKfztOEMRlUHsAfsmSRfP/BZVV2VzqryqOpjwNewCyqYP7X7A6EA1wMPtK0PYNl1U5lPEMkO/AmbvqOpHIcFTLA+d56qPh6zgSJ/GC2gtXNH2rHtD73HAIeJyAvKmfe/niuxZzJgtx1fSWhLDE7HrvAAt6vqL1MaU4WQOBxL65w4S1WXjVCk1qjqUixBAvNp2OSoKahpzd6ZrWKi6lHJHTxD9nl5WF1f1vki7EF6+4iP8XQh+xTj7SJySPhpcxE5RER27HTbnaIt+wT496ZmnRkh+7whrDY568y4HlOgh2ZnnRnZq1H9NDvrzMj62I9iZ51Q/FWlS8JyfTt2KTaVwQ/D+tPAtdi/9J1mJhbczw/rG2CzUza9U18Zlt9OakU8LgVoctaZERKIP4fvjc06M0L2CXBPk7PODG3NdPC9TtRfVEl+BTzbaZ6DmurzPwOIiGJTAcypZGFOVPUREbkO2Cf8NBZ7NvuFbrTfKVR1wGbWZXVqWyLRC9OqtNMrx6Wd/tQGRKYjfa7XXpI/gdZtFFhW8LNEtjiO08P0VPBU1d/RmhZ3LQ1/6O04Tn3pqeAZOCEsl+JZp+M4HaLngmfIPgEWetbpOE6n6Lng2Ub0VxMcx3Eyejl41gIRGScih4rINqltiYGITBeRD4rItNS2xEBEthaR94nImNG3rj8i8noR2S+1HTEQkTEi8k4ReUVqW4bDg2eHCEHzfcAjwDeAoxKbVIkQNE/H/DkHmJ3YpEqEoHkx8AfgqzR8apIQNG/FnvOfP9r2dSYEzYOx4a8XYYpPtaPoe57/2xmfQ8tzLKYz+RlMUCTTAh0vIrWaxyWnLuk04OOY4o9gGqBPAZNK+pN7moOi5PRnK+BUYH/sWI3F3tXcQESeHqHc+uhYApLTn12xoYc7AZMzm+rW14CxOfzpA96OCXZvREsLdFwN/fHgmZO1YbmkQh0fDp+ylDmx18ftwH6Uf8F7Q1qjzerA/cC2lPdnIoNFPooyMPomhVAsyJf15/nYBa4sKyuUXR9VRvkcQDV/Yp47z+LBMweq+rSIjKc13cNItGeeU2llnudgYrllWKuqMTvAadjVPQ/TgI8Cx2CZQZZ5fgC4umT7z4y+SSG2J39ffjGWef4TgzPP7Sl/cYytOTCOwYryI/FabBTdK7FjI8CjQOln7KoaO3iOI1+GLphG6xcwPdAs8/wBpplRhtjnTouCGnmzyal5h109r+2mht+Qts9N0XabDeMwJfpHscz1rJT2RPBnOqaYtQKT+JqT2qaK/myFPU/rx/RAN0ltU0V/dgVuwbKsxantqehLHyZr+VDoa5d2oc3Cep4SCuZCRGYDP1fVUTOwMLb9R6r61twNRCK0fZ6qJv+TRkTGYWrVt6nq/antqUr4l/2dwPe0B8QwROTFwJ7At1V1XVprqiMiuwIbq+o1qW2pSpgH6UDgXlX9fYfbsgiaI7Zl+G17h1HVNViG0xOEgPnV1HbEQlUfAi5IbEY0VPW21DbEQk2S8fJRN0yEv6rkOI5TAg+ejuM4JfDg6TiOUwIPno7jOCXw4Ok4jlOCov+2Fx0ilftv/w4wI2HbjuM0izUUjFdFg+d2Bbe/afRNOsbLE7btOE6zGEvB4Fn0tn1Nwe1jj/ktQuNfeHYcp2sUFk73Z56O4zgl8ODpOI5TAg+ejuM4JfDg6TiOUwIPno7jOCXw4NkFRGS2iPxYRHZJbUsMRORYEblARLZMbUsMROQiETlFRDZMbUtVRGQjEfmJiLxbRBqvmiYiO4vI9SKyZ2pbhuLBszu8DtgHuFFEFvZAEN0PU/a+p0eC6AHAicCiHgiiMzDR8q8AD4vIexoeRF+J+XO1iPy2TkG0qBjy0cDZBcSQT1TVMyrYV4rQ9h9UdYdI9QlwIfCqklVsjc2Tk+23tdi88guBpTnK96nqkSXbfg4iMgW4BtikZBWz2r4PYO/I3Qf8gnxTbCxS1c+WbPs5iMgs4GLyT10xlHZ/1mL+3AH8d46yY4CFqvqdkm0/B5G+94IeT7n3pKcC7RezNeFzC/nmaZoOnBRTuFtEvgzsVbL4FthUNtkFYC12zvycfNOmTAKOUNW1I20kIuuw8yz/i/IFpeqPptg0HCckkvFX4K6I9Y0Pdab87BLRn3mp/Yl8vJf0mD+p+9rDPebPqHEIG1RT6Dj6bXsxJqqqFP1gt4TrsKvmamzmyW0KlAfIMw1tXpYBlPEl2JMNu12NzbT4OWzqhzxlF0T0I2NdRX+yCc9WAI9hk9tNKFg2NheV9GX74AfBtruANxOyqpx9bVoH/NmppD+HYxPsDWDzTV0N7Fjw3JmyHpsq0eRnIU1iNfZ8+XtEviVKxErMpy8CC1T1ycT2VCW7CJwEXKiqsWf37CbPYMHiLuBY4Ceqllo1lNXYfPRXA3NV9a7E9jyLB8/ucA5wmaouTm1IJA4CxvdA0Mx4ObC04UETAFV9UEReiM2g2eSgmXEJcJOqLkptyFA8eHYBtYmseiVwojavd6duV7uOqj6W2oaY1DHQlCVcAGrpT0eeeYbpdgGi/NvtOI5TNzr1h9GhYbmfiJR9fcRxHKe2RA+eIes8PaxOAt4Ruw3HcZzUdCLzPBT7dwzs/cgzPPt0HKfXiBo827LOqW0/b4Rnn47j9BixM8+tgU0ZPERvErBv5Haeg4g8X0SWikh/+GmWiPSLyCURqn9eWJYdzlg3toRnh532AhNTG9ABNk9tQAxEJLsL3SKpIaNT+FyIGjxV9V5sbGw2c+XPsMBzeMx21sMTwNMMPpEU+H2Euo8Py+Mi1JWUEDAPDKtvS2lLDERkGmFWVxGZmdicyojI68LX1/XI466jwvIjSa0YARHZhBA8RWTTvOWiP/NU1adUdVlYfVpVl3Xj5WNVXQN8itbQNLBhe1+uUm/YsZkoxweK7Nya8kbs7gBgQQ9knx9r+35qKiMiMj8sx2ODERpLyDo/HVb3EpGXpLRnBOa2ff9k3kK9Nrb9QmwcLNg48i+o6vKKdc6lldL3MXhHN4oQKBdgJybYrWFjs8+QdX6C1vF5V5Ozz5B1ZspdY2n+n61H0RqIMxaIpqQVi5AcHdX205F5E6SeCp5t2SfYLXuMrPODtB4FTAA+2ODs843AS9vWp9Ls7PNjDJah6wNOSWRLDM7E/iPI2JiGZp9tWWcmyjEGOKCG2eeJDH7e2Ufe7FPzSUq9B7gTe66o4fvvgO1HKKPAtXnqj/kBxoW2fxmhrg9gWohLQ51Lw/oHu+DHq8N+vjO0/Zfw/YAKdX4PE1p4us0fBXbogj9HBvtXtvWhO4CtKtT5SFt9CiwPn7Fd8Oe8tmOT+fNTbMx/mfqeH+pZFpbrMBWhm7rgywTgxmH8+Y8Kde5PS3szOzbPAKd1wZ+XhPh0Z1vbd2K6nu3bCfD/sEd9md8rQpyT0drJO7Z9ApaxZLJoO4aDWztUdU1IpP4nQnUXAX/EsoHrgIOxkzWPSG5V+oCX0brteRHW+aroERwD/DumhvR/MEX4dcAfKtSZl0nAdrQeGexIS3i4LHthQSeTyHsr8KSOInwbiZkMHn68I/Ao5QSMUdW/i8hrsP10ExZ0DsAump1mABNHmdH22yygivrXtdidDpg/C7A/kGP8gZuHWbTuSqYyOH4BNm5eRHbH3qLJ+tBbgCUaouuI5Izk44G/04rOa4Hvj1ImSebZ1va5EeubEurcoMt+3EhLqV2BBzFdxqr1XkVkAd8cbU4CnmzzZQ3wrUh1FxayjdDmDtjz9fas99BIdSvwQJf9OSL4kPmzihHuLEv4c2CX/flu6GOZP/8P02gdycZCfSjXM0+1f8tPwm4jwDKg3P9KOaU5ntY+XwUcp6bQ1DhUtR84DXtkAHYBPjmdRdVQ05X8Ka3MeQUQbSqOBLT/2arAdap6T0J7qnIS1sfAHlWdrKpPj7B9YYr8YXQhrY7/44bv2Eagqr8Bfh1WnwC+n9CcGJxPq0NfpqoPJbQlBu1vXpyo3Xlc0BH0uX+2NvatEgA1wfEfhNU1wNdjt5E7eIbsMzPAs87ukb2g//mmZp0ZIfu8NKw2NuvMCNln9k5zk7POjAvD8skeSY5OCsuLY2edUPxVpR8B9MiObQQh+wT4cVJD4nEdQA9knRn3AzQ568wI2SfAn5IaEgltTXfTkXOnaPCs5T/s/0vohSkVoHXb3iv04jnRaz51xJ+eeknecRynW3jwdBzHKYEHT8dxnBJ48HQcxymBB0/HcZwS1C54ish0ETldRE5IbUsMRGRrEblQRPZJbUsMROT1InKZiGyV2paqiMgYEXmniFwsIhNGL1FvRGSKiJwoIp9PbUsMRGSmiJwnIgentmU4ahM8s6CJKeUcBxyS2KRKhKB5MSa68W7gzYlNqkQImrcCPwHmYEIYjSQEzYOBB7CBH+8gqNE3kSxoAosxQehDRy5Rb0LQPB+4Dxtzf+AoRZJQRaEnD1vnyCDHAv8AvB4L5plNO4rItZ00rgQfFZHRVPE3wOZs2gHzJ7tAvVdEXtZJ4wqQTTmQJ7vfHFP32QyT+wPTNvi8iHywRNvPL1EmFzn9mYWpSU2mpfAEcKmIrC7R7HYlyuRhWg5/BJiNqRf10To+m1U8dzqRVB0gItuOss0kTNXoVZhvmSrS22oYC/KpKrUpj8wmp/IIJmihaT6iwLERFVrGpPPl2c/GEf15Z2JfBiIr6Hw/qT/Sd39kf5ak80WUiIpk2qZYlPDzmrw2FvKr4E6YXbSBnPVuCnwJU3VZHRy5M3Y73fpgOpw/CP5kmpVnpbarpC+C3RncQUs0dhkwJ7VtJf0ZBxyGaW9mEmz9wCapbSvpz3RgXjg2mUTe4tR2VfBnK0xHt5+WpNylXWi3cPCsxTNPVX1cVT+GTYl7LnZb2PFJ4zqFqv5RVQ/A1OCvxg7MmhEL1RQ1fgrsjKmD/w7YkIYOs1TVNar6TeCFwNGYTu0ESooYp0ZVl6rqydjUvguwoNPkc+dBVX0XJs78Pey41PLckRB1820sMhv4uap2dM6bMEeQqOpjnWynW4jINsAjqlrmmVqtCPMdvRy4W4t0npoiIuOArdWmzW48IjIdmKqqj6S2JQYi8mJM2X3FaNtWbMfSzwKxrZbB03Ecp5uUCZ61uG13HMdpGh48HcdxSlDL4BleYv7HXhj1ASAiG4rIXg2eH30QIrKViLwitR2xEJHXisjmqe2IgRj7iMiU0beuPyIySUTeJCK1i1W1MyiwGzZ16SIROaoHguhRwA3A3SKybw8E0W8Ct4vIdSLyytTGROAG4AEROVdENkttTEVeiqn1LxaR43ogiB6MzWDxgIi8o05BtDaGDMNy4HnAfHojiK4DtgcuoflBVLC+szdwaw8EUcFeVzqc3giiK7DXyU7FgujxDQ+i/dhrjN+gRkG01L/tHbPGcRwnHU+p6rS8GxcNnoK9KP2iEoYVYTvgvdj44wHshezfYFPXPhyxHcHG0Xb6he/9sPHHY7AXmBW4ChtJsWyEckWZjI2pjz5T4BCOB2aG76uxl5gvwHyK+ULz5thY+E6/wH522/d+bJrnrwC/iNzONkCnM8DNgY9imbRi/e0P2OCTP0duayeCTkIHmQ28FZiInacDwPXYo6MlEdsZi41qfDx3idTDsdYzVGoP7KCvAq4EXp7apor+nBAO/CpsGOqmqW2q6M/CcHyewEbpTExtU0V/VgArseCyPyGpaOIHSzyyc+enwKtT21TRn8OwC3I/poC1RWqbsk+hzLNbiMjGwFzgm6p6d2p7qiIiszCJvS9pkStbTRGRA7Hhjedrb4yaOh64F7hK63hCFCA82zwZGw/+36ntqUoYYXQUcI7WbNRULYOn4zhO3Un+j5XjOE4T8eDpOI5TAg+ejuM4JfDg6TiOUwIPno7jOCXw4Ok4jlMCD56O4zgl8ODpOI5Tgv8PUzu6VjhR1SoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "plt.axis(\"off\")\n",
    "ax.set_aspect('equal')\n",
    "\n",
    "def draw_invP_2(x):\n",
    "    a = 0.2; oma = 1.0 - a; l = 1.0 - 2*a; \n",
    "    tma = 2.0 - a; tl = 1.0 - a;\n",
    "    draw_U([x[0]+oma,x[1]+a],[-l,0.0])\n",
    "    draw_U([x[0]+a,x[1]+tma],[l,0.0])\n",
    "    draw_U([x[0]+oma,x[1]+tma-tl],[0.0,-tl])\n",
    "    draw_U([x[0]+oma,x[1]+tma],[0.0,-tl])\n",
    "    draw_U([x[0]+a,x[1]+a],[0.0,tl])\n",
    "    draw_U([x[0]+a,x[1]+tl+a],[0.0,tl])\n",
    "\n",
    "nx = 6\n",
    "ny = 4\n",
    "plt.xlim(-0.1,nx + 0.1)\n",
    "plt.ylim(-0.1,ny + 0.1)\n",
    "\n",
    "for x in range(nx):\n",
    "    draw_U([x,0],[1,0])\n",
    "    draw_U([x+1,ny],[-1,0])\n",
    "for y in range(ny):\n",
    "    draw_U([0,y+1],[0,-1])\n",
    "    draw_U([nx,y],[0,1])\n",
    "    \n",
    "for x in range(nx):\n",
    "    for y in range(ny):\n",
    "        if x != 0 or y >= 2:\n",
    "            draw_invP([x,y])\n",
    "draw_invP_2([0,0])\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ccf89e4",
   "metadata": {},
   "source": [
    "Combining the factorial and $1/N_c$ powers, we find\n",
    "for the lowest-order $\\beta$ contribution for a $n_s \\times n_t$ Wilson loop\n",
    "$$\n",
    "\\begin{align}\n",
    "\\langle W_C \\rangle\n",
    "&= \n",
    "\\left(\\frac{\\beta}{2N_c^2}\\right)^{n_s n_t} \\textrm{Tr}\\:[\\mathbb{1}]( 1 + O(\\beta)) \\\\\n",
    "&= \n",
    "N_c\\exp\\left(n_s n_t\\log\\left(\\frac{\\beta}{2N_c^2}\\right)\\right)( 1 + O(\\beta)) \\,,\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c829ed8b",
   "metadata": {},
   "source": [
    "If we now again identify the Wilson loop with the static quark potential $V(n_s)$, we find\n",
    "$$\n",
    "\\begin{align}\n",
    "\\langle W_C \\rangle \\propto \\exp(-n_t V(n_s))\n",
    "\\end{align}\n",
    "$$\n",
    "for sufficiently large $n_t$.  We can therefore at leading order in $\\beta$\n",
    "identify\n",
    "$$\n",
    "\\begin{align}\n",
    "V(n_s) = -n_s\\log\\left(\\frac{\\beta}{2N_c^2}\\right) \\,.\n",
    "\\end{align}\n",
    "$$\n",
    "We indeed find a linear growth of the potential in distance $n_s$.  The strong-beta expansion\n",
    "converges quite poorly such that the numerical value of the resulting string tension is not\n",
    "significant, however, we can read off\n",
    "$$\n",
    "\\begin{align}\n",
    "\\sigma = -\\log\\left(\\frac{\\beta}{2N_c^2}\\right)\n",
    "\\end{align}\n",
    "$$\n",
    "which for our coarsest simulation at $\\beta=5.5$ corresponding to $a=0.25$ fm would yield\n",
    "$$\n",
    "\\begin{align}\n",
    "\\sigma = -\\log\\left(\\frac{\\beta}{2N_c^2}\\right) \\approx 1.2 \\times \\frac{0.1973 \\:\\textrm{GeV}\\:\\textrm{fm}}{(0.25 \\:\\textrm{fm})^2} \\approx 3.8 \\frac{\\textrm{GeV}}{\\textrm{fm}}\n",
    "\\end{align}\n",
    "$$\n",
    "quite a bit larger than the approximate $1$ GeV / fm found numerically in the continuum limit of pure QCD in the last chapter."
   ]
  }
 ],
 "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
}
