{ "cells": [ { "cell_type": "markdown", "id": "ccd7827a", "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 14: Hadron Spectroscopy\n", "\n", "In this section we discuss the spectrum of full QCD. To this end we first need to identify the good quantum numbers.\n", "Based on previous chapters, we need to consider space-time, flavor, and charge-conjugation symmetries. On the lattice the Hamitlonian can be classified using the irreducible representations of the full octahedral group, as discussed in chapter 11. In the continuum, however, the irreducible representations of $O(3)$ become relevant and are typically used to classify the states.\n", "\n", "## Quantum numbers\n", "\n", "We use the angular momentum Casimir $J$ of O(3), parity $P$, and charge conjugation $C$ to classify the different space-time and charge conjugation behavior. We amend this discussion by also using the Casimir $I$ of the **isospin** flavor SU(2) symmetry in the up-down-quark sector. One typically writes the states as $J^{PC}$.\n", "\n", "\n", "## $N_f=2$ meson spectrum\n", "\n", "We first discuss a theory with only up and down quarks with degenerate masses $m_u=m_d=m_l$ such that isospin becomes an exact symmetry.\n", "\n", "The up and down quark fields live in an SU(2) doublet with $I=1/2$. For mesonic states\n", "we can combine two such fields to $I=0$ and $I=1$ irreducible representations.\n", "\n", "For $I=0$, i.e., singlets in isospin, operators that create these states are of the form\n", "$$\n", "\\begin{align}\n", "O_{I=0} &= \\frac1{\\sqrt{2}}(\\bar{u} \\Gamma u + \\bar{d} \\Gamma d) \\,,\n", "\\end{align}\n", "$$\n", "where $\\Gamma$ can be a matrix in spin and space-time. A two-point correlation function\n", "between this operator and another one, say\n", "$$\n", "\\begin{align}\n", "O^\\prime_{I=0} &= \\frac1{\\sqrt{2}}(\\bar{u} \\Gamma^\\prime u + \\bar{d} \\Gamma^\\prime d)\n", "\\end{align}\n", "$$\n", "can then be used to measure the $I=0$ spectrum. The Wick contractions\n", "$$\n", "\\begin{align}\n", "C(t) &= \\langle O^\\prime_{I=0}(t) O_{I=0}(0) \\rangle \\\\\n", "&=-\\langle \\text{Tr}[ D^{-1}[U,m_l] \\Gamma D^{-1}[U,m_l] \\Gamma' ] \\rangle_U\n", "+2\\langle \\text{Tr}[ D^{-1}[U,m_l] \\Gamma] \\text{Tr}[ D^{-1}[U,m_l] \\Gamma' ] \\rangle_U \\,.\n", "\\end{align}\n", "$$\n", "The first Wick contraction is called **quark-connected** or sometimes only **connected** and the\n", "second is called **quark-disconnected** or just **disconnected** diagram. For $I=0$ we note that we\n", "have both types of contraction. This typically leads to $I=0$ correlation functions to be particularly\n", "noisy since the disconnected diagrams have constant noise for $t\\to\\infty$ but exponentially falling signal.\n", "For the connected diagrams both the noise and the signal fall exponentially for $t\\to \\infty$, albeit at different rates. The **Lepage** argument of chapter 3 explains this behavior.\n", "\n", "For $I=1$, i.e., triplets in isospin, operators are\n", "$$\n", "\\begin{align}\n", "O_{I=1,I_3=+1} &= \\bar{u} \\Gamma d \\,, \\\\\n", "O_{I=1,I_3=0} &= \\frac12(\\bar{u} \\Gamma u - \\bar{d} \\Gamma d) \\,, \\\\\n", "O_{I=1,I_3=-1} &= \\bar{d} \\Gamma u \\,.\n", "\\end{align}\n", "$$\n", "Contractions between two such operators only have connected diagrams.\n", "\n", "We next discuss the $J^{PC}$ quantum numbers and give the corresponding names and experimentally measured\n", "masses" ] }, { "cell_type": "markdown", "id": "e8967c42", "metadata": {}, "source": [ "| I, $I_3$ | $J^{PC}$ | Name | Mass |\n", "| --- | --- | --- | --- |\n", "| 1,0 | $0^{-+}$ | $\\pi^0$ | 134.9766(6) MeV |\n", "| 1,$\\pm$ 1 | $0^{-}$ | $\\pi^\\pm$ | 139.57018(35) MeV |\n", "| 0,0 | $0^{-+}$ | $\\eta$ | 547.853(24) MeV |\n", "| 1,$\\pm$ 1 | $1^{-}$ | $\\rho^-$ | 775.4(4) MeV (center), this state is a resonance |\n", "| 1,0 | $1^{-+}$ | $\\rho^0$ | 775.49(34) MeV (center), this state is a resonance |\n", "| 0,0 | $1^{--}$ | $\\omega$ | 782.65(12) MeV (center), this state is a narrow resonance |\n", "\n", "Some of these states are marked as resonances, i.e., they are not themselves asymptotic states but can decay to lighter particles. For physical masses, e.g., the $\\rho^0$ can decay in two pions $\\pi^+ \\pi^-$ with relative momentum such that the combined operator has the quantum numbers of the $\\rho$.\n", "\n", "The lowest-lying two-pion state that the $\\rho$ can decay into has approximately the energy of two non-interacting pions with smallest allowed momentum in the finite volume. When the $\\rho$ becomes a stable particle for the lattice Hamiltonian therefore depends both on the masses and on the volume." ] }, { "cell_type": "code", "execution_count": 1, "id": "ff6dff75", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "For volumes with mpi L = 4\n", "mpi = 140 MeV rho stable if mrho ~< 521.3868489532042 MeV\n", "mpi = 175 MeV rho stable if mrho ~< 651.7335611915053 MeV\n", "mpi = 200 MeV rho stable if mrho ~< 744.8383556474347 MeV\n", "mpi = 250 MeV rho stable if mrho ~< 931.0479445592932 MeV\n" ] } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "def energy_lowest_moving_pion(mpi, L):\n", " return (mpi**2. + (2.0*np.pi/L)**2.)**0.5\n", "\n", "# for mpiL = 4\n", "print(\"For volumes with mpi L = 4\")\n", "for mpi in [140, 175, 200, 250]:\n", " print(\"mpi =\",mpi, \"MeV\", \"rho stable if mrho ~<\", 2.*energy_lowest_moving_pion(mpi, 4 / mpi), \"MeV\")" ] }, { "cell_type": "markdown", "id": "b9a73c97", "metadata": {}, "source": [ "In order to study the $N_f=2$ meson spectrum from first principles lattice QCD, we generate a small ensemble using the HMC from chapter 13. We use the following parameters: Shamir DWF fermions ($b=1$, $c=0$)\n", "and the Iwasaki gauge action ($c_1=-0.331$) with\n", "$$\n", "\\begin{align}\n", "V &= 16^3 \\times 32 \\,, \\\\\n", "L_s &= 16 \\,, \\\\\n", "m_u &= m_d \\in \\{0.01, 0.02, 0.04\\} \\,,\\\\\n", "M_5 &= 1.8 \\,, \\\\\n", "\\beta &= 2.13 \\,, \\\\\n", "\\tau &= 0.6 \\,,\n", "\\end{align}\n", "$$\n", "where $\\tau$ is the MD trajectory length. Let us first study the convergence of the plaquette." ] }, { "cell_type": "code", "execution_count": 134, "id": "82fb291a", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.11694571822930958\n", "0.20272987289475408\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEWCAYAAABv+EDhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABK5klEQVR4nO29e5gcZZX4/zndczeEyzCQkMtMNIAJQUCiCChIcPcHmIUNRIFMhGjYkWTdrLpe0Kx+gRhQwAtCgkZkJcwo2cSoAV0gOIFdXEYJCwhJBgnJJAQmF8I1l8lMT5/fH1U9U9PT1V19m+6ZOZ/nqaer6/LWW9XV57zvOec9r6gqhmEYxvAjVOgKGIZhGIXBFIBhGMYwxRSAYRjGMMUUgGEYxjDFFIBhGMYwxRSAYRjGMMUUgDGsEJHHROSaQtfDMIoBUwDGkERE2kTkoIjsE5FdIvILERlR6HoBiMjHRWRH3LbrRaSxUHUyhiemAIyhzD+o6gjgg8BU4N8LXB/DKCpMARhDHlV9FfgvYIp3u4i8T0SaRWSviLwuIk0icoRn/2ki8n8i8q6IrBCR+0XkO+6+OSLyRFx5KiIT3fVyEblNRLa7PZCfiEiliLzHrctxbu9kn4jMAr4JXO5+f84t43AR+bmItIvIqyLyHREJ5/FRGcMMUwDGkEdExgEXAc/E7wJuBo4DJgHjgOvdc8qA3wL3AUcBK4HL0rjsd4ETgFOBicAY4Nuquh+4EHhNVUe4yy+Bm4AV7vdT3DJ+AUTc808D/h4w/4WRM0wBGEOZ34rIW8ATwOM4QrYHVd2sqmtV9ZCq7gF+AJzr7v4IUAr8SFW7VHUV8FSQi4qIAA3Al1T1DVV91732FUErLiLH4iitL6rqflXdDfwwnTIMIxUlha6AYeSRf1TVR70bHNncs34scDvwMeAwnAbRm+7u44BXtW+2xG0Br1sDVAFPe64nQDrmm1ocBdTuKSMEvJJGGYaRFOsBGMOZmwAFTlbVkcBsHEEN0A6MEa/GgPGe9f04Qh4AERnl2fc6cBA4SVWPcJfDXYc07jXjid/2CnAIONpTxkhVPSnNezQMX0wBGMOZw4B9wNsiMgb4qmffkzj29wUiUioilwIf9ux/DjhJRE4VkQpc3wGAqkaBnwE/FJFjAERkjIj8f+4hu4BqETncU94uoE5EQm4Z7cAjwPdFZKSIhFyn9bkYRo4wBWAMZ27ACRF9G/g9sDq2Q1U7gUuBOcAbwOVx+/8G3Ag8CryE42fw8nVgM9AiIu+4x53ontsK/ArYIiJvichxOE5mgL0i8n/u+lVAGbARxzS1Chidg/s2DADEJoQxjGCIyC+AHapq4wmMIYH1AAzDMIYppgAMwzCGKWYCMgzDGKZYD8AwDGOYMqgGgh199NFaV1dX6GoYhmEMKp5++unXVbUmfvugUgB1dXWsX7++0NUwDMMYVIhIwlHsZgIyDMMYppgCMAzDGKaYAjAMwximDCofgGEYQ5Ouri527NhBR0dHoasyqKmoqGDs2LGUlpYGOt4UgGEYBWfHjh0cdthh1NXV0TcBqxEUVWXv3r3s2LGDCRMmBDpn2JiA1m1dx8l3nczOfTv7rBuGUXg6Ojqorq424Z8FIkJ1dXVavahh0QNYt3Ud0381nUORTk78Rj37jmxBSjpZ9PgilnxySaGrZxgGmPDPAek+w2HRA1jw0AI6I510a4R3RrYQDR+gWyOs3Lgy9cmGYRhDlEAKQEQuEJEXRWSziFyXYP8cEdkjIs+6yzWefbeIyAYR2SQiP47NsCQil4vIX91938vdLfXnpW89QuTlc6CzCsoOOBu7Ktnzi6VUVubzyoZhGA5bt27ljDPOYOLEiVx++eV0dnYmPO7mm29m4sSJnHjiiTz88MMAvPLKK5x33nlMnjyZk046idtvvz0ndUqpAEQkDCwBLgQmA1eKyOQEh65Q1VPd5W733LOAs4EPAFOADwHnikg1cCtwvjvF3SgROT8nd5SAxodbCde19Ap/AOlmwrRmtm7N11UNw8gn7e1w7rmwc5C48r7+9a/zpS99ic2bN3PkkUfy85//vN8xGzdu5P7772fDhg089NBDzJ8/n+7ubkpKSvj+97/Pxo0baWlpYcmSJWzcuDHrOgXpAXwY2KyqW9xZku4HLglYvgIVOLMaleNMcr0LeC/wkqrucY97FLgsnYqnww1/WUBUXG3bVQmRMijppP2oVYwalfxcwzCKk0WL4Ikn4MYbsy+rra2N97///cyZM4cTTjiB+vp6Hn30Uc4++2yOP/54/vKXv2RVvqrS3NzMzJkzAbj66qv57W9/2++43/3ud1xxxRWUl5czYcIEJk6cyF/+8hdGjx7NBz/4QQAOO+wwJk2axKuvvppVnSCYAhiDM0F1jB3utnguc006q0RkHICqPgmsw5lgux14WFU34UyVd6KI1IlICfCPwLhEFxeRBhFZLyLr9+zZk+iQlKz9zFrq9jYQ6qjhzJ3Led87cwl11DB1y4qMyjMMo3BUVoII3HUXRKPOpwhZm3M3b97Mv/3bv9Ha2kprayu//OUveeKJJ7jtttu46aab+h3/4osvcuqppyZc3nrrrT7H7t27lyOOOIKSEifuZuzYsQkF+Kuvvsq4cb2iMNFxbW1tPPPMM5xxxhnZ3TC5iwJ6APiVqh4Skc8D9wLTRGQiMAkY6x63VkQ+pqr/IyLzgBVAFPhf4H2JClbVZcAygKlTp2Y0ecGoEaPYcscSHEsWwExgaSZFGYZRYLZsga98BX77WzhwAKqqYMYMuO227MqdMGECJ598MgAnnXQS559/PiLCySefTFtbW7/jTzzxRJ599tnsLpom+/bt47LLLuNHP/oRI0eOzLq8ID2AV+nbOh/rbutBVfeq6iH3693A6e76DKBFVfep6j7gv4Az3XMeUNUzVPVM4EXgb5nfRjBsLIBhDH5Gj4aRI6GjAyoqnM+RI8nanFteXt6zHgqFer6HQiEikUi/49PpAVRXV/PWW2/1lLNjxw7GjOlvSBkzZgyvvNJrcPEe19XVxWWXXUZ9fT2XXnppdjcbu88AxzwFHC8iE0SkDLgCWOM9QERGe75eDGxy17fjOH1LRKQUODe2T0SOcT+PBObjKI68ERsL0Lqnlfpf1/esL3p8UT4vaxhGHti1C669FlpanM9COIJjPYBEyxFHHNHnWBHhvPPOY9WqVQDce++9XHJJf1fqxRdfzP3338+hQ4fYunUrL730Eh/+8IdRVebOncukSZP48pe/nLN7SKkAVDUCfAF4GEd4/6eqbhCRG0XkYvewBW4453PAAmCOu30V8DLwPPAc8JyqPuDuu11ENgJ/Ar6rqnntAcTGAkQ0QsuOFg50HSBiYwEMY1CyejUsWQKnnOJ8rl5d6Bql5nvf+x4/+MEPmDhxInv37mXu3LkArFmzhm9/+9uAY3r69Kc/zeTJk7ngggtYsmQJ4XCYP/3pT9x33300Nzf39DL+8Ic/ZF2nQTUn8NSpUzXdCWHWbV3HgocW0Dijkc/+7rM8u/NZFOeeK0sqWT5jOTMnz8xHdQ3DCMimTZuYNGlSoasxJEj0LEXkaVWdGn/skB4JHDP7bNrTytm3fpZndj7TI/wBurWb5q3NBayhYRhG4RjSCsCbAmL/iGf77CsLl9HZ3cmqjasKUznDMIwCM6QVQJ8UECG35a8Q/vOXmXvaXGqqalgx08YCGIYxPBnSCiBRCoiQljH7swdZ+sml7P7qbs6bcF4Ba2gYhlE4hrQCSJQCIhrq5A9tZvYxDMMY0goglgKiIlrDrWct56TOuZRFzOxjGIYBQ1wBxFJAHLxhN1+5aCYvfG8phxaZ2ccwjIEnm3TQMbq7uznttNOYPn16Tuo0pBWAYRhDlKYmqKuDUMj5bGoqdI1Skk066Bi33357TsdLmAIwDGNw0dQEDQ2wbRuoOp8NDVkpgWJPBw1OXqDf//73XHPNNf3OyxRTAIZhDC4WLnTSgHo5cMDZngXFng76i1/8IrfccguhUO7E9rCYFN4wjCHE9u3pbQ9IMaeDfvDBBznmmGM4/fTTeeyxx3JWrikAwzAGF+PHO2afRNuzIJN00JdffnnCsh577LE+GUG96aBLSkrSTge9Zs0a1qxZwx/+8Ac6Ojp45513mD17No2NjZnernNvWZ1tGIYx0Cxe7MwC46Wqytk+gAxkOuibb76ZHTt20NbWxv3338+0adOyFv5gCsAwjMFGfT0sWwa1tc5ckLW1zvf6+kLXLCnZpIPOF0M+HTRNTY5zaPt2OOooZ9sbbzjdxcWLi/6lMYzhgKWDzh3ppIMe2j6AWLhYLGJg797efbHQMTAlYBjGsGRom4AShYt5yUHomGEYxmBlaCuAIGFhWYaOGYZhDFaGtgIIEhaWZeiYYRjGYGVoK4BE4WJeChA6ZhiGUSwMbQUQHy5WXe0sIkTG1HLjuGXsPN8cwIZhDE+GtgIARwm0tUE0Cq+/7izRKAsubuOGl+q58cZCV9AwjOFAtumgP/e5z3HMMccwZcqUnNVp6CuAOCornc7AXXc5OuGuu5zvlZWFrplhGEFZt3UdJ991Mjv37eyzXsxkmw56zpw5PPTQQzmt07BTAFu2wKxZva6Bqiqnk7B1a2HrZRhGMNZtXcf0X02ndU8r9b+u71lf9PiijMscDOmgzznnHI6KDWbNEcNOAYweDSNHQkcHVFQ4nyNHwqhRha6ZYRhBWPDQAjojnUQ0QsuOFg50HSCiEVZuXJlVucWeDjofDO2RwD68sG8dR3xzAf85fS13/GoTjRUL+Pa+tYwaYVrAMIqdR2Y/wuzVsx3hH3EGelaWVLL0k0uzKreY00Hni2GnANZtXcf/TZpOZ6STm1rraTmmhc5IJ4seX8SSTy4pdPUMw0hB6+uttLzaK/wBurWb5q3NzJw8M+NyizkddL4YdiagfHUfDcMYGGL/YXBa/mXhMjq7O1m1cdWA1mMg00Hni2GnAB6Z/QjnVJxAVZf0dh+lLOvuo2EYA8Paz6yl4fQGaqpqWD5jOXNPm0tNVQ0rZq4odNWSkm066CuvvJIzzzyTF198kbFjxyaMIkqXQOmgReQC4HYgDNytqt+N2z8HuBWIeSvuVNW73X23AJ/EUTZrgX9VVRWRK4FvAgq8BsxW1deT1SOjdNBxrPvZQqZvu4kDpb3byiIw96jzWfqlR7Mq2zCMzLB00LkjnXTQKXsAIhIGlgAXApOBK0VkcoJDV6jqqe4SE/5nAWcDHwCmAB8CzhWREhyFcp6qfgD4K/CFNO4xYxZsuJVO964ruxzh31kCq3atG4jLG4ZhFA1BTEAfBjar6hZV7QTuB/obrxKjQAVQBpQDpcAuQNzlPSIiwEicXkDeWfvzLhqehpr9sPw3MPcZZ/3+/4wOxOUNwzCKhiAKYAzwiuf7DndbPJeJyF9FZJWIjANQ1SeBdUC7uzysqptUtQuYBzyPI/gnAwkNWiLSICLrRWT9nj17gt5Xf5qaoK6OUftgyR9g960wcyMs/b2zfsqe2szLNgwjawbT7ITFSrrPMFdO4AeAOtecsxa4F0BEJgKTgLE4SmOaiHxMREpxFMBpwHE4JqBvJCpYVZep6lRVnVpTU5NZ7WIzg23blnD3fqr4l32LLSWEYRSIiooK9u7da0ogC1SVvXv3UlFREficIOMAXgXGeb6PpdfZG7uwZ65F7gZucddnAC2qug9ARP4LOBPocM972d3+n8B1gWudLj4zgynwitRynS7md1X11M+A227LWy0Mw/Bh7Nix7Nixg6x6+QYVFRWMHTs28PFBFMBTwPEiMgFH8F8BzPIeICKjVbXd/XoxsMld3w78k4jcjGPzPxf4kVvOZBGpUdU9wN95zsk9vrN+CROkjbJy6LSUEIZRMEpLS5kwYUKhqzHsSKkAVDUiIl8AHsYJA71HVTeIyI3AelVdAywQkYuBCPAGMMc9fRUwDcfWr8BDqvoAgIjcAPy3iHQB2zzn5J7x4xOaf/ZUjufazzrWoWXLoL09wbmGYRhDlEDjAIqFjMcBxHwAXjNQVZUj9ettQhjDMIY2GY8DGBLEzwxWW2vC3zCMYc/wSQZXX28C3zAMw8Pw6AEYhmEY/TAFYBiGMUwxBeDDYJxz1DAMIx1MASQgH3OOGoZhFBumABLgnTTmsZdt0hjDMIYmpgA8rNu6jgm3T6Cjq4P3H/1+UCEazt2co4ZhGMXE8AkDTcG6reu4sOlCDnUfcsYsg5O8wuVgRzdXfKOZyO8yn3PUMAyjmBjWPYB1P1vIyV8sY+dhwueW/J0j/KF3toIYCpR0csRZAzvnqGEYRj4ZtgogNjVk68gu6i+FXeXdvS1/L90l8PLfUxGtYeWni3vOUcMwjHQYtiagBRtupXMkRMLwp3FwqDTxcaFwlEnHvY8Tnn+Y8yxZoWEYQ4hh2wN45J4uztkG5V0e4Z+gBxAlyu6jV7F69YBWzzAMI+8M2x5A6+RjaRm7y7flXxoqJSxhykvKWTHTTD+GYQw9hq0CWHBJKZ0HnfVQFKIex29lSSXLZyxn5mSL+DEMY+gybE1Aa//lKRqOPJ/DDwklUfpE/XRrN81bmwtWN8MwjIFg2CqAUSNGseRLjzJu3ElES52OUGVJJWXhMjq7O1m10UI+DcMY2gxbBRBj7WfW0nB6AzVVNSyfsZy5p82lpqrG7P6GYQx5hseUkGnQ3g5XXAErVtgE8YZhDA2G95SQabBoETzxBNx4Y6FrYhiGkV9MAbhUVjrTBd91F0SjzqeIs90wDGMoYgrAZcsWmDULqqqc71VVzhTCW7cWtl6GYRj5whSAy+jRMHIkdHRARYXzOXIkbDpoM4MZhjE0MQXgYdcuuPZaaGlxPp9/12YGMwxj6GJRQEk4+a6Tad3TSkQjVJVUcSDiTA5TU1XD7q/uHrB6GIZhZINFAWXAI7Mf4Zzac/oI/3zMDGYT0BuGUQhMAXhpaoK6OgiFoK6O1vvvpOXVlh7hD7lPE2ET0BuGUShMAcRoaoKGBti2DVRh2zYWPHcznV0dQP7SRHgnoG/ZYRPQG4YxcARSACJygYi8KCKbReS6BPvniMgeEXnWXa7x7LtFRDaIyCYR+bE4HOY59lkReV1EfpTD+0qbyNcXwoEDfbatvVdp2FSV1zQRA2VmMgzDiCelE1hEwsDfgL8DdgBPAVeq6kbPMXOAqar6hbhzzwJuBc5xNz0BfENVH4s77mngS6r638nqkk8nsEoISTQjjIgzMixPxExAB7p6lU9ZuIy5p801JWAYRk7Ixgn8YWCzqm5R1U7gfuCSgNdVoAIoA8qBUmBXXMVOAI4B/idgmTklNgJ4G+MTHzDeZ3uOiJmAwLKRGoYxsARRAGOAVzzfd7jb4rlMRP4qIqtEZByAqj4JrAPa3eVhVd0Ud94VwAotUDxqbATwDaWL2U9V351VVbB4cV6vb9lIDcMoFLlyAj8A1KnqB4C1wL0AIjIRmASMxVEa00TkY3HnXgH8yq9gEWkQkfUisn7Pnj05qm4vsRHAy7vrmV+yjG3UogjU1sKyZU4+iDwyasQolnxyCbu/upuZk2ey9JNL2f3V3Zw34by8XtcwDCOIAngVGOf5Ptbd1oOq7lXVQ+7Xu4HT3fUZQIuq7lPVfcB/AWfGzhORU4ASVX3a7+KqukxVp6rq1JqamgDVTZ/YCOAvr6/nlvltXDYjCm1teRf+hmEYhSSIAngKOF5EJohIGU6LfY33ABEZ7fl6MRAz82wHzhWREhEpBc717AO4kiSt/4Fi9WpYsgROOcX5XL06+fE2cMswjKFASgWgqhHgC8DDOML7P1V1g4jcKCIXu4ctcEM9nwMWAHPc7auAl4HngeeA51T1AU/xn6YIFEA65HLglikSwzAKieUCAmcQ2MKFsH27E/WzeLGv+SdX+YFiiqQz0sk5tefQ8moLnZFOGk5vYMknl+TktgzDMMByAfmTYAQwDQ3O9gTkauCWjQA2DKPQmAJY2H8EMAcOONtd2tvh3HNh505ofb01J/mBbASwYRiFxhTA9u2Jt2/b1pMYruyEOsb+dxM33pi7gVu5UiSGYRiZYgrAb6SvSI9ZqHrfNpbRwFt3NfHC19ai67MfuGUjgA3DKDTmBI75ALxmIBHHHxDHdqnlm7PauO02GDUqu8vu3LeTRY8vYuXGlSz95FKatzazauMqVsxcYYPADMPIKeYE9qO+3hnxW1vrCP7a2oTCH2CsbmfkyOyFP/iPAAYsNNQwjAHBFAA4SqCtzcn62dbmKIEEvDliPDvzKIttchjDMAYSUwAJePOri+kIVfXbXl2+j9WX9YaH5nogl4WGGoYxkJgCSMDCDfXMjS7j3fLqvjv27u0ZI5CP1rqFhhqGMZCYAvAQmxvgrrvgl9Sz99CInn3r6uDkebAzdIB1d/wbFzZdyKHIoZy21hOFhkaiEe577j7zBRiGkXNMAXjYsgXuPLOJbVJHNyFq2QY4wn/6LGg9Guovhenn7yLSHWHUiFGBW+sxE9GvN/6aUbeNovw75dz9f3cz4fYJTLpzEjv37ewTGhqjK9rFmr+tMV+AYRg5xxSAh9HNTVzzVAPjdRshZ1YAABZcCJ1hiIShZSwcKINuunn13VcTDuTyjhwGWPX0Oj7xH9PZuHsjn1r5KXbt30Vndyf/9MA/0fZWGy/ufZFFjy/qmRwmLGFC0venSdS7sGRyhmFkg40D8FJX5wz+iqN9BMy+tFf4AwiCiBDVKETKINwJAoeXH0543zje/PFarp09iqVLofrfT+aNUCuEI76X9iaTa3+3ndmrZzumJU/vYvmM5cycPBOwZHKGYQTHxgEEwSctROvRfYU/gEZKiO6cBAcPdzcIREp5e38Hb4Ra0Y8t4q67HJ/CG7c/AtvOgc7KhOWXhcv6mI6CpImwiCHDMLLFFIAXn7QQCy4ppbO8BOhN20BJF+XVu5F940CiEFLCoVIoOeS09E9aSWWl06koG9MKY1ug7GDC8iPRSELh7r1efJoIixgyDCNbTAF4WbzYmQjeS1UVa0/7YcKJ2//+nRVw3yPIK+dAZxXdIbfFHhVKHl9Mx+h1tM84mc5PzHdMRAAJLG5RjfYR7kEmik/VSzD/gGEYqTAfQDxJJodZt3UdCx5awNrPrGXTnk1c8tMFnNV5PY8fdRUd3Z5cQgo1lWN48+CbRLSTEzpmc+rkKn63ZQUj9p/MWyP/h6hGOX306byw+wXKS8r5zeW/SSsHkHdimsqSSrq1m87uzh5FYf4BwzBi+PkATAEExM/pOrJ8JO8ceoeIRhAETdDEr66o4fWvB58tLAjJkskteGhBTmYtMwxjaGBO4Czxc7qKSI+55qfTf8qoqjF9zTxdlZy+MzO7fHs7nHLJOib9uL8pxy+Z3HkTzjP/gGEYgTAFEBA/ofqT6T/pEcQTj5rI6/vfpGcAAUCoi0fKPo8ctpOyE9OzxV/73XX8dcp0XnwjeLqJdVvXcdY9Z/HkjidtshnDMJJiCiAgqZyulZUw7dYFRKKus7erErpLnIigyjc4dn49JVcFyxcUS0mxpmsBhDtRidD80p96eh1NzzclVCQxM1XbW210RDp6tockZJPNGIbRD1MAAUkVmrllC8zYt5bwcw2wv4bQ75bDoZGOOUhgb0ULByPBYvW3bIFZs6BylTt+IFIOpYcACBHiYORgQkXiraN3JHFZqCyjWcsMwxjamAIISKrQzNGj4dj3jEIfXELFHbuJvjCTE/74Ah+qmUZJtIqIGyIa6q7kprP72+K96SNGj4aRI6HjsFYY1+KMLXCJEqWzuzNhT8BrpurWbsBRVvddel+Pf8AwDCOGKYCAjBoxin//4BImPbibjx7V1+kaY9cuuPZaaGmB+fPhmEmtbHi7pUf4A0S1mx/9rr8tftEieOIJuPHG3rKOnLWAcFls/EDfnypRT8AmmjcMIx0sDDQN5s+Hn/4UPv95WBogoMYbq09XJUg3lHTC/hq4dTcVFc5xHR39z62ogK17nFDP+55rYt+hg2ios/+B9IZ3JhsbYOGfhjF8sTDQLPDOExCN0pPjJxTqzfiZKGSzq7uL+g/Uc1RFDR/ds5zw83Nhfw3lD6ygvh62bu2198cGIFdV0bMvFupZemAcqlHAafn3qZsnvDPICGLDMIwYpgACkEhIH3+8sx4z2cRCNv+2dyP1/zaB6XdP4+U9L/KerTvY+/XdTAnNRB9YSsUdu+n623k9k8v32Ps7nFZ/Rwc9+2KK543b18LTDXDwcKKRkj5185p4ko0NyBRLKWEYQxdTAAHwCmmAAwfgpZdAtbc3sKZrARLuIBqK0nJ0BwfKIBKClTv/CE1NffwD115Ln8nl/fb1KJ7oKPjDEmTfOCTs9AT8ksTlEpuk3jCGNuYDCMillzqKYMYMxxewZQt0dzut9GOPhZ372jntovfy3NiOnrTRlV2w/Dcwc38ttLVldN1582DZMigrg0OlO5k8bxG7j+6f/iEfET5en4KllDCMwUtWPgARuUBEXhSRzSJyXYL9c0Rkj4g86y7XePbdIiIbRGSTiPxYRMTdXiYiy0TkbyLSKiKXZXOD+Wb1aliyBE46Cfbtc3wBFRVw6JBjEjo0srWP8AfoFmiegO88A9A3/DORuWXb3p09vYN5nxnFCS8tYcXMFdzw+A18+9xv9+T+yYdZxlJKGMbQJqUCEJEwsAS4EJgMXCkikxMcukJVT3WXu91zzwLOBj4ATAE+BJzrHr8Q2K2qJ7jlPp7tzQwEixY5Qnvy5F6TzZtvOiGbna55vrILyiLQWQKrJuM7z0CsvCeegM/fnNjcMuHqRSxZAqec4iigf/l+fs0yXoVkYaWGMbQJ0gP4MLBZVbeoaidwP3BJwPIVqADKgHKgFNjl7vsccDOAqkZV9fV0Kj7QeCOBADZsgFNPhXvugddegw1fX0vDEedTc8Ax+8x9Bmr2w4oHyp2U0knKi0YdH8KBjtQzfCWbCSwXDlvveIQgE9MYhjF4CaIAxgCveL7vcLfFc5mI/FVEVonIOABVfRJYB7S7y8OquklEjnDPWSQi/yciK0Xk2EQXF5EGEVkvIuv37NkT8LZyT7JwTXAjcL70KLs/0MjM/bV8aiMc21nKpOu+z7qzjusnjOPLq1z1CMceOofKFOYWP7PM/A/Nz6pnkCjU9YWvrUXXW1ipYQxVchUF9ABQp6ofANYC9wKIyERgEjAWR2lME5GPASXutv9V1Q8CTwK3JSpYVZep6lRVnVpTU5Oj6qZPsnDNPtTXs27dfzD9mipaj1LqWZ1QGPcr77BW9lY6+YJiJDK3+Jll7vjzHVnNEZxQwV0yih3LchtWahhG8RBEAbwKjPN8H+tu60FV96pqLGHN3cDp7voMoEVV96nqPuC/gDOBvcABYLV73ErggxndwQCSLJTTS9AJ273lHTlrAd2kNrf0SfjWXUlZyDnurY63OGPsGZSHyzNy2AZWcIZhDBmCKICngONFZIKIlAFXAGu8B4jIaM/Xi4FN7vp24FwRKRGRUhwH8CZ1Yk8fAD7uHnc+sDHjuxggYpFAMYfs6tWJjwsaPeMtb8PX1zLvQ6nNLbHRvhXRGnT1co47cEHPvid3PMmh7t7Ecek6bIMquHTxOpYNwygeAo0DEJGLgB8BYeAeVV0sIjcC61V1jYjcjCP4I8AbwDxVbXUjiJYC5+A4hB9S1S+7ZdYC9wFHAHuAz6qqf7wkhc8F1A+f+YNjA6gOdPWaacrCZcw9bW7WIZSVlXG5g+adDEe3OvMOJCBXMfvx8yHH1keNSN1FSDeHkmEYucXmBM41TU3Q0OAMC45RVQXLlnHyO9/NOClbeztccQWsWJHY/NLeDl/5Cvz2t86lK2vaGfnZ2bx92JN0dB8EnPz/H6/7OM/sfIbrPnod//Hsf2QkuGP4zYecapL5fsrKpaICDh4MfHnDMLLEksHlgqYmqKtzssBdfXVf4Q/O96uvZm3VtRknZYtPC+0lphxKSvo7j2PCHwCB9x31PlbMXMG31n0r6zEDfj6Nn/3vyn5mHW/46X1PrOPwb5xMZY1zUHzklGEYhcV6AEGZPx9+8hMnAVAq3J4A9fUJd7e3OyklROA3v+lN/OZtLV9JEzexkPFsJ1TrmJfm/6men/4UamvhwgudDsi0VSfzZrgVlf69jWNHHJuTVA7t77Yze/VsR/i7ZYS1ku7Hv8ZRZ/+aDdc5vYvPrfkcu/btoqu7q6encLCzE13fQMUfl9DZaWYgwygE1gPIhqam4MIfnJ7AwoW+uxctgj//2XG2xlr63jDMK2niZzRQxzZCKGzbxv7ZDbx1VxPRqNOCXroUPvKR5M7jXKVySBh6Gu2Cjy3mjVAro/+5nmk/752L2NtTUIlQcfrKnDuWDcPIHusBBKGuDrZtS+8cEWdElQc/mzg45pw5c5yOw5ZoHbX0v952qaVW26iqcnoQt92WPEwzV85ob1K4inAlnV3dREOdPfMdh6NVdLuznpWGSikNlfZROMtnLGfm5JmBr2cYRm6xHkA2JEnm5kso5PQcPGzZ4gjucLh3WzjsZBrdurU3DHO8JL7eWN2eVox+rlI5eCeaue/S5bz/4Fw4cBSy6zTo7BX+ZaEyRMRyBxnGIMEUQBCSJHPzpbsbPvMZpydQVwdNTc7E8cc6u7yHHXusI8xj4wLeeE/i6705YnxappREM4QdXn447yl7T1r5guInmjlx81IuObSK8jEvQlmvsO+MdtLV3QUMrdxBNo7BGKqYAgjC4sW9ORLSIWZe27YNPvtZOPpolvwkxCvhOm4/o4lPfxomTOgVLLF8PP+ybzH7ibteVRXVP1mcchBan8u/O4oXblnCXz/jCO5PTf4UXdEudry9o29U0M9m90Y3ucoqGatXw8snLiASdXoXZaHeHNiloVLCEubHF/54yOQOShaZZRiDGfMBBCU26CtdX4AfCSKFvDH+lxxo4ruykLG6nejY8ZR8d7FvVJEf8QOwfCd4OQC7b/GcKOIor9ransFt8ezc50xY3/R8EwcjB+nq7qIkVEJIQnRHu1OOERgM2DgGY6hgPoBsqa93ZvVqbOzfG6iqgurq9MpLECnkzcfzm4p6JkgbX5gXpeSVtrSEv98k9i99K0FUUERY+mBcAd6eS0NDwh5BzCw07vBxRKNRFCVEiEPdh4hohKbnmwb9/MGpMsAaxmDHFEC61Nc7LffaWkeq1tY632+/PX0zUQLnci7y8fgJrqaHE4Rzos6sZX7EKar4OQe6urucJHShcg5FnTxEIUIcjBwc9PMHW4I8Y6hjJqBc4jUTxcwoyajNfK7gVHjnEo4NwGp+/8m8+LqbN6irEqQbSjqp2Q+7b4V1dbDgQlh7H2w6unf92P2CRKMJU0J0dHUQCoWIRD25iNzwUAg+6CxVCoxCEZsLuqHBeZ7t7cH8L4ZRTJgJaCCImYlU4b77ensJ1dWOJPZSVZVwprBsiUWsbNvWvyfxkb+thacbCHfUUPb75fDMXMr3H8Z9vy5nXR1MnwWtR0P9pb3ri87pjUqa98ACDhzqmxIiSrRH+IcI9RH+6Qw6K0ZHa3s77N0L3/pW6gywRhzetCkBAguMwmA9gIHCJ3NorkmUeTPZALRQCH52XhM/nPJZWkd2EQlDVSc9k9sfvR/+7tZGfkU9ocPbif7jbEpqW4i4sf8V4QrOrTuXh19ogdKDUNLZW3ikjPBf5xL5nb8SKGZHq2UxzZAkiRLz8c4bqbEeQKGJ9Q6iUeczR3+EWIu/oiKx47eysr9PIBSC44+HRx91egcPjqznkUXbOOd906iSsh7hX9ElnPj7BfwKp67RI1thTK/wB4gS5b1HvpcTjxuHhN2Rz12VECmDkk6OOCv5GIBidLT6OdErK/N3zSE11mDhwsSJEpOkRzEKgymAYidFVzpmOrn8cn9BGu/MBPjEJ+D88x2zxr98fx1n3XMWT+54kgPa24LvkFJaJnYh4igM+eQCCDv7w9He2chWbVzFY59by+QDDbC/htIHHfNSRbSGlZ/uHQMQE3LPPdcr7IrR0VoIpVSMJrCM8Rs5n8mIeiOvmAIoZmJd6W3bHL+CJywzvpW6fDn88pdOQytekLa3w8qVzsDkPtFFTU2s+8gopt89jbY32+iIeGwuGoKSTkZ+ZBXz5kFXF7B8LeFnHSF//jvLmfvB3oFeo0aM4oTNS5i/fzdP3TuT+bVLufD5vvMHx4RcfX1fYZevmcgyZSCVUiF6G3nHb+R8JiPqjfyiqoNmOf3003VYUVur6oj+vkttrb72muqsWapVVc6mqirVsWNVr75a9dlnVefPV50xwylm3jxVeW+zHrVwira/267NW5p1yk1jtb2mUqfMQ0u+hXI9GnY/uR6tWFShNbfUaPOWZlV1ypo/3yn7kn9t1sO+7ilrqbPuR0VF4tuILRUV+Xl8r72mes45qu3+VfPFe7/eZ5nLa8TOj/8d6+szL68oaGzsvaHYUlXlbB8uNDY6/18R57PA944ze2M/mVpwoZ7OMuwUgEhiiSmiqqrXXqsaCjkCNBRyBL1qr1AqL3dPqWtWvlmlfKtEuWqasrBKS76Nzr8IfW0EOu0qtOqbvcK/8juVunLDyoRVat7SrFWLq7TkhhKd9otpPevzH5zf77iYYli5vlkP/8YUrTi6vc9tVFbmV9jNm5f4ueTyevHXyAS/33FQU2QCcEApQgXopwAsCqiY8UtD7Y4f8ItRj0WvzJ4NkQj86siT0aOc+H9v6ubDD8KRB2HXYXCwtLf4ZCmjfdNJeOL9E40XiE0ME35oCd3djksjGnUmVvvFL3L72Pwii8Jh59+Yi6ieXEYvFWKsQbGOuxgSpPjfFgKbE3gwkmY4nW+454h2ZOZs9LiWnuydZd2AQmfYCdtXcb6HFKIh/wFc7e+286n7Z/PkKy1Ew4lz/vspiYpoDZOeX8HGugUcsWYtu7o3UXrJArZ/J705ilPR3g6rL2viH1qcXErbGc83WdwTzRQjm1DT+LmZg87RUCxYiGseCYUSDwJNMEfIQGFhoIMRv7QTPiGkiaJXxo6Fv69vpfx9LX1TN4ehKwQIhGLvpDjZPGuqarjuo9f1SfkQW299vZU/v9or/KF/zn+/mci+ft58npk8nUMjWtl1dj3Mmk7XyFZGX76ox+mZi3DI0c1NXPNUA+PVmVGtjm3cLQ3MKW3qeS7ZRvXkylE80OGffk7nWHSYkQOCOsGLYLCcKYCBJt0fPY3xA4mE0j/8A7x2Sm/q5sqSSkpCJQCEQyHKItDtTlAT7irj3Ld+6TuZ/NhZi5h2a29ZsXj/+Jz/CaeQ1G7u+PMdhMs6nVQU41yFFI5Q/sGVPcI4J+GQCxdSHukbh16lB7i+a2FOhXUuopcGOvwzvpEQm5zoiivydMEiEHIDTqL08fEj/5NE+A0kpgAGklz96An+VMlSQHgnhvna2V9zylCIRJXOkt5iuwUeeamZabcu4EBH35QPEY1w+NkrmbFvLeHnnFDQst8vZ8Jbcyk5VMNd03rj/f1mInur4y0+MvYMQloOpa6A7qrk3H1LmTAhh+GQPvHm49meU2Edm8AnkzQRyQbu5ZNYIyFmVYxNTnTvvXm4fpEIuYLgfZDV1f177sUyWC6RZ7hYl0EfBZQkrDMlsaiKWBRQXITBz89vDBRBMmXpFC25ocSJ+Pl/vZE/Petfq9FLr3pNz142Tau+U9UvMig+YuWkk/pHrrS/267zH5yvNbfU6MoNK/XiX16sXI+Grg9p6IZQ7zWvR0PfLtPyS+fps8/2DYcMhVQvvTTDiJ1snnMC/MJYswlfveoqp4ySkt4gkYEK/5wxwwkXvvDCPF8/6O8wlCKGgkYApYjwyzX4RAFZD2AgyXSEpLclBf0dTAcOMO2PCwO1JB+Z/QhnjD4H6a7sSdpGpAxe/nvYX4OsWkHkiFae2dPfhNO8tbnH7BGNOsuGDf1bsN4pJKsrq3n45YcJS5goUaLa1wkWDXVy6H2r+OlPe81X4bBT5osvZuhQDdIFT4NcjgyO2eCXL3e+R9wkqgcO9DdLpe0fCGhuWb3aibyqrXWec94GuwV534daLyFoy75IBsuZAhhIMv3RE71U8UXg/KlSCaeYE1fDveEvEoIjeR/Tnt7NvAvOY12l/2TyMbNHW1tqoRgLB410Rzis/LC+FYmUwWZH6bByBXfdBT/5iSOQYmaJDRsyNEuk6TxPRS5HBieywV90kRMOGy/o401OSRVCBoI07yOwg7zvxWIKyYRECjdoIy/HjZSMSdQtKNZl0JuAEnUPY13BZF1fv+6iZ2mjNuVAoooKVeZNcQaEXY+ysFL59zLlerTmlpqe4+JNOPMenNdnVHCMVAOY+pib4payG8t04hfn9RkBe+mljnkim1Gx+RjspRp8ZHAQkj23115ztiX6mcPhJL9vjs1eOSGIOWSATSE5w+/eqquD/w4DaPrCRgIXCSls+QlfAr8/t7tEQd+tqNa2mxpTpi2Y8Zl2DV88X/lqjZadulKP/+I8rf5uf+EehFRC8ZnNr+kRC6apXC99hH/JjY5SqPhWTT9BmO2o2EQjc3OiFHL4Z0323ObNcy5x/PH95UtSH0QuBGk+BFKqMotRcSUi/j78BH11dTAfwAD7PUwBFBt+L351df8XI1FrI9ESYLh5EAGbVGCm8eJevMBNQREn/E9acpLW3FKjH53d3E8QZtrSrqhQvZJG3UqtdiO6lVq9kkatqMhBuoYBGNqfLF9SvEJI2DPKVpD63eO8eYl/71wJsCJMm9CPoP+/2I+V6tkU4J6zUgDABcCLwGbgugT75wB7gGfd5RrPvluADcAm4Mf0jj5+zC0zds4xqeoxpBRAALNOnxejsdG/1ZHGHz6IgPXLofPmkmAvbo8wC2Bu6kPcH+eNOxsDt9zfuLNRO8J96xYF3U21Xkmjf8s5jti9PvusRwkOQCs1PilcKOQI/UcfdX6nuroUijuBUDkYqnJ+syD43WOiXuq8ecnfg3SVQ7FHAaXogSd9JxLdWwF6PRkrACAMvAy8FygDngMmxx0zB7gzwblnAX9yywgDTwIf114FMDXV9b3LkFIAmbxUQc5J1eVP8mfza4XGbM+vj/C5ftyLGxNmlTXtykXzla/V6Mc+v1Kv/s/EvoSeesUJlY6SKp1FY7CWe5Jns48qvZLGQD6FmPLrE946QHbqZL0zr+K++mqnLdDvPjy/7esjavs9u6Q9u6ANktgL4fceDIYWfYygiifos4m1/r3lJ3oWmf53syAbBXAm8LDn+zeAb8Qd46cAzgSeBiqBKmA9MElNAaTfrVQN9iIma0X4OaFdKRHfCo1fuvEXhPHCJW1bvo8A30ptz9dkLfdoimfTRm1qB7nP6VtJXLdMW2x+gjio+SuZSSvZmIWkprB0GiTJ3tPBZNMPqqjSeTZBzkumQOPrmKOeUTYKYCZwt+f7Z+KFvasA2oG/AquAcZ59twFvAW8Diz3bHwOed80/34qZhhJcv8FVHOvHjx+f8QMoKpJ1A5O8GG1+gii2lJY6TcN0HW6elotXcMfbnreJz/m1tf2ES9q2fB8BHgXdJrV651mNSVvuvr2TWDkiKR3ks2Y5Kaq9p1ZWqt55VqN2V/Y3L71ZUh3cxOIhiE8ikZKICfd4X8dVJY19zoufX8AvqqiPQk0WoZaOAMt1bylf5qF0FFXQxlr8uckaJamUT457UvlWANVAubv+eaDZXZ8I/B4Y4S5PAh9z941xPw8DHgGuSlWXQdkDiH+BE9lPky2eHz2xnVs0Co7gLytL/sIkeyHdP9evTzxWp8xDt713rH71sm9q6b9OUTmsXa8qadTdVDvX8iwx80p8cbNo1K6xtcH/uCkUYkdJ4pffKxT3keS5BmiBxpRfTL71Cbv08cH41SuRAO+ZnyGZIHbxi2a688z+99ld2bcO8b2vmTNVa2p6lZuvKSzIu+rnA4i9W0Fbt0HIpzkpXUWVygeXqF7JlMwAR0fl1QQUd3wYeNtd/yrwLc++bwNfS3BOQhNS/DLoFEA6rapkL4qHn5/fqG1u66+NWv35+e7+IC9MCiHb/P5yrfqmM0PYtKvQqoVo6NshPa/hfO2uSCBcq6v1jTsb+7Q455Q26l7pryhS/nGDtLISvPzeFu+VNOqeBEoq0bUTCehYr2XaNMcHMG1aXO8liZkqXognEuBB0j+kSjsRxA8T3/uaMKFvOYl6H77+AT9BlSycOdnzT6dFn09zUqZle+87puz87iPRO52qlx4jxz2pbBRACbAFmOBxAp8Ud8xoz/oMoMVdvxx41C2jFPgj8A/u96PdY0pds9G1qeoy6BRApnbVJMLS17QS5IVpbEzarfdOD+mdIazm6z42BPfPEmtxpmyFh8PJX/wUprFo/MvvHh91leFVJU4+pJ+fn1rIZBQa6vPsupEewZlq+sv4pc/1Gxu1a4xzP9uktsdxfemlqh/5iCOck/lh4knm1I83hQV+HkHj4RP91um26PPpfM9l7yKZUvPuC9JLj1EsPQDnXC4C/oYTDbTQ3XYjcLG7fjNOqOdzwDrg/e72MPBTnBDQjcAP3O3vwXEO/9U973YgnKoeg04BpBtZkY2dM9m4Ai+xkUbxLyE+00MuRFdOTlLvWqc38vqIWj00MkCYapA/m8+9vD6itveYBH/g7ooqvecTjUn9DX5CsbxcM+6We+uVyAYfP8I5HFa96CInmqenrgnuZx9OFFSfiKQ0BEOQ+YbTSnSXSeBCgOfnK9Ty7VDOhX/B75lUV2dn1ikWH0AxLYNOAaQTW52tXbOx0elexl+rrCzYKMTaWm2u6yv8uR4t+3d03uUjgt1HuovfH9lHGMZ8DRUVSZ6tW2afsQtxYwtmzOh1jIbDTs/l7bIEyiuAYy6K6O8n9G02J4qAShkVFSAKCnx6WUnen1TXjfkVtomjyJM63LOJh1fNzO6erRBMR8hnohCSPZNs014UQxRQMS2DTgH4vcB+oyuzJZ3ueIK6Tvln6TEBVS50hD/XozU3Hpae4zrokqwr75pDuuPMIXee5TqXU5Q5b57jiO4o6f/8v3Jcr9M6beexXw/K81wTmem82+75RKPuqqzt+5skMS/Ft+DjlVqy9ydlNFZj/2fUUeLzjqY7eDGeTFr02QjBdBRIpsom1TMJ4oMbgBBZUwCFIl9hbIkI8gdN8lK3/+JOnX/5CK35KrrywyN03qXlWvNVtPmMY3sFQkzIpiPo/eIQ04jMqajwEehxS5unxewXvx9rVV9Jo3bhE7Xip6SSxXb7PNc+DlY/BeKjvCMS1lluSoug/o20SKeXmiwKJpVDVHXgB4mlI3AzFc6pekXxPrgCDZIzBTAcCNpFDxLp4GfXTJWOIkAuo+Y6dMo/i7b/4k5t3tKsU5ZO0fZ3Ewf5e1uwqWL9FTRyZLXeeabTW/BzmHYjqVv+fs8qVWy3588cE/xXX+0R3n7nJ0oi5i4dJVXadlOj3vOJRj0YzrEASUeZJ6mjQvIIl3SiZ3JFOiaXVCHS6eT1Sfb+DGSD0IMpgOFAOk66TMLvUi0BBEDzBLRqoWjJ9SGd9otpWrW4SktuKNH5D87vc1yiMvwEenzYZyyFhN/AuTZqUw+q8xOuqZ6N5w8fHxLvO6I4JpQaG5PH0efDhJDOb+1NdBbkeG8eq0K0fHPRAwjir/MbI1BEKTBMAQxl4kPNYkLYT5ikejmzde4muYZ3jgDvlJM1t9SkFBRdY2oDX/v1EbXadlN/k9HBsBMplCx1RBS0u+o9iZVZKiUr4htZ4xvC6RVKyVqtAcJQM3p3go5V8QrOoO9IvhRXpveWjg8gyHOIL6MArfsgmAIYqiR7ydPtnsbItAcQ4BqvvfOa0/JPMN9wSkGRwGHZb9CX988bez6J/pQ+1+oirOsmB8h2maSlniinUjicpAfgTSKW7Bn47GujNr3xDPHPJN7hmyrbZzrvSBLFlc/kZ773mk4UUKp3axBhCmAwkMswNI/QTPtFTqU4vHb+ZL2MBNdo3tKsVYv7zhFQtqhM5z2YwDmaoIw7mZe8JZ1KuSW5x9g4gl2VKZ6p3zPyCMmY8zrm//7UpzSxDd+TjM+33Nhzqa7WDvoOJNpHld7BvJ7cQG2keG8CRDH11CMb23cAxZVVD2AgWtsFjNrJNaYAip1ch6F5BW+m4XdB7JrphMFpXxNQ5XcqtWyRZ46AAPUMZAYKKtD8hEjQFqvn/K4xtXrDib3x877hl0EEl9fOHleX7tJSfbusuicR3B3MCz4uIFkoZyZCLd706DfKNZc+gIGytw8Cu346mAIodnIdhuZtpaczBD2eVAIrWVc5wTWSzjecSFDEO5aTCf50W6uJRmtm+FtkPetYInzqEQsNTWpWShR9kqy3lqv0Cn7vSi5a7JmaNHN1Hb/3ZRBgCqDYydRO6ic04wV+0CRU6ZKLP4ufEzuR4sqkBZumkkrnntJKpZAuSXpXHSVV+p3JjcFyAwUx1xTCrJGuUkgn5j4bhpDpJ4YpgGInm5cu/o/kF6ufrxc43T9yUNNBNqF5XtI0U/XUMYAJIFG+nZ7Rytkq2yAhp0HemxTldCPp/Wa5aEBkYhZK53f0q2+Q+yik0zpPmAIodnJpJy3mFzgd52Gmg3PiybTlGFAppxytnGn6jyDPKomSig1ESx7yKrrv6hR2q3zE8WfS4Anak/Orb5DoJtWBb0DFkwcHtymAwUCufvhi7sKmEz7od2yStAsJydR2HFCRBhqtnG6vxVt3P9t9XJlR0Hcres1UMb+EX526CPfOJ5GMfLxPmTRSUpnmkvUaky3e96mxMXhSxXyQD2WragpgWJGnlygnpDOAKJngziQrZLpRHRkIvlRzE/cTXLFrxKdHSGUm87lOfPZQSJzwLkjq7JS/WTY9ykyVSjr+onSWVGbH+LTqucR7T8kc9Vk0Ck0BDDcGIk46E1LlEooXyslawNmGLgYZGJSmIg2SryjlEsSJ73NuN6LgzAKWM79EvuL4s2mkBDUlprMkMzvmy3ya7n1k2JAzBWAUHr/udTicNEIpWVK3AalzAIVRUeG0tBPNm5xTARXDRygn6gFkHZ6arx5lNo2UTEw9qZZkyjVf5tNM7iODupgCMApPht1rv8FfXWNr/eeyHWDeuLNRO+JH+eZDQMVIIJQPhqp0tjT2CP3jj1d99FGfeQBiZQSNlCm2HmU6prZUo9W9gnWgzaeZ5N3KoDdiCsAoPFmMdUg4aUljY34GYGVCOi25WI8n25ZfnFD++fmNyWcdiz83m0iZfBLvr6mu7h9FlY4pMYiZJd7sOFDKLlmgQ9D3IACmAIzCk0X3+gdTG/X1EbUaFWf+4Vk0JiwqJwOwUpCw1xG0JRdEKCXyAQQQwiln//KSruAZiAgyP0d90GcUZD4CvwR4+VZwyXpb6SRzNB+AMWjJ4QsdZMLzfJGw15GsB+A3AYo36iRZFFA+BFS6pod8jyHJxKmbaPKhYiTVe5/st87Re2AKwCgO/ML4Mni5U060nmOSpn0o5tDbRBS6BxAv2DIxiRXDwMYgFMG4HD8FEMIw8klTE9TVQSjkfAK0tcF998HBg7B3r/N32LYNGhqc4wOyaxdcey20tDifO3fm4wZ62bIFZs2Cqirne1UV1NfD1q04K8uWQW0tiDify5Y524uRxYt7byRGVZXzGyTavnhx7q7d1ORcZ9u23t9+7970yxk/3r987zuXxjuVF7ZvT2/7QJJIKxTrYj2AQUayVnERtIoyYaB7HXklm3w52ZCLEM50ZvYqdE+sCN51zARkDDjJXvzBMktUHGk5WgtJsYVteknX/5AoCiiHo7fzThEoJVMAxsCTTMgX6o9aBH/GvFPs95hsPEi2SqtYEyEWWCGbAjAGnmRCvlBCqhhbiLmm0NksU5Hr3z5ILp1iufdk5FFJmAIwBp5swt/yRbG2EHNFY2Pi+yu2e8zVb5/uIK9iJc8NIj8FIM6+wcHUqVN1/fr1ha6GkQ5NTbBwoRPxMH68E01SyMiYujon6iSe2lonOmmw43d/MHTu0Yvf/YbDEI0WxzsXhDy/lyLytKpOjd9uYaBGfqmvd17gaNT5LPQf0S/80RvmWGxhhOmQLLQwl6GcxYLf/UajxfPOBaFAoaKBFICIXCAiL4rIZhG5LsH+OSKyR0SedZdrPPtuEZENIrJJRH4sIhJ37hoReSH7WzGMAKSK108Uo57m+ISC4hcbX109OARhuvjdr9/2RBSDws/FfWRCIruQdwHCwMvAe4Ey4Dlgctwxc4A7E5x7FvAnt4ww8CTwcc/+S4FfAi+kqoeaD8AYCAa7k7jYI4ByTT7mFSjE8yqQDyBID+DDwGZV3aKqncD9wCVB9QtQ4SqOcqAU2AUgIiOALwPfCViWYeSfTLrixdCCjDHYRiRnS7b3u3AhHDjQd9uBA872gaRAv1tKJ7CIzAQuUNVr3O+fAc5Q1S94jpkD3AzsAf4GfElVX3H33QZcAwhOL2Ghu/2HwH8DzwAPquoUn+s3AA0A48ePP32bn4PLMHJBus64mMnIK0Sqqoa20B1KhEJOezseEceHMETItxP4AaBOVT8ArAXudS86EZgEjAXGANNE5GMicirwPlX9TaqCVXWZqk5V1ak1NTU5qq5h+BDESeylWFqQRmYUyvZeJARRAK8C4zzfx7rbelDVvap6yP16N3C6uz4DaFHVfaq6D/gv4Ex3mSoibcATwAki8limN2EYOSPdrngxJ/oyUpOuwh9iBFEATwHHi8gEESkDrgDWeA8QkdGerxcDm9z17cC5IlIiIqXAucAmVb1LVY9T1Trgo8DfVPXj2d2KYeSIdEJXi70FWUz+iWJkuPlM4ihJdYCqRkTkC8DDOJE896jqBhG5EcezvAZYICIXAxHgDZyoIIBVwDTgeRyH8EOq+kDub8MwCsTixYl9AMXQgoz3T8RCWmHYCLhA1NcP2+dhI4ENI1uKbbRzjKE+6tkIjJ8T2BSAYQxVhkmEi5EaSwVhGMONYvdPGAXHFIBhDFWGeYSLkRpTAIYxVBnmES5GalJGARmGMYgZxhEuRmqsB2AYhjFMMQVgGIYxTDEFYBiGMUwxBWAYhjFMMQVgGIYxTBlUI4FFZA/gHdt+NPB6gapT7NizSY49H3/s2fgzWJ9Nrar2y6c/qBRAPCKyPtHwZsOeTSrs+fhjz8afofZszARkGIYxTDEFYBiGMUwZ7ApgWaErUMTYs0mOPR9/7Nn4M6SezaD2ARiGYRiZM9h7AIZhGEaGmAIwDMMYpgxaBSAiF4jIiyKyWUSuK3R9BhoRGSci60Rko4hsEJF/dbcfJSJrReQl9/NId7uIyI/d5/VXEflgYe8g/4hIWESeEZEH3e8TROTP7jNYISJl7vZy9/tmd39dQSueZ0TkCBFZJSKtIrJJRM6098ZBRL7k/p9eEJFfiUjFUH5vBqUCEJEwsAS4EJgMXCkikwtbqwEnAvybqk4GPgL8s/sMrgP+qKrHA390v4PzrI53lwbgroGv8oDzr8Amz/fvAT9U1YnAm8Bcd/tc4E13+w/d44YytwMPqer7gVNwntGwf29EZAywAJiqqlOAMHAFQ/m9UdVBtwBnAg97vn8D+Eah61XgZ/I74O+AF4HR7rbRwIvu+k+BKz3H9xw3FBdgLI4gmwY8CAjOCM6S+HcIeBg4010vcY+TQt9Dnp7L4cDW+Puz90YBxgCvAEe578GDwP83lN+bQdkDoPeHirHD3TYscbuepwF/Bo5V1XZ3107gWHd9uD2zHwFfA2Kzn1cDb6lqxP3uvf+eZ+Puf9s9figyAdgD/IdrHrtbRN6DvTeo6qvAbcB2oB3nPXiaIfzeDFYFYLiIyAjg18AXVfUd7z51mibDLs5XRKYDu1X16ULXpQgpAT4I3KWqpwH76TX3AMP6vTkSuARHSR4HvAe4oKCVyjODVQG8CozzfB/rbhtWiEgpjvBvUtXV7uZdIjLa3T8a2O1uH07P7GzgYhFpA+7HMQPdDhwhIrFpUL333/Ns3P2HA3sHssIDyA5gh6r+2f2+Ckch2HsDnwC2quoeVe0CVuO8S0P2vRmsCuAp4HjXO1+G46hZU+A6DSgiIsDPgU2q+gPPrjXA1e761Ti+gdj2q9yojo8Ab3u6/EMKVf2Gqo5V1Tqcd6NZVeuBdcBM97D4ZxN7ZjPd44dkC1hVdwKviMiJ7qbzgY3YewOO6ecjIlLl/r9iz2bovjeFdkJk4bC5CPgb8DKwsND1KcD9fxSnm/5X4Fl3uQjHBvlH4CXgUeAo93jBiZx6GXgeJ9Kh4PcxAM/p48CD7vp7gb8Am4GVQLm7vcL9vtnd/95C1zvPz+RUYL377vwWONLem55ncwPQCrwA3AeUD+X3xlJBGIZhDFMGqwnIMAzDyBJTAIZhGMMUUwCGYRjDFFMAhmEYwxRTAIZhGMMUUwDGoEVEbnUzN94qIteKyFUDcM05InKc5/vdxZyIUETqROSFQtfDKE5KUh9iGEVLA068encuCxWRcJIy5+DEiL8GoKrX5PLaxYaIlGhvHhxjiGE9AGPAEZGr3Nzyz4nIfe62OhFpdrf/UUTGu9t/4eaj/18R2SIiM93ta4ARwNMicrmIXC8iX3H3fcgt51m3d/CCu32OiNzpqceDIvJxd32fiHxfRJ4DzhSRb4vIU25e+GXuSNiZwFSgyS27UkQeE5GpbhlXisjz7jnf81xnn4gsdu+3RURiida8z+R6EbnHLW+LiCzwPJcXPMd9RUSud9cfE5Efish6cfL6f0hEVouT0/87nuJLRKTJPWaViFS5558uIo+LyNMi8rAnFcRjIvIjEVmPk1LbGKKYAjAGFBE5Cfh3YJqqnkKvgLkDuFdVPwA0AT/2nDYaZ+TzdOC7AKp6MXBQVU9V1RVxl/kP4POqeioQtHfwHuDPqnqKqj4B3KmqH1InL3wlMF1VV+GMoK13r3vQc1/H4eSDn4Yz0vZDIvKPnrJb3Pv9b+CffOrwfpz0wx8G/p84uZ5S0amqU4Gf4KQo+GdgCjBHRGKZKU8ElqrqJOAdYL5b9h3ATFU9HbgHWOwpt0xVp6rq9wPUwRikmAIwBpppwEpVfR1AVd9wt58J/NJdvw9H4Mf4rapGVXUjvWmKEyIiRwCHqeqT7qZfJjncSzdOYr0Y54kzy9Pzbp1PSnH+h4DH1EkkFsFRYue4+zpxcsuDk164zqeM36vqIffZ7CbFvbrEcmA9D2xQ1XZVPQRsoTeJ2yuq+id3vRHn2Z6IoyjWisizOEp5rKfceKVqDEHMB2AMBg551iWLciL0bfRUeNY7YnZ/EakAluLkvXnFNbl4j02XLu3NudKN///Oe5+x45LV2XtONO78qOc68fleFOc5blDVM33qst9nuzGEsB6AMdA0A5+KmSdE5Ch3+//iZO4EqAf+J5PCVfUt4F0ROcPddIVndxtwqoiERGQcjqklETEh+7o48y3M9Ox7FzgswTl/Ac4VkaPFmbL0SuDxTO4hjl3AMSJSLSLlOGawdBkvIjFBPwt4Amdmr5rYdhEpdc1zxjDCegDGgKKqG0RkMfC4iHQDz+BE1vwLzixVX8WZseqzWVxmLvAzEYniCOG33e1/wpkOcSPOPLj/51PHt0TkZzjRPjtx0o/H+AXwExE5iGO2ip3TLiLX4aQOFhxzzu/IElXtEpEbcRTMqziZKtPlRZw5o+/Bufe7VLXTdWr/WEQOx5EFPwI2ZFtnY/Bg2UCNIYeIjFDVfe76dThz2Fo0i2HEYT0AYyjySRH5Bs77vQ2nh2EYRhzWAzAMwximmBPYMAxjmGIKwDAMY5hiCsAwDGOYYgrAMAxjmGIKwDAMY5jy/wOhe7o9xFf9ZQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plaquette_0p01_conf = [int(a.split(\" \")[0]) for a in filter(lambda x: x!=\"\",open(\"plaquette_0p01.dat\").read().split(\"\\n\"))]\n", "plaquette_0p01 = [float(a.split(\" \")[1]) for a in filter(lambda x: x!=\"\",open(\"plaquette_0p01.dat\").read().split(\"\\n\"))]\n", "\n", "plaquette_0p02_conf = [int(a.split(\" \")[0]) for a in filter(lambda x: x!=\"\",open(\"plaquette_0p02.dat\").read().split(\"\\n\"))]\n", "plaquette_0p02 = [float(a.split(\" \")[1]) for a in filter(lambda x: x!=\"\",open(\"plaquette_0p02.dat\").read().split(\"\\n\"))]\n", "\n", "plaquette_0p04_conf = [int(a.split(\" \")[0]) for a in filter(lambda x: x!=\"\",open(\"plaquette_0p04.dat\").read().split(\"\\n\"))]\n", "plaquette_0p04 = [float(a.split(\" \")[1]) for a in filter(lambda x: x!=\"\",open(\"plaquette_0p04.dat\").read().split(\"\\n\"))]\n", "\n", "fig, ax = plt.subplots()\n", "\n", "#plt.xlim()\n", "\n", "def ac(x, dt):\n", " m = sum(x) / len(x)\n", " var = sum([ (x[i] - m)**2. for i in range(len(x)-dt) ])\n", " cov = sum([ (x[i] - m)*(x[i+dt] - m) for i in range(len(x)-dt) ])\n", " return cov / var\n", "\n", "print(ac(plaquette_0p04[60:],2))\n", "print(ac(plaquette_0p02[60:],2))\n", "\n", "plt.title(\"Plaquette\")\n", "plt.xlabel(\"configuration number\")\n", "ax.plot(plaquette_0p02_conf[6:],plaquette_0p02[6:],marker=\"*\",ls=\"\",color=\"blue\",label=\"m = 0.02\")\n", "ax.plot(plaquette_0p04_conf[6:],plaquette_0p04[6:],marker=\"o\",ls=\"\",color=\"red\",label=\"m = 0.04\")\n", "ax.plot(plaquette_0p01_conf[6:],plaquette_0p01[6:],marker=\"X\",ls=\"\",color=\"green\",label=\"m = 0.01\")\n", "plt.legend()\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "61853c2f", "metadata": {}, "source": [ "We note that currently the $m=0.01$ ensemble is not yet thermalized. We therefore need to generate more configurations before obtaining reliable results for this mass.\n", "\n", "For the other ensembles, one may perform an analysis for configuration numbers beyond 400. The code to compute the mesonic correlation functions is given below. We note, that $\\Gamma=\\gamma_5$ has the correct quantum numbers, however, couples strongly to the excited states. We therefore combine the $\\Gamma$ with a **covariant Gaussian smearing** operator $S_\\Delta$,\n", "$$\n", "\\begin{align}\n", "\\Gamma^\\prime = \\gamma_5 S_\\Delta \\,,\n", "\\end{align}\n", "$$\n", "where\n", "$$\n", "\\begin{align}\n", "\\Delta \\psi(x) &= \\sum_{\\mu \\in \\{0,1,2\\}} (U_\\mu(x) \\psi(x+\\hat{\\mu}) + U_\\mu^\\dagger(x-\\hat{\\mu}) \\psi(x-\\hat{\\mu}) - 2 \\psi(x)) \\,, \\\\\n", "S_\\Delta(\\sigma,n) &= \\left( 1 + \\frac{\\sigma^2}{4 n} \\Delta \\right)^n \\,.\n", "\\end{align}\n", "$$\n", "One finds in practice that this leads to a suppressed coupling to the excited state relative to the ground state.\n", "The code to perform this contraction is shown below." ] }, { "cell_type": "markdown", "id": "06a73834", "metadata": {}, "source": [ "Another very useful method is the use of random insertions of the identity matrix as an outer product of two vectors.\n", "Consider a vector $\\eta_x$ with elements in the complex numbers satisfying\n", "$$\n", "\\begin{align}\n", "\\langle \\eta_x \\eta_y^* \\rangle &= \\delta_{xy} \\,,\n", "\\end{align}\n", "$$\n", "where $\\langle \\cdot \\rangle$ denotes the statistical average.\n", "We can then write this as\n", "$$\n", "\\begin{align}\n", "\\langle \\eta \\eta^\\dagger \\rangle &= \\mathbb{1}\n", "\\end{align}\n", "$$\n", "such that\n", "$$\n", "\\begin{align}\n", "C(t) &= -\\langle \\text{Tr}[ D^{-1}[U,m_l] \\Gamma D^{-1}[U,m_l] \\Gamma' ] \\rangle_U \\\\\n", "&= -\\langle \\text{Tr}[ D^{-1}[U,m_l] \\eta \\eta^\\dagger \\Gamma D^{-1}[U,m_l] \\Gamma' ] \\rangle_U \\\\\n", "&= -\\langle \\text{Tr}[ D^{-1}[U,m_l] \\eta \\Gamma \\eta^\\dagger D^{-1}[U,m_l] \\Gamma' ] \\rangle_U \\\\\n", "&= -\\langle \\text{Tr}[ D^{-1}[U,m_l] \\eta \\Gamma \\eta^\\dagger \\gamma_5 (D^{-1}[U,m_l])^\\dagger \\gamma_5 \\Gamma' ] \\\\\n", "&= -\\langle \\text{Tr}[ \\phi \\Gamma \\gamma_5 \\phi^\\dagger \\gamma_5 \\Gamma' ] \\rangle_U\n", "\\end{align}\n", "$$\n", "with\n", "$$\n", "\\begin{align}\n", "\\phi = D^{-1}[U,m_l] \\eta \\,.\n", "\\end{align}\n", "$$\n", "In this way we can factor the correlation function and statistically generate a full volume average. This comes at the cost of introducing additional noise due to the $\\eta$ fluctuations and depending on if the gauge noise or this $\\eta$ noise is larger, it may be a better statistical estimator.\n", "\n", "For the mesons to be computed in this chapter, this method performs well and for concreteness, we use $\\eta_x \\in [1,-1]$, where both are drawn with equal probability. We call this a $Z_2$ source. The corresponding GPT code is given below." ] }, { "cell_type": "code", "execution_count": 102, "id": "d262a550", "metadata": {}, "outputs": [], "source": [ "import gpt as g\n", "\n", "g.default.set_verbose(\"io\", False)\n", "g.default.set_verbose(\"cg_convergence\", False)\n", "g.default.set_verbose(\"cg\", True)\n", "\n", "def mesons(mass, config):\n", " mass_fmt = str(mass).replace(\".\",\"p\")\n", " U = g.load(f\"/glurch/home/lec17310/configurations/DWF/IWASAKI/2f/16c32/ls16/mu{mass_fmt}/ckpoint_lat.{config}\")\n", " D_m = g.qcd.fermion.mobius(U, M5=1.8, mass=mass, Ls=16, b=1., c=0.,\n", " boundary_phases=[1,1,1,-1])\n", "\n", " grid = U[0].grid\n", "\n", " inv = g.algorithms.inverter\n", " pc = g.qcd.fermion.preconditioner\n", " cg = inv.cg({\"eps\": 1e-7, \"maxiter\": 1000})\n", " invD = D_m.propagator(inv.preconditioned(pc.eo1_ne(), cg))\n", "\n", " src = g.mspincolor(grid)\n", " rng = g.random(str(mass) + \"-\" + str(config))\n", " t0 = rng.uniform_int(min = 0, max = grid.gdimensions[3] - 1)\n", " g.create.wall.z2(src, t0, rng)\n", "\n", " g.message(\"smear start\")\n", " smear = g.create.smear.gauss(U, sigma=3.0, steps=30, dimensions=[0, 1, 2])\n", "\n", " src = g(smear * src)\n", " g.message(\"smear end\")\n", "\n", " prop_xy = g(invD * src)\n", " prop_yx = g(g.gamma[5] * g.adj(prop_xy) * g.gamma[5])\n", "\n", " pion_correlator = np.array(g.slice(g.trace(prop_xy * g.gamma[5] * prop_yx * g.gamma[5]),3))\n", " pion_correlator = np.roll(pion_correlator, -t0)\n", "\n", " rho_correlator = sum([np.array(g.slice(g.trace(prop_xy * g.gamma[i] * prop_yx * g.gamma[i]),3)) for i in range(3)])\n", " rho_correlator = np.roll(rho_correlator, -t0)\n", "\n", " return {\n", " \"pion\" : pion_correlator,\n", " \"rho\" : rho_correlator\n", " }\n", "\n", "#mesons_0p04 = {}\n", "#for c in range(100,260,10):\n", "# g.message(c)\n", "# mesons_0p04[c] = mesons(0.04, c)" ] }, { "cell_type": "code", "execution_count": 4, "id": "02c03b16", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "dict_keys([400, 410, 420, 430, 440, 450, 460, 470, 480, 490, 500, 510, 520, 530, 540, 550, 560, 570, 580, 590])\n" ] } ], "source": [ "array = np.array\n", "# only use eval for trusted files!\n", "mesons_0p04=eval(open(\"mesons_0p04.json\").read())\n", "mesons_0p04={x: mesons_0p04[x] for x in mesons_0p04 if x >= 400} # keep only thermalized\n", "\n", "mesons_0p02=eval(open(\"mesons_0p02.json\").read())\n", "mesons_0p02={x: mesons_0p02[x] for x in mesons_0p02 if x >= 400}\n", "\n", "print(mesons_0p02.keys())" ] }, { "cell_type": "markdown", "id": "7a46bfb1", "metadata": {}, "source": [ "Let us now compute the pion and $\\rho$ energies for different quark masses." ] }, { "cell_type": "code", "execution_count": 104, "id": "e83dd050", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "45\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAX70lEQVR4nO3dfXAd9X3v8ffHj0I2TxfUjm3ZknJrwMY8JNa4askNpqQMpIAbmulARMhtCZ4EaNLAJLhxkt5w67lA79Cna9oK0oQLimlKXGpap3RaHnMHMpbDw8V23PhS2ZbtgHGNG5CNsf29f+yRLYkj6Ug65+ye1ec1c0Znf7s6+9Xq6LN7frv7kyICMzOrfZPSLsDMzMrDgW5mlhMOdDOznHCgm5nlhAPdzCwnHOhmZjnhQDczywkHutkIJH1S0nZJ70h6TNJ/GmbZCyVtlNRb+HphkWWmSdoiqaeihduE40A3G4akc4G/BD4F/DzQC9w3xLLTgL8DHgZOBx4E/q7Q3t+XgL2VqtkmLge61QxJ3ZK+JOmVwtHyNyX9vKTvS/qZpH+WdHqZV9sOPB4Rz0bE28DXgGsknVxk2aXAFOCPI+LdiPhTQMCv9PsZWoDrgf9R5jrNHOhWc34D+FXgLOAq4PvAV4AGkvfz54t9k6R5kt4a5vHJIdZ3LvBy30RE/D/gcGH9xZZ9JQaOp/FKob3PnxXqPVjCz2o2KlPSLsBslP4sIl4HkPQc8EZEvFiY/lvg0mLfFBE7gNPGsL6ZwIFBbQeAYkfowy4r6ePA5Ij4W0lLx1CL2bAc6FZrXu/3/GCR6ZllXt/bwCmD2k4BfjaaZSXNAO4BPlbm+syOc5eLTQiFLpe3h3m0D/Gtm4AL+r3OB4DpwL8Osez5ktSv7fxC+3ygGXhO0k+BtcAsST+V1DzuH9AMH6HbBFHochnL0Xsn8Lyk/wL8CLgTWBsRxY7QnwaOAp+X9BfATYX2J4FjwNx+y/4y8L+AD+ErXqxMfIRuNoyI2AR8liTY3yDpD7+5b37hCpuvFJY9DPw6cAPwFvDbwK9HxOGIOBIRP+17AP8OHCtMH63mz2T5Jf+DCzOzfPARuplZTowY6JL+StIbkl4dYr4k/amkbYUbPj5U/jLNzGwkpRyhfxu4fJj5V5CcwZ8PLAf+fPxlmZnZaI0Y6BHxLMkJnKEsA/53JF4ATpM0q1wFmplZacpx2eIcYGe/6Z5C257BC0paTnIUz4wZMxafc845ZVi9mdnEsXHjxjcjoqHYvKpehx4RHUAHQGtra3R1dVVz9WZmNU/S9qHmleMql10MvGGisdBmZmZVVI5AXwfcULjapQ04EBHv624xM7PKGrHLRdIaknGezyz8h5XfB6YCRMRfAOtJBhzaRjL4/29VqlgzMxvaiIEeEdeNMD+AW8pWkZmZjYnvFDUzywkHuplZTjjQzcxywoFuZpYTDnQzs5xwoJuZ5YQD3cwsJxzoZmY54UA3M8uJ2gr0zk5oboZJk5KvnZ1pV2RmlhlVHT53XDo7Yfly6O1NprdvT6YB2tvTq8vMLCNq5wh95coTYd6ntzdpNzOzGgr0HTtG125mlkVLlyaPCqidQJ83b3TtZmYTTO0E+qpVUF8/sK2+Pmk3M7MaCvT2dujogKYmkJKvHR0+IWpmVlA7V7lAEt4OcDOzomrnCN3MzIblQDczywkH+lj4jlUzy6Da6kPPAt+xamYZ5SP00fIdq2aWUQ700fIdq2aWUQ700fIdq2aWUQ700fIdq2aWUQ700fIdq2Y2Vp2d8MIL8MwzFblCzle5jIXvWDWz0eq7Qu7dd5PpClwh5yN0M6usCg4XW1OqcIWcA90srxyk2VKFK+Qc6LXMd6yalS7tHVwVrpBzoNeqvv647dsh4kR/3EQO9c5OqKtLTlZP9B1chU++1VwdWVCNK+QiIpXH4sWLw8ahqSkiifKBj6amtCtLx8MPR9TXD9wW9fVJ+0STlW2RlTr6XHxx8kjTww9HTJ9+4m91DNsC6IohclXJ/OprbW2Nrq6uVNadC5MmJX8ig0lw7Fh1a+nshBtvTM7ez54Nt90GV11V3RouuQR2735/++zZ8NRT1a0lbVnZFlmpo8+nPpV8feih6q97cB1S8sllDCRtjIjWYvN82WKtmjcv6WYp1l5Ngy/F2r0bvvY1mD4drrmmenXs2TN0+8yZ1atj7Vq4/XY4fBjmzIEVK6q7HSA72yIrdfSZPDn5msa6B9dx9GhFXtp96LUqK3esFrsU6+BBuOuu6tYxe/bo2ith7Vr48peTMAfYtSuZXru2ejVANrZFluqA5HewcSM8/zwsWVL930mVlBToki6XtFXSNkkrisyfJ+kpSS9KekXSx8pfqg2QlTtWh7rkqthH7UpasQJOOmlg20knJe3Vctddyc6svzR2blnYFlmqIys72ioYMdAlTQZWA1cAC4HrJC0ctNhXge9GxAeBa4H7yl2oFdHeDt3dSZ95d3c6d68O1cVT7aOwa66Be+5Jujmk5Os991S3u2OonVi1d25922LatGQ6jW2RpTqysqOtglL60JcA2yLiNQBJjwDLgM39lgnglMLzU4Eqv4MtNatWDfyHH5BcOvjFL8Lbb1e3lssuSx79VbOGWbOKh/esWelsiwsvTJ73nQSsdg1ZqWO4HW0a2+To0eSgowJKCfQ5wM5+0z3ALw5a5r8B/yTpd4AZwEeLvZCk5cBygHkebjYf+j4VrFyZdL/Mm5eE/EQc6+aee96/c6uvT9rPOqv69fR1d6Sx7izVMdwFBGn+XiqgXCdFrwO+HRGNwMeAhyS977UjoiMiWiOitaGhoUyrttRloesnC/rOa0yfnkynPRLn008nj4kuKxcQVEEpR+i7gLn9phsLbf3dCFwOEBHPS6oDzgTeKEeRZjXDI3G+X9o7lb7fR9+9Ek1Nuf0UWUqgbwDmS2ohCfJrgU8OWmYHcCnwbUkLgDpgbzkLNTMbs/Z2uP/+5HnaO5gKGrHLJSKOALcCTwBbSK5m2STpTklXFxa7HbhJ0svAGuC/Rlq3oJqZTVAl3SkaEeuB9YPavt7v+WbgovKWZmZmo+Fb/83MqqmCXT4OdDObGHLcd97HY7mYmeWEA93MLCcc6GZmOeFANzPLCQe6mVlOONDNzHLCgW5mlhMOdDOznHCgm5nlhAPdzCwnHOhmZjnhQDczywkHuplZTjjQzcxywoFuZpYTDnQzs5xwoJuZ5YQD3cwsJxzoZmY54UA3M8sJB7qZWU440M3McsKBbmaWEw50M7OccKCbmeWEA93MLCcc6GZmOeFANzPLCQe6mVlOONDNzHLCgW5mlhMOdDOznCgp0CVdLmmrpG2SVgyxzG9K2ixpk6TvlLdMMzMbyZSRFpA0GVgN/CrQA2yQtC4iNvdbZj7we8BFEbFf0s9VqmAzMyuulCP0JcC2iHgtIg4DjwDLBi1zE7A6IvYDRMQb5S3TzMxGUkqgzwF29pvuKbT1dxZwlqT/I+kFSZcXeyFJyyV1Serau3fv2Co2M7OiynVSdAowH1gKXAfcL+m0wQtFREdEtEZEa0NDQ5lWbWZmUFqg7wLm9ptuLLT11wOsi4j3IuLfgH8lCXgzM6uSUgJ9AzBfUoukacC1wLpByzxGcnSOpDNJumBeK1+ZZmY2khEDPSKOALcCTwBbgO9GxCZJd0q6urDYE8A+SZuBp4AvRcS+ShVtZmbvV1IfekSsj4izIuI/R8SqQtvXI2Jd4XlExG0RsTAizouIRypRbGcnNDfDpEnJ187OSqzFzKw2jXgdelZ0dsLy5dDbm0xv355MA7S3p1eXmVlW1Myt/ytXngjzPr29SbuZmdVQoO/YMbp2M7OJpmYCfd680bWbmU00NRPoq1ZBff3Atvr6pN3MzGoo0NvboaMDmppASr52dPiEqJlZn5q5ygWS8HaAm5kVVzNH6GZmNjwHuplZTjjQx8B3rJpZFtVUH3oW+I5VM8sqH6GPku9YNbOxqvSnewf6KPmOVbPak4Vu0r5P99u3Q8SJT/flrMWBPkpZumM1C2/SLNVhVkw1grQUVfl0HxGpPBYvXhy16OGHI+rrI5K3RvKor0/aXUd6dfTV0tQUISVf06jBBsrC76SpaeD7s+/R1FTdOqTidUijex2gK4bIVQf6GPhNmr06srJjycJ7Iyuy8jspV5COV7n+VoYLdCXzq6+1tTW6urpSWXceTJqUvB0Gk+DYsWzU8eMfV6+OSy6B3bvf3z57Njz1VHVqePxx+OpX4dChE211dfAHfwBXXVWdGvrXcu+9sGcPzJoFt91W/Rqy8DvJUh3F3h/19aMfwkTSxohoLTbPfeg1Kit9+UOtb9as6taxZ8/o2ivh3nsH/rFCMn3vvdWrAU4Ex+7dyc529+5k+vHHq1tHFn4nkOzM6uoGttXVJe3VdNVVyc599uwKjkc11KF7pR+13OWSBVn5OJuVOrLQ9ZO3j/Z5qSMiX11huA89n7LyJs1CHVnYsWQlwLKyY8nC7ySPHOg2IaS9Y8lKgGVlxxKR/u8kj4YLdJ8UNSujzs7kuuIdO5LzC6tWVX9IiMHDU8DYTr5ZNg13UtRjuZiVURbG7O9bf9o7Fqs+B7pZDmVhx2LV58sWzcxywoFuZpYTDnQzs5xwoJuZ5YQD3cwsJxzoZmY54UA3M8sJB7qZWU440M3McsKBbmaWEw50M7OcKCnQJV0uaaukbZJWDLPcb0gKSUVHAjMzs8oZMdAlTQZWA1cAC4HrJC0sstzJwBeAH5a7SDMzG1kpR+hLgG0R8VpEHAYeAZYVWe6/A3cDh4rMMzOzCisl0OcAO/tN9xTajpP0IWBuRPzDcC8kabmkLklde/fuHXWxZmY2tHGfFJU0CbgXuH2kZSOiIyJaI6K1oaFhvKs2M7N+Sgn0XcDcftONhbY+JwOLgKcldQNtwDqfGDUzq65SAn0DMF9Si6RpwLXAur6ZEXEgIs6MiOaIaAZeAK6OCP/DUDOzKhox0CPiCHAr8ASwBfhuRGySdKekqytdoJmZlaak/ykaEeuB9YPavj7EskvHX5aZmY2W7xQ1M8sJB7qZWU440M3McsKBbmaWEw50M7OccKCbmeWEA93MLCcc6GZmOeFANzPLCQe6mVlOONDNzHLCgW5mlhMOdDOznHCgm5nlhAPdzCwnHOhmZjnhQDczywkHuplZTjjQzcxywoFuZpYTDnQzs5xwoJuZ5YQD3cwsJxzoZmY54UA3M8sJB7qZWU440M3McsKBbmaWEw50M7OccKCbmeWEA93MLCcc6GZmOeFANzPLCQe6mVlOlBToki6XtFXSNkkrisy/TdJmSa9I+hdJTeUv1czMhjNioEuaDKwGrgAWAtdJWjhosReB1og4H3gUuKfchZqZ2fBKOUJfAmyLiNci4jDwCLCs/wIR8VRE9BYmXwAay1ummZmNpJRAnwPs7DfdU2gbyo3A94vNkLRcUpekrr1795ZepZmZjaisJ0UlXQ+0An9YbH5EdEREa0S0NjQ0lHPVZmYT3pQSltkFzO033VhoG0DSR4GVwMUR8W55yjMzs1KVcoS+AZgvqUXSNOBaYF3/BSR9EPhL4OqIeKP8ZZqZ2UhGDPSIOALcCjwBbAG+GxGbJN0p6erCYn8IzAT+RtJLktYN8XJmZlYhpXS5EBHrgfWD2r7e7/lHy1yXmZmNUkmBXi3vvfcePT09HDp0KO1SKqauro7GxkamTp2adilmljOZCvSenh5OPvlkmpubkZR2OWUXEezbt4+enh5aWlrSLsfMciZTY7kcOnSIM844I5dhDiCJM844I9efQMwsPZkKdGDUYb50afKoFXndWZlZ+jIX6GZmNjY1HeidnfDCC/DMM9DcnExXymc+8xk2b95cuRWYmY1Tpk6KjkZnJyxfDu8W7kndvj2ZBmhvL//6HnjggfK/qJlZGdXsEfrKldDbO7CttzdpH4/u7m7OOecc2tvbWbBgAZ/4xCfo7e1l6dKldHV1AbBmzRrOO+88Fi1axB133HH8e2fOnMnKlSu54IILaGtr4/XXXx9fMWZmo1Czgb5jx+jaR2Pr1q3cfPPNbNmyhVNOOYX77rvv+Lzdu3dzxx138OSTT/LSSy+xYcMGHnvsMQDeeecd2traePnll/nIRz7C/fffP/5izMxKVLOBPm/e6NpHY+7cuVx00UUAXH/99fzgBz84Pm/Dhg0sXbqUhoYGpkyZQnt7O88++ywA06ZN48orrwRg8eLFdHd3j78YM7MS1Wygr1oF9fUD2+rrk/bxGnxpYamXGk6dOvX4spMnT+bIkSPjL8bMrEQ1G+jt7dDRAdOnJ9NNTcl0OU6I7tixg+effx6A73znO3z4wx8+Pm/JkiU888wzvPnmmxw9epQ1a9Zw8cUXj3+lZmbjVLOBDkl4t7XBxRdDd3f5rm45++yzWb16NQsWLGD//v187nOfOz5v1qxZ3HXXXVxyySVccMEFLF68mGXLlg3zamZm1aGISGXFra2t0XfVSJ8tW7awYMGCVOrp093dzZVXXsmrr75asXVk4ec0s9okaWNEtBabV9NH6GZmdoIDfZDm5uaKHp2bmVWKA93MLCcc6GZmOVH7gV5r4+eamVVI7Qe6mZkBtR7oFR4/t7u7m0WLFpX1Nc3MKqV2A32o8XPLGOoRwbFjx8r2emZmlVS7gV6h8XO7u7s5++yzueGGG1i0aBEHDx7kpptu4txzz+Wyyy7j4MGDALz00ku0tbVx/vnn8/GPf5z9+/ePa71mZuNVu4FewfFzf/KTn3DzzTezadMmdu7cyS233MKmTZs47bTT+N73vgfADTfcwN13380rr7zCeeedxze+8Y1xr9fMbDxqN9ArOH5uU1MTbW1tALS0tHDhhRcCJ4bEPXDgAG+99dbxQbk+/elPHx9C18wsLbUb6BUcP3fGjBnHn0/vG84RD4lrZtlWu4FeyfFzR3Dqqady+umn89xzzwHw0EMPeQhdM0tdzf6TaCAJ775/8/b001Vd9YMPPshnP/tZent7+cAHPsC3vvWtqq7fzGwwD5+bgonyc5pZ+Xn4XDOzCcCBbmaWE5kL9LS6gKol7z+fmaUnU4FeV1fHvn37cht6EcG+ffuoq6tLuxQzy6FMXeXS2NhIT08Pe/fuTbuUiqmrq6OxsTHtMswshzIV6FOnTqWlpSXtMszMalJJXS6SLpe0VdI2SSuKzJ8u6a8L838oqbnslZqZ2bBGDHRJk4HVwBXAQuA6SQsHLXYjsD8ifgH4I+DuchdqZmbDK+UIfQmwLSJei4jDwCPAskHLLAMeLDx/FLhUkspXppmZjaSUPvQ5wM5+0z3ALw61TEQckXQAOAN4s/9CkpYDywuTb0vaOpaigTMHv/YE5+0xkLfHCd4WA+VhezQNNaOqJ0UjogPoGO/rSOoa6tbXicjbYyBvjxO8LQbK+/YopctlFzC333Rjoa3oMpKmAKcC+8pRoJmZlaaUQN8AzJfUImkacC2wbtAy64BPF55/Angy8np3kJlZRo3Y5VLoE78VeAKYDPxVRGySdCfQFRHrgG8CD0naBvw7SehX0ri7bXLG22Mgb48TvC0GyvX2SG34XDMzK69MjeViZmZj50A3M8uJmgv0kYYhmCgkzZX0lKTNkjZJ+kLaNWWBpMmSXpT092nXkjZJp0l6VNKPJW2R9Etp15QWSV8s/J28KmmNpFwOeVpTgV7iMAQTxRHg9ohYCLQBt0zgbdHfF4AtaReREX8C/GNEnANcwATdLpLmAJ8HWiNiEcnFHZW+cCMVNRXolDYMwYQQEXsi4keF5z8j+WOdk25V6ZLUCPwa8EDataRN0qnAR0iuQCMiDkfEW6kWla4pwEmF+2Tqgd0p11MRtRboxYYhmNAhBlAY3fKDwA9TLiVtfwx8GTiWch1Z0ALsBb5V6IJ6QNKMtItKQ0TsAv4nsAPYAxyIiH9Kt6rKqLVAt0EkzQS+B/xuRPxH2vWkRdKVwBsRsTHtWjJiCvAh4M8j4oPAO8CEPOck6XSST/ItwGxghqTr062qMmot0EsZhmDCkDSVJMw7I2Jt2vWk7CLgakndJF1xvyLp4XRLSlUP0BMRfZ/aHiUJ+Inoo8C/RcTeiHgPWAv8cso1VUStBXopwxBMCIXhib8JbImIe9OuJ20R8XsR0RgRzSTviycjIpdHYaWIiJ8COyWdXWi6FNicYklp2gG0Saov/N1cSk5PEGfqX9CNZKhhCFIuKy0XAZ8C/q+klwptX4mI9emVZBnzO0Bn4eDnNeC3Uq4nFRHxQ0mPAj8iuTrsRXI6BIBv/Tczy4la63IxM7MhONDNzHLCgW5mlhMOdDOznHCgm5nlhAPdzCwnHOhmZjnx/wGapT1qOEzXIwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "(0.4980639522900839, 0.0026241899215518967) (0.8165568756531584, 0.011563011895523894)\n", "20\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAYBElEQVR4nO3de5Bc5Z3e8e+DbuNB3IKVLdBIM3IFdEFcbE0pk2AbsXgpcDBas9QW7GCcXYzKBtZOTGy0VrJZE5MCb8pxNiXWO2CzBAYRL5ZZOcHLVgISdgp5NTKXIMnCijySRhAYtIIYRiAEv/xxeqSe0Vx6Zrr7XPr5VE11n4v6/Pqo+5l33vP224oIzMws/05IuwAzM6sOB7qZWUE40M3MCsKBbmZWEA50M7OCcKCbmRWEA93MrCAc6GbjkPR7kvZIekvSo5L+wRj7XiBpq6SB0u0FZdu+IukFSb+W9CtJX6nLE7CG4UA3G4Okc4C/AD4D/AYwANw9yr4zgb8GHgROA+4H/rq0HkDA9aVtlwG3SLqmpk/AGooD3XJDUm+plft8qbX8XUm/IenHpVbv/5B0WpUP2wn8KCKeiog3gX8DXCXppBH2XQFMB74dEe9ExJ+RhPhvAkTENyPi5xFxJCJ2koT/hVWu1xqYA93y5neA3wLOBj4F/Bj4GjCH5PX8xZH+kaT5kl4f4+f3RjneOcBzgwsR8X+Aw6Xjj7Tv8zF0Po3nS+uH1yPgY8C2sZ+uWeWmp12A2QT954h4BUDST4BXI+KZ0vIPgUtG+kcRsRc4dRLHmw28MWzdG8BILfSJ7PsnJL+A7ptETWYjcqBb3rxSdv/QCMuzq3y8N4GTh607Gfj1ZPeVdAtJX/rHIuKdKtVp5i4XawylLpc3x/jpHOWfbgPOL3ucDwGzgBdH2fe8UnfKoPMo61aR9AfAauCSiOib6vMyK+cWujWEUpfLZFrv3cDTkj4G/By4HVgfESO10DcC7wFflPQd4MbS+icASr80/j1wcUTsnkQtZmNyC91sDBGxDfg8SbC/StIfftPg9tIIm6+V9j0M/DZJd8rrwB8Av11aD/AN4HRgS9lfBt+p01OxBiB/wYWZWTG4hW5mVhDjBrqk70l6VdILo2yXpD+TtKv0gY+PVL9MMzMbTyUt9L8k+ZjyaC4Hzir9rAL+fOplmZnZRI0b6BHxFPD3Y+yyEvgvkdgMnCrpjGoVaGZmlanGsMW5wL6y5b7SupeH7yhpFUkrnhNPPHHZokWLqnB4M7MK7NyZ3C5cmOs6tm7d+lpEzBlpW13HoUdEF9AF0N7eHj09PfU8vJk1shUrktuNG9OsYsp1SNoz2rZqjHLZD8wrW24prTMzszqqRqBvAK4vjXbpAN6IiOO6W8zMrLbG7XKRtI5knucPSuoD/i0wAyAivgM8BnwS2EUy+f/v16pYMzMb3biBHhHXjrM9gJurVpGZmU2KPylqZlYQDnQzs4JwoJuZFYQD3cysIBzoZmb10t0NmzfDpk3Q1pYsV5ED3cysHrq7YdUqeKf0NbJ79iTLVQx1B7qZWT2sWQMDA0PXDQwk66vEgW5mVg97905s/SQ40M3M6mH+/ImtnwQHuplZPdxxBzQ3D13X3JysrxIHuplZPXR2QlcXzJqVLLe2JsudnVU7RF3nQzcza2idnXDPPcn9GszL7ha6mRVfjcd/Z4UD3cyKrQ7jv7PCgW5mxVaH8d9Z4UA3s2Krw/jvrHCgm1mx1WH8d1Y40M2s2Oow/jsr8hXo3d3Q1ARSoa9Um1kV1WH8d1bkZxz6aFeqoZD/MWZWRTUe/50V+WmhN9CVajOzychPoDfQlWozs8nIT6A30JVqM7PJyE+gN9CVajOzychPoA9eqW5tTUa5FPhKtZnZZORnlAsk4e0ANzMbUX5a6GZmNiYHuplZQTjQzcwKwoFuZlYQDnQzs4JwoJuZFYQD3cysIBzok+FpfM0sgxzoE9VAXzhrZvlSUaBLukzSTkm7JK0eYft8SU9KekbS85I+Wf1SM8LT+JpZRo0b6JKmAWuBy4ElwLWSlgzb7V8D34+IDwPXAHdXu9DM8DS+ZpZRlbTQlwO7ImJ3RBwGHgZWDtsngJNL908BXqpeiRnjaXzNLKMqCfS5wL6y5b7SunJ/AlwnqQ94DPjDkR5I0ipJPZJ6+vv7J1FuBngaXzPLqGpdFL0W+MuIaAE+CTwg6bjHjoiuiGiPiPY5c+ZU6dB15ml8bTwrViQ/ZnVWyfS5+4F5ZcstpXXlbgAuA4iIpyU1AR8EXq1GkZnjaXzNLIMqaaFvAc6StEDSTJKLnhuG7bMXuARA0mKgCchpn4qZWT6NG+gRcQS4BXgc2EEymmWbpNslXVna7VbgRknPAeuAfx4RUauizczseBV9Y1FEPEZysbN83R+X3d8OXFjd0ix39uw59oGrtHzmM8ntAw+kV8OhQ8ntiy+mV4MdLyv/L4cOJdffaiBfX0Fn2bViRfJC/eEP061j2rTkdvbsxq4hS66+Orl95JF068jK/8u0afDeezV5aH/0P+88osLMShzoedbdDZs3w6ZN6U4SNljH3/0dLF8O69enU0cWrF8PW7fC00/7XFjducslr0abJAzqO6RyeB3798NXv5rcv+qq+tWRBevXJ8/98OFkuZHPRRal3eVTB26h51VWJgkbqY5Dh+DOO+tbB6TfOr7zzmMX3galdS6sITnQ8york4SNdryX6jydz2it43qG+mjPud7nwhqWAz2vsjJJ2GjHO/PM+taRhdbxaM+53ufCGpYDPa+yMknYSHV84AOw+rhp82srC63j1auT514ujXMx6Oqrjw0ZtIbgQM+rwUnCZs1KltOaJGx4HXPnwje/Wf+LgFloHV91VfLcZ85MltM6F9awPMolz7IySdhgHS++mN6HNlavTvrMy7td0mgdX3UVPPRQcr8BRlVYtjjQrRgGW8G33ppcGJ07Nwlzt46tgTjQrXpmzYI330zv+JdeChdckNwfnMsljXoGP9ad5rlwHdn13nuey8VyoLU17QqOXZQ8++zGrsF1ZNfwC+dV5IuiZlY7WZmeokE40M2KKAtBOtr0FI0e6hs3Jj814EA3q7YavmErkpUgzcr0FA1EaX2xUHt7e/T09KRybLNCa2tLQny41lbo7a1fHSecACPliwTvv1+/OgpG0taIaB9pm1voZkWTlXl+sjI9RQNxoJsVTVaCNCvTUzQQB7pZ0WQlSLMyPUUD8Th0s6IZDMwbbkgujLa2JmGeRpB2dsI99yT307xQ3CAc6GZF5CBtSO5yMTMrCAe6mVlBONDNzArCgW5mVhAOdDOzgvAoF7Oi8uiWhuMWuplZQTjQzcwKwoFuZlYQDnQzs4JwoJuZFYRHuZhZbXm0Td1U1EKXdJmknZJ2SVo9yj6/K2m7pG2SHqpumWZmNp5xW+iSpgFrgd8C+oAtkjZExPayfc4C/gi4MCIOSvqHtSrYzMxGVkkLfTmwKyJ2R8Rh4GFg5bB9bgTWRsRBgIh4tbplmpnZeCoJ9LnAvrLlvtK6cmcDZ0v6X5I2S7pspAeStEpSj6Se/v7+yVVsZmYjqtYol+nAWcAK4FrgHkmnDt8pIroioj0i2ufMmVOlQ5uZGVQW6PuBeWXLLaV15fqADRHxbkT8CniRJODNzKxOKgn0LcBZkhZImglcA2wYts+jJK1zJH2QpAtmd/XKNDOz8Ywb6BFxBLgFeBzYAXw/IrZJul3SlaXdHgcOSNoOPAl8JSIO1KpoMzM7niIilQO3t7dHT09PKsc2M8srSVsjon2kbf7ov5lZQTjQzcwKwoFuZlYQDnQzs4JwoJuZFYQD3cysIBzoZmYF4UA3MysIB7qZWUE40M3MCsKBbmZWEA50M7OCcKCbmRWEA93MrCAc6GZmBeFANzMrCAe6mVlBONDNzArCgW5mVhAOdDOzgnCgm5kVhAPdzKwgchXo3d3Q1gYnnJDcdnenXZGZWXZMT7uASnV3w6pVMDCQLO/ZkywDdHamV5eZWVbkpoW+Zs2xMB80MJCsNzOzHAX63r0TW29m1mhyE+jz509svZlZo8lNoN9xBzQ3D13X3JysNzOzHAV6Zyd0dUFrK0jJbVeXL4iamQ3KzSgXSMLbAW5mNrLctNDNzGxsDnQzs4JwoJuZFYQD3cysIBzoZmYFUVGgS7pM0k5JuyStHmO/35EUktqrV2L2dHdDU1MyfNKThFk5TyBnaRo30CVNA9YClwNLgGslLRlhv5OALwE/q3aRWTI4Sdg77yTLg5OE+Y2bvrR/0Q6+NvbsgQi/Nqz+KmmhLwd2RcTuiDgMPAysHGG/fwfcBbxdxfoyx5OEHS8LrdIs/KL1a8PSVkmgzwX2lS33ldYdJekjwLyI+O9jPZCkVZJ6JPX09/dPuNgs8CRhQ2WlVZqFMPVrw9I25Yuikk4AvgXcOt6+EdEVEe0R0T5nzpypHjoVniRsqCwEKWQjTP3asLRVEuj7gXllyy2ldYNOApYCGyX1Ah3AhqJeGPUkYUNlIUghG2GapddGFrrBsiTt6yt1ExFj/pDM97IbWADMBJ4Dzhlj/41A+3iPu2zZssirBx+MaG2NkJLbBx9Mu6L0tLZGJJ0tQ39aW+tbx4MPRjQ3D62hubn+/zdZeG1k5VxkRdHOB9ATo+Sqku1jk/RJ4NvANOB7EXGHpNtLD7xh2L4bgX8VET1jPWZ7e3v09Iy5i+XA8K8GhKQl9I1vwKc+Vd9afvQj+Na34OWX4Ywz4Mtfrn8NWXDxxfDSS8evP/NMePLJ+teTtiyej1mzkhljJ0PS1ogYsQekokCvBQd6cXR3J33me/c2dpBmxaJFSTt0OAl+8Yv615O2LJ6PWgV6rqbPtaHKg3T+/KSvNo3phT2tcbbMn5+MNhpp/dln17+etGXpfHR3ww03JMNrW1ur/571R/9zKivDBS17fHF2qKycj7p8VmK0zvVa/+T5omgWZOVipGWTL84eX0va56Na71mmelG0FtyHPjUnnDB6v+D779e/HrPh2tpG7upobYXe3npXk75qvWfH6kN3l0tOZWHctdlYsvIZhayox3vWgZ5TWekXNBuNGx1D1eM960DPqc5O6OpK/nyVktuuLo82sexwo2Ooerxn3YduZjWTlaG1ReJx6GaWCn9Gob7c5WJmVhAOdDOzgnCgm5kVhAPdzKwgHOhmZgXhQDczKwgHuplZQTjQzcwKwoFuZlYQDnQzs4JwoJuZFYQD3cysIBzoZmYF4UA3MysIB7qZWUE40M3MCsKBbmZWEA50M7OCcKCbmRWEA93MrCAc6GZmBeFANzMrCAe6mVlBONDNzArCgW5mVhAVBbqkyyTtlLRL0uoRtn9Z0nZJz0v6n5Jaq1+qmZmNZdxAlzQNWAtcDiwBrpW0ZNhuzwDtEXEe8AjwzWoXamZmY6ukhb4c2BURuyPiMPAwsLJ8h4h4MiIGSoubgZbqlmlmZuOpJNDnAvvKlvtK60ZzA/DjkTZIWiWpR1JPf39/5VWamdm4qnpRVNJ1QDvwpyNtj4iuiGiPiPY5c+ZU89BmZg1vegX77AfmlS23lNYNIekTwBrgooh4pzrlmZlZpSppoW8BzpK0QNJM4BpgQ/kOkj4M/AVwZUS8Wv0yzcxsPOMGekQcAW4BHgd2AN+PiG2Sbpd0ZWm3PwVmA38l6VlJG0Z5ODMzq5FKulyIiMeAx4at++Oy+5+ocl1mZjZBFQV6vbz77rv09fXx9ttvp11KzTQ1NdHS0sKMGTPSLsXMCiZTgd7X18dJJ51EW1sbktIup+oiggMHDtDX18eCBQvSLsfMCiZTc7m8/fbbnH766YUMcwBJnH766YX+C8TM0pOpQAcmHOYrViQ/eVHUX1Zmlr7MBbqZmU1OrgO9uxs2b4ZNm6CtLVmulc997nNs3769dgcwM5uiTF0UnYjubli1Ct4pfSZ1z55kGaCzs/rHu/fee6v/oGZmVZTbFvqaNTAwMHTdwECyfip6e3tZtGgRnZ2dLF68mKuvvpqBgQFWrFhBT08PAOvWrePcc89l6dKl3HbbbUf/7ezZs1mzZg3nn38+HR0dvPLKK1MrxsxsAnIb6Hv3Tmz9ROzcuZObbrqJHTt2cPLJJ3P33Xcf3fbSSy9x22238cQTT/Dss8+yZcsWHn30UQDeeustOjo6eO655/j4xz/OPffcM/VizMwqlNtAnz9/YusnYt68eVx44YUAXHfddfz0pz89um3Lli2sWLGCOXPmMH36dDo7O3nqqacAmDlzJldccQUAy5Yto7e3d+rFmJlVKLeBfscd0Nw8dF1zc7J+qoYPLax0qOGMGTOO7jtt2jSOHDky9WLMzCqU20Dv7ISuLpg1K1lubU2Wq3FBdO/evTz99NMAPPTQQ3z0ox89um358uVs2rSJ1157jffee49169Zx0UUXTf2gZmZTlNtAhyS8Ozrgoougt7d6o1sWLlzI2rVrWbx4MQcPHuQLX/jC0W1nnHEGd955JxdffDHnn38+y5YtY+XKlWM8mplZfSgiUjlwe3t7DI4aGbRjxw4WL16cSj2Dent7ueKKK3jhhRdqdowsPE8zyydJWyOifaRtuW6hm5nZMQ70Ydra2mraOjczqxUHuplZQTjQzcwKIv+Bnrf5c83MaiT/gW5mZkDeA73G8+f29vaydOnSqj6mmVmt5DfQR5s/t4qhHhG8//77VXs8M7Naym+g12j+3N7eXhYuXMj111/P0qVLOXToEDfeeCPnnHMOl156KYcOHQLg2WefpaOjg/POO49Pf/rTHDx4cErHNTObqvwGeg3nz/3lL3/JTTfdxLZt29i3bx8333wz27Zt49RTT+UHP/gBANdffz133XUXzz//POeeey5f//rXp3xcM7OpyG+g13D+3NbWVjo6OgBYsGABF1xwAXBsStw33niD119//eikXJ/97GePTqFrZpaW/AZ6DefPPfHEE4/enzU4nSOeEtfMsi2/gV7L+XPHccopp3Daaafxk5/8BIAHHnjAU+iaWepy+yXRQBLeg1/ztnFjXQ99//338/nPf56BgQE+9KEPcd9999X1+GZmw3n63BQ0yvM0s+rz9LlmZg3AgW5mVhCZC/S0uoDqpejPz8zSk6lAb2pq4sCBA4UNvYjgwIEDNDU1pV2KmRVQpka5tLS00NfXR39/f9ql1ExTUxMtLS1pl2FmBZSpQJ8xYwYLFixIuwwzs1yqqMtF0mWSdkraJWn1CNtnSfqvpe0/k9RW9UrNzGxM4wa6pGnAWuByYAlwraQlw3a7ATgYEf8I+I/AXdUu1MzMxlZJC305sCsidkfEYeBhYOWwfVYC95fuPwJcIknVK9PMzMZTSR/6XGBf2XIf8I9H2ycijkh6AzgdeK18J0mrgFWlxTcl7ZxM0cAHhz92g/P5GMrn4xifi6GKcD5aR9tQ14uiEdEFdE31cST1jPbR10bk8zGUz8cxPhdDFf18VNLlsh+YV7bcUlo34j6SpgOnAAeqUaCZmVWmkkDfApwlaYGkmcA1wIZh+2wAPlu6fzXwRBT100FmZhk1bpdLqU/8FuBxYBrwvYjYJul2oCciNgDfBR6QtAv4e5LQr6Upd9sUjM/HUD4fx/hcDFXo85Ha9LlmZlZdmZrLxczMJs+BbmZWELkL9PGmIWgUkuZJelLSdknbJH0p7ZqyQNI0Sc9I+m9p15I2SadKekTSLyTtkPRP0q4pLZL+Zel98oKkdZIKOeVprgK9wmkIGsUR4NaIWAJ0ADc38Lko9yVgR9pFZMR/Av4mIhYB59Og50XSXOCLQHtELCUZ3FHrgRupyFWgU9k0BA0hIl6OiJ+X7v+a5M06N92q0iWpBfhnwL1p15I2SacAHycZgUZEHI6I11MtKl3TgQ+UPifTDLyUcj01kbdAH2kagoYOMYDS7JYfBn6Wcilp+zbwVeD9lOvIggVAP3BfqQvqXkknpl1UGiJiP/AfgL3Ay8AbEfG36VZVG3kLdBtG0mzgB8C/iIj/l3Y9aZF0BfBqRGxNu5aMmA58BPjziPgw8BbQkNecJJ1G8pf8AuBM4ERJ16VbVW3kLdArmYagYUiaQRLm3RGxPu16UnYhcKWkXpKuuN+U9GC6JaWqD+iLiMG/2h4hCfhG9AngVxHRHxHvAuuBf5pyTTWRt0CvZBqChlCanvi7wI6I+Fba9aQtIv4oIloioo3kdfFERBSyFVaJiPi/wD5JC0urLgG2p1hSmvYCHZKaS++bSyjoBeJMfQXdeEabhiDlstJyIfAZ4H9Lera07msR8Vh6JVnG/CHQXWr87AZ+P+V6UhERP5P0CPBzktFhz1DQKQD80X8zs4LIW5eLmZmNwoFuZlYQDnQzs4JwoJuZFYQD3cysIBzoZmYF4UA3MyuI/w+jPzO2reuBhQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "(0.3791607121124395, 0.0029252620443503334) (0.7470131943371477, 0.027296044975130426)\n" ] } ], "source": [ "def jackknife_variance_estimator(f, x):\n", " N = len(x)\n", " X = sum(x) / N\n", " mean = f(X)\n", " return sum([(f((N*X - x[j])/(N-1)) - mean)**2 for j in range(N)])*(N-1)/N\n", "\n", "def bin_data(res, bin_size):\n", " assert len(res) % bin_size == 0\n", " n = len(res) // bin_size\n", " return [ sum(res[bin_size*i:bin_size*(i+1)])/bin_size for i in range(n) ]\n", "\n", "# C(t) = A*(exp(-t E) + exp(-(T-t)*E))\n", "def emp_fold(C):\n", " Cf = [(C[t] + C[(32 - t) % 32])/2. for t in range(17)]\n", " E = [np.arccosh((Cf[t].real + Cf[t+2].real) / 2 / Cf[t+1].real) for t in range(2,12)]\n", " return np.array(E)\n", "\n", "def plateau(C, t0, t1, E, emp):\n", " C = emp(C)\n", " return np.array([sum([C[t]/E[t] for t in range(t0,t1)]) / sum([1.0/E[t] for t in range(t0,t1)])])\n", "\n", "def mesons(data, nbin, tag):\n", " corrs = bin_data(np.array([data[x][\"pion\"] for x in data]),nbin)\n", " print(len(corrs))\n", " E_pi = emp_fold(np.mean(corrs,axis=0))\n", " E_pi_err = jackknife_variance_estimator(emp_fold, corrs)**0.5\n", " \n", " t0 = 4\n", " t1 = 10\n", " E_pi_plateau = plateau(np.mean(corrs,axis=0), t0, t1, E_pi_err, emp_fold)\n", " E_pi_plateau_err = jackknife_variance_estimator(lambda c: plateau(c, t0, t1, E_pi_err, emp_fold), corrs)**0.5\n", "\n", " corrs = bin_data(np.array([data[x][\"rho\"] for x in data]),nbin)\n", " E_rho = emp_fold(np.mean(corrs,axis=0))\n", " E_rho_err = jackknife_variance_estimator(emp_fold, corrs)**0.5\n", "\n", " t0 = 3\n", " t1 = 10\n", " E_rho_plateau = plateau(np.mean(corrs,axis=0), t0, t1, E_rho_err, emp_fold)\n", " E_rho_plateau_err = jackknife_variance_estimator(lambda c: plateau(c, t0, t1, E_rho_err, emp_fold), corrs)**0.5\n", "\n", " fig, ax = plt.subplots()\n", "\n", " plt.ylim(0.0,1.0)\n", " plt.title(tag)\n", "\n", " X=len(E_pi)\n", " ax.errorbar(np.arange(X),E_pi,E_pi_err,marker=\"o\",ls=\"\",color=\"blue\",label=\"pion\")\n", " ax.errorbar(np.arange(X),E_rho,E_rho_err,marker=\"o\",ls=\"\",color=\"red\",label=\"rho\")\n", " ax.fill_between(np.arange(t0,t1),[E_pi_plateau[0]-E_pi_plateau_err[0]]*(t1-t0),\n", " [E_pi_plateau[0]+E_pi_plateau_err[0]]*(t1-t0),alpha=0.1,color=\"blue\")\n", " ax.fill_between(np.arange(t0,t1),[E_rho_plateau[0]-E_rho_plateau_err[0]]*(t1-t0),\n", " [E_rho_plateau[0]+E_rho_plateau_err[0]]*(t1-t0),alpha=0.1,color=\"red\")\n", " plt.legend()\n", "\n", " plt.show()\n", "\n", " return ((E_pi_plateau[0], E_pi_plateau_err[0]),(E_rho_plateau[0], E_rho_plateau_err[0]))\n", "\n", "mpi_0p04, mrho_0p04 = mesons(mesons_0p04, 1, \"m = 0.04\")\n", "print(mpi_0p04, mrho_0p04)\n", "\n", "mpi_0p02, mrho_0p02 = mesons(mesons_0p02, 1, \"m = 0.02\")\n", "\n", "print(mpi_0p02, mrho_0p02)" ] }, { "cell_type": "markdown", "id": "c831be33", "metadata": {}, "source": [ "## $N_f=2$ baryon spectrum\n", "\n", "Next, we consider baryons, i.e., states made of three quarks with $N_f=2$. The most prominent example is the $I=1/2$ and $J^{P}=(1/2)^{+}$ which is the proton/neutron. Note that for degenerate up and down quarks the proton and neutron have the same mass.\n", "\n", "One can show that a baryonic operator transforming in the $J=1/2$, $I=1/2$ representation is given by\n", "$$\n", "\\begin{align}\n", "O_N &= \\varepsilon_{abc} u_a (u_b^T C \\gamma_5 d_c) \\,,\n", "\\end{align}\n", "$$\n", "where $C$ is the charge conjugation matrix of chapter 11 and $\\varepsilon_{abc}$ is the total antisymmetric tensor of rank 3.\n", "\n", "Since we want a state that is a parity eigenstate with $P=1$, we also consider the parity transformation (see chapter 11) of this operator and find\n", "$$\n", "\\begin{align}\n", "O_N^P(\\vec{x},t) &= \\gamma_4 O_N(-\\vec{x},t) \\,.\n", "\\end{align}\n", "$$\n", "If we also project to zero momentum by summing the operators over all spatial positions, we find that a parity\n", "definite operator is given by\n", "$$\n", "\\begin{align}\n", "O_N^+(t) = \\sum_{\\vec{x}} P_+ O_N(\\vec{x},t)\n", "\\end{align}\n", "$$\n", "with\n", "$$\n", "\\begin{align}\n", "P_+ &= \\frac12 (1 + \\gamma_4) \\,.\n", "\\end{align}\n", "$$\n", "\n", "The two-point correlator for fixed gauge field is then given by\n", "$$\n", "\\begin{align}\n", "C(t) &= \\langle O_N^+(t) (O_N^+(0))^\\dagger \\rangle_\\psi \\\\\n", "&=\\varepsilon_{abc} \\varepsilon_{a'b'c'} (Cg5)_{\\alpha'\\beta'} (Cg5)_{\\alpha\\beta} (P_+)_{\\gamma\\gamma'} \n", " \\langle d \\bar d \\rangle_{\\beta'\\beta;b'b} ( \\langle u \\bar u \\rangle_{\\alpha'\\alpha;a'a} \\langle u \\bar u \\rangle_{\\gamma'\\gamma;c'c} - \\langle u \\bar u \\rangle_{\\alpha'\\gamma;a'c} \\langle u \\bar u \\rangle_{\\gamma'\\alpha;c'a}) \\,,\n", "\\end{align}\n", "$$\n", "where the space-time argument of all propagators is $\\vec{x},t;0$ and a sum over $\\vec{x}$ is implied.\n", "\n", "Note that there are also baryons with negative parity and with higher $J$ and $I$ quantum numbers, e.g.,\n", "\n", "| I| $J^{P}$ | Name | Mass |\n", "| --- | --- | --- | --- |\n", "| 1/2 | $1/2^{+}$ | proton | 938.2720813(58) MeV |\n", "| 1/2 | $1/2^{+}$ | neutron | 939.5654133(58) MeV |\n", "| 1/2 | $1/2^{-}$ | N(1535) | $\\approx 1535$ MeV |\n", "| 3/2 | $3/2^{+}$ | $\\Delta$ | 1232(2) MeV |\n", "\n", "We conclude this chapter by discussing code to compute the baryon correlation function." ] }, { "cell_type": "code", "execution_count": 135, "id": "e95cad3a", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "GPT : 61722.610912 s : Initializing gpt.random(d,vectorized_ranlux24_389_64) took 0.000323057 s\n", "4.073750915814535e-31\n" ] } ], "source": [ "def proton(Q1, Q2):\n", " #\n", " # eps(a,b,c) eps(a',b',c') (Cg5)_alpha'beta' (Cg5)_alphabeta Pp_gammagamma' D^-1_beta'beta_b'b\n", " # (U^-1_alpha'alpha_a'a U^-1_gamma'gamma_c'c - U^-1_alpha'gamma_a'c U^-1_gamma'alpha_c'a)\n", " #\n", " # =\n", " #\n", " # eps(a,b,c) eps(a',b',c') (Cg5)_alpha'beta' (Cg5)_alphabeta Pp_gammagamma' D^-1_beta'beta_b'b\n", " # U^-1_alpha'alpha_a'a U^-1_gamma'gamma_c'c\n", " # -\n", " # eps(a,b,c) eps(a',b',c') (Cg5)_alpha'beta' (Cg5)_alphabeta Pp_gammagamma' D^-1_beta'beta_b'b\n", " # U^-1_alpha'gamma_a'c U^-1_gamma'alpha_c'a\n", " #\n", " # =\n", " #\n", " # eps(c,a,b) eps(c',a',b') (U^-1 C g5)_alpha'beta_a'a (C g5 D^-1)^Tspin_betaalpha'_b'b Tr_S[Pp U^-1_c'c]\n", " # +\n", " # eps(a,b,c) eps(a',b',c') Pp_gammagamma' \n", " # (U^-1 Cg5)_gamma'beta_c'c (Cg5 D^-1)^Tspin_betaalpha'_b'b U^-1_alpha'gamma_a'a\n", " #\n", " # =\n", " #\n", " # Tr_S[diquark(U^-1 C g5, C g5 D^-1)_{cc'}] Tr_S[Pp U^-1_c'c]\n", " # +\n", " # diquark(U^-1 C g5, C g5 D^-1)_{gamma'alpha', aa'} Pp_gammagamma' \n", " # U^-1_alpha'gamma_a'a\n", " #\n", " # =\n", " #\n", " # Tr_C[Tr_S[diquark(U^-1 C g5, C g5 D^-1)] Tr_S[U^-1 Pp]]\n", " # +\n", " # Tr[ U^-1 Pp diquark(U^-1 C g5, C g5 D^-1) ]\n", " # \n", " #\n", " # with\n", " #\n", " # diquark(Q1,Q2)_{gammagamma', cc'} = eps(c,a,b) eps(c',a',b') (Q1)_gammabeta_a'a (Q2)^Tspin_betagamma'_b'b \n", " #\n", " C = 1j*g.gamma[1].tensor()*g.gamma[3].tensor()\n", " Gamma = C*g.gamma[5].tensor()\n", " \n", " Pp = (g.gamma[\"I\"].tensor() + g.gamma[3].tensor())*0.5\n", " \n", " dq = g.qcd.baryon.diquark(g(Q1*Gamma),g(Gamma*Q2))\n", " q1p = g(Q1 * Pp)\n", " \n", " return g(\n", " g.color_trace(g.spin_trace(dq) * g.spin_trace(q1p) + g.spin_trace(dq*q1p))\n", " )\n", "\n", "rng = g.random(\"d\")\n", "U = rng.cnormal(g.mspincolor(grid))\n", "D = rng.cnormal(g.mspincolor(grid))\n", "proton1 = proton(U, D)\n", "\n", "V = rng.element(g.mcolor(grid))\n", "\n", "# test gauge invariance\n", "proton2 = proton(g(V * U * g.adj(V)),g(V * D * g.adj(V)))\n", "\n", "print(g.norm2(proton1 - proton2)/g.norm2(proton1))" ] }, { "cell_type": "code", "execution_count": 188, "id": "489be7a9", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "98\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEICAYAAABRSj9aAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAa+klEQVR4nO3dfZRU9Z3n8ffH5kkaE4n2ushDN7PjKkoUoQ6aUSPGkWDGiDLmBNIm6mp619U8TM4xEp0dz5JjNHGPcZM1YxhlfRjEGQkYskej7hpUdiWhQWMUoiGI2B0NLRgDAvLgd/+4t7Boqruru6u7ituf1zl1qu7vPtS3Svz0rd/93XsVEZiZWXYdVukCzMysbznozcwyzkFvZpZxDnozs4xz0JuZZZyD3sws4xz0ZmYZ56A36yFJX5D0uqT3JD0i6WOdLDtJ0mpJO9LnSUWWGSJpnaSWPi3cBhwHvVkPSDoJ+DHwReAYYAfwow6WHQL8FPhnYCRwH/DTtL3QdUBbX9VsA5eD3g55kjZKuk7Si+ne9T2SjpH0mKRtkv63pJFlfttG4GcR8UxEbAf+CzBL0hFFlp0GDALuiIj3I+IHgIBPFXyG8cClwC1lrtPMQW+Z8bfAecC/Bz4LPAbcANSR/Dv/arGVJI2T9KdOHl/o4P1OAn6dn4iI3wO70/cvtuyLceD1Rl5M2/N+mNa7s4TPatYtgypdgFmZ/DAi/ggg6Vlgc0Q8n04vBc4ttlJEbAKO7MH7jQDebdf2LlBsj77TZSVdDNRExFJJ03pQi1mnHPSWFX8seL2zyPSIMr/fduAj7do+AmzrzrKSaoHvAZ8pc31m+7nrxga0tOtmeyePxg5WfRk4pWA7fwEMBV7tYNmTJamg7eS0/TigAXhW0lvAEmCUpLckNfT6A5rhPXob4NKum57s7S8EnpN0FrAGmAcsiYhie/TLgX3AVyXdBXw5bX8K+AAYW7DsXwH/A5iMR+BYmXiP3qwHIuJl4D+RBP5mkv72/5yfn474uSFddjdwEfAl4E/AfwAuiojdEbE3It7KP4CtwAfp9L7+/EyWXfKNR8zMss179GZmGddl0EsaK+kXktZKelnS14osI0k/kLQ+PWllcsG8yyT9Ln1cVu4PYGZmneuy60bSKGBURKxJz/pbTdK/uLZgmc8AXyEZInYa8N8j4rT02h/NQA6IdN0pEfFOn3waMzM7SJd79BHxZkSsSV9vA9YBo9stNhO4PxIrgSPTPxCfBp6MiK1puD8JzCjrJzAzs051a3hlOq73VOCX7WaNBt4omG5J2zpqL7btJqAJoLa2dsoJJ5zQndLMzAa01atXvx0RdcXmlRz0kkYAPwG+HhF/LldxeRExH5gPkMvlorm5udxvYWaWWZJe72heSaNuJA0mCfmFEbGkyCKtHHjSx5i0raN2MzPrJ6WMuhFwD7AuIm7vYLFlwJfS0TenA+9GxJvA48B0SSPTy8ROT9vMzKyflNJ1cwbJzRV+I+mFtO0GYBxARNwFPEoy4mY9yQ0YrkjnbZX0bWBVut68iNhaturNzKxLXQZ9RKwguUlCZ8sEcE0H8xYAC3pUnZkZsGfPHlpaWti1a1elS6m4YcOGMWbMGAYPHlzyOr6omZlVvZaWFo444ggaGho48CKgA0tEsGXLFlpaWhg/fnzJ6/kSCGZW9Xbt2sVRRx01oEMeQBJHHXVUt3/ZOOjN7JAw0EM+ryffg4PezDJp2rTkYQ56M7N+8Z3vfKdi7+2gN7PMWbgQVq6Ep5+GhoZkuj/s29fxvWIc9GZmZbJwITQ1wfvvJ9Ovv55M9zbsN27cyAknnEBjYyMTJkzgkksuYceOHTQ0NHD99dczefJkHn74YRYtWsTHP/5xJk6cyPXXXw/A3Llz2blzJ5MmTaKxMbkN8e23387EiROZOHEid9xxx/73mDBhAl/+8pc56aSTmD59Ojt37uxd4ZAM16m2x5QpU8LMLG/t2rUlL1tfHwEHP+rre1fDa6+9FkCsWLEiIiKuuOKKuO2226K+vj6++93vRkREa2trjB07NjZv3hx79uyJc845J5YuXRoREbW1tfu31dzcHBMnTozt27fHtm3b4sQTT4w1a9bEa6+9FjU1NfH8889HRMTnPve5eOCBBw6qpdj3ATRHB5nqPXozy5RNm7rX3h1jx47ljDPOAODSSy9lxYoVAHz+858HYNWqVUybNo26ujoGDRpEY2MjzzzzzEHbWbFiBRdffDG1tbWMGDGCWbNm8eyzzwIwfvx4Jk2aBMCUKVPYuHFjr+t20JtZpowb17327mg/tDE/XVtb2/uNp4YOHbr/dU1NDXv37u31Nh30ZpYpN98Mw4cf2DZ8eNLeW5s2beK5554D4MEHH+TMM888YP7UqVN5+umnefvtt9m3bx+LFi3i7LPPBmDw4MHs2bMHgLPOOotHHnmEHTt28N5777F06VLOOuus3hfYAQe9mWVKYyPMnw/5HeP6+mQ6PQbaK8cffzx33nknEyZM4J133uHqq68+YP6oUaO49dZbOeecczjllFOYMmUKM2fOBKCpqYmTTz6ZxsZGJk+ezOWXX87UqVM57bTTuOqqqzj11FN7X2AHurxnbCX4xiNmVmjdunVMmDChW+vkT5Zavrw8NWzcuJELLriAl156qTwb7IVi34ek1RGRK7a8L2pmZplUroDPAnfdmJmVoKGhoSr25nvCQW9mh4Rq7GauhJ58Dw56M6t6w4YNY8uWLQM+7CO9Hv2wYcO6tZ776M2s6o0ZM4aWlhba2toqXUrF5e8w1R1dBr2kBcAFwOaImFhk/nVAfuDSIGACUBfJ/WI3AtuAfcDejo4Im5l1ZvDgwd26o5IdqJSum3uBGR3NjIjbImJSREwCvgU8HQfeAPycdL5D3sysAroM+oh4Btja1XKpOcCiXlVkZmZlVbaDsZKGk+z5/6SgOYAnJK2W1FSu9zIzs9KV82DsZ4H/267b5syIaJX0b4AnJf02/YVwkPQPQRPAuHJcfcjMzIDyDq+cTbtum4hoTZ83A0uBqR2tHBHzIyIXEbm6uroylmVmNrCVJeglfRQ4G/hpQVutpCPyr4HpwKF5WpmZ2SGslOGVi4BpwNGSWoCbgMEAEXFXutjFwBMR8V7BqscAS9PrNQ8CHoyIn5evdDMzK0WXQR8Rc0pY5l6SYZiFbRuAU3pamJmZlYcvgWBmlnEOejOzjHPQm5llnIPezCzjHPRmZhnnoDczyzgHvZlZxjnozcwyzkFvZpZxDnozs4xz0JuZZZyD3sws4xz0ZmYZ56A3M8s4B72ZWcY56M3MMs5Bb2aWcQ56M7OMc9CbmWVcl0EvaYGkzZJe6mD+NEnvSnohffxDwbwZkl6RtF7S3HIWbmZmpSllj/5eYEYXyzwbEZPSxzwASTXAncD5wInAHEkn9qZYMzPrvi6DPiKeAbb2YNtTgfURsSEidgMPATN7sB0zM+uFcvXRf0LSryU9JumktG008EbBMi1pW1GSmiQ1S2pua2srU1lmZlaOoF8D1EfEKcAPgUd6spGImB8RuYjI1dXVlaEsMzODMgR9RPw5Iranrx8FBks6GmgFxhYsOiZtMzOzftTroJf0byUpfT013eYWYBVwnKTxkoYAs4FlvX0/MzPrnkFdLSBpETANOFpSC3ATMBggIu4CLgGulrQX2AnMjogA9kq6FngcqAEWRMTLffIpzMysQ0oyubrkcrlobm6udBlmZocMSasjIldsns+MNTPLOAe9mVnGOejNzDLOQW9mlnEOejOzjHPQm5llnIPezCzjMhP0CxdCQwMcdljyvHBhpSsyM6sOXZ4ZeyhYuBCammDHjmT69deTaYDGxsrVlTdtWvK8fHklqzCzataXOZGJPfobb/ww5PN27EjazcwGukwE/aZN3WvvTwsXwsqV8PTT7lIys8rIRNCPG9e99v6S71J6//1kOt+l5LA3s/6UiaC/+WYYPvzAtuHDk/ZKqvYupWnTPuwXNLPsykTQNzbC/PlQXw9S8jx/fuUPxFZzl5KZDRyZGHUDSahXOtjbGzcu6a4p1m5m1l8ysUdfraq1Swl8kNhsIHHQ96F8l9LQocl0tXQp+SCx2cDiO0z1g2o7YaqhoXiXUn09bNzY39UcrNq+L7O+tnAhXHllsvNVX5/86u/uDmFnd5gq5Z6xC4ALgM0RMbHI/EbgekDANuDqiPh1Om9j2rYP2NtREVlXbYHlg8Rm1aOjX9hQvl//pXTd3AvM6GT+a8DZEfFx4NvA/Hbzz4mISQM15KtRtZ53YDYQ9ccw7C6DPiKeAbZ2Mv//RcQ76eRKYEyZarM+Us0Hic0Gmv74hV3ug7FXAo8VTAfwhKTVkpo6W1FSk6RmSc1tbW1lLssKVetBYvBoIBt4+uMXdtmCXtI5JEF/fUHzmRExGTgfuEbSJztaPyLmR0QuInJ1dXXlKss60NgIp58OZ5+dHICtlpD3aCAbaPrjF3ZZgl7SycDdwMyI2JJvj4jW9HkzsBSYWo73s2yq9ktGVCtfyuLQ1h+/sHt9ZqykccAS4IsR8WpBey1wWERsS19PB+b19v2sfDwayKw6NDbCP/1T8rov/r8sZXjlImAacLSkFuAmYDBARNwF/ANwFPAjSfDhMMpjgKVp2yDgwYj4efk/gmWFLxmRLT4fonp0GfQRMaeL+VcBVxVp3wCc0vPSbKC5+eYD7xQGHg1kVg6+BIJVjWoeDVStPErJSpGZq1daNuT7KnftgvvvT9pefbXzdQaqn/0M/v7vDxyldNVV8NZb8NnPVrY2gJ07k2f/9yvNzp3JZdb7goPeqlIEjBhR6Sqq2/e/n/xBLLRrV9I+p9MO1763ZAm88ALs3g2f+hTMnQuzZlW2prxLLkmeFy+ubB3t1dTAvn19s20HvVWd5cu9F1iKP/yhe+39ZckS+OY3k5AHaG1NpqF6wn6gcR+92SHq2GO7195fbr31w26bvJ07k3arDAe92SFq7lw4/PAD2w4/PGmvpGr9pVHtFi+GBx7om2076M0OUbNmwfe+B0OGJNOjRyfTle4eqdZfGpB0K61eDc89B1OnJtMDgYPerARLlsD48UmYVlNAzJoFU6bAJz4Bv/pV5UMeqveXRkfHDqrlv2Vf8sFYq0pDh8L27ZWuIpEfxlgYENddlwxrrIZhjPmRGtXyfU2fnnw3N9yQfGfHHgvf+EbSXskab7ml+LGDW25JaqsG+XNIys23EjTrgm+92DPVVtdhhyXDdtuT4IMP+r+ecuvVrQTNBrpqv9hatQRptRvI11JyH71ZF3zrxWwYyHdWc9CbdWEgB0RvLF9eXb82BvK1lNx1Y9aFfBDceGPSXTNuXBLyAyEgsqavr/terRz0ZiVobHSw26HLQW9mA8pA2pPPcx+9mVnGOejNzDKupKCXtEDSZkkvdTBfkn4gab2kFyVNLph3maTfpY/LylW4mZmVptQ9+nuBGZ3MPx84Ln00Af8IIOljJDcTPw2YCtwkaWRPizUzs+4rKegj4hlgayeLzATuj8RK4EhJo4BPA09GxNaIeAd4ks7/YJiZWZmVq49+NPBGwXRL2tZR+0EkNUlqltTc1tZWprLMzKxqDsZGxPyIyEVErq6urtLlmJllRrmCvhUYWzA9Jm3rqN3MzPpJuYJ+GfCldPTN6cC7EfEm8DgwXdLI9CDs9LTNzMz6SUlnxkpaBEwDjpbUQjKSZjBARNwFPAp8BlgP7ACuSOdtlfRtYFW6qXkR0dlBXTMzK7OSgj4i5nQxP4BrOpi3AFjQ/dLMzKwcquZgrJmZ9Q0HvZlZxjnozcwyzkFvZpZxDnozs4xz0JuZZZyD3sws4xz0ZmYZ56A3M8s4B72ZWcY56M3MMs5Bb2aWcQ56M7OMc9CbmWWcg97MLOMc9GZmGeegNzPLOAe9mVnGlRT0kmZIekXSeklzi8z/vqQX0serkv5UMG9fwbxlZazdzMxK0OU9YyXVAHcC5wEtwCpJyyJibX6ZiPi7guW/ApxasImdETGpbBWbmVm3lLJHPxVYHxEbImI38BAws5Pl5wCLylGcmZn1XilBPxp4o2C6JW07iKR6YDzwVEHzMEnNklZKuqijN5HUlC7X3NbWVkJZZmZWinIfjJ0NLI6IfQVt9RGRA74A3CHp3xVbMSLmR0QuInJ1dXVlLsvMbOAqJehbgbEF02PStmJm067bJiJa0+cNwHIO7L83M7M+VkrQrwKOkzRe0hCSMD9o9IykE4CRwHMFbSMlDU1fHw2cAaxtv66ZmfWdLkfdRMReSdcCjwM1wIKIeFnSPKA5IvKhPxt4KCKiYPUJwI8lfUDyR+XWwtE6ZmbW93RgLleHXC4Xzc3NlS7DzOyQIWl1ejz0ID4z1sws4xz0ZmYZ56A3M8s4B72ZWcY56M3MMs5Bb2aWcQ56M7OMc9CbmWWcg97MLOMc9GZmGeegNzPLOAe9mVnGOejNzDLOQW9mlnEOejOzjHPQm5llnIPezCzjHPRmZhlXUtBLmiHpFUnrJc0tMv9ySW2SXkgfVxXMu0zS79LHZeUs3szMutblzcEl1QB3AucBLcAqScuK3OT7XyLi2nbrfgy4CcgBAaxO132nLNWbmVmXStmjnwqsj4gNEbEbeAiYWeL2Pw08GRFb03B/EpjRs1LNzKwnSgn60cAbBdMtaVt7fyvpRUmLJY3t5rpIapLULKm5ra2thLLMzKwU5ToY+zOgISJOJtlrv6+7G4iI+RGRi4hcXV1dmcoyM7NSgr4VGFswPSZt2y8itkTE++nk3cCUUtc1M7O+VUrQrwKOkzRe0hBgNrCscAFJowomLwTWpa8fB6ZLGilpJDA9bTMzs37S5aibiNgr6VqSgK4BFkTEy5LmAc0RsQz4qqQLgb3AVuDydN2tkr5N8scCYF5EbO2Dz2FmZh1QRFS6hoPkcrlobm6udBlmZocMSasjIldsns+MNTPLOAe9mVnGOejNzDLOQW9mlnEOejOzjHPQm5llnIPezCzjHPRmZhnnoDczyzgHvZlZxjnozcwyzkFvZpZxDnozs4xz0JuZZZyD3sws4xz0ZmYZ56A3M8s4B72ZWcaVFPSSZkh6RdJ6SXOLzP+GpLWSXpT0fyTVF8zbJ+mF9LGs/bpmZta3urw5uKQa4E7gPKAFWCVpWUSsLVjseSAXETskXQ18D/h8Om9nREwqb9lmZlaqUvbopwLrI2JDROwGHgJmFi4QEb+IiB3p5EpgTHnLNDOzniol6EcDbxRMt6RtHbkSeKxgepikZkkrJV3U/RLNzKw3uuy66Q5JlwI54OyC5vqIaJX0F8BTkn4TEb8vsm4T0AQwbty4cpZlZjaglbJH3wqMLZgek7YdQNJfAzcCF0bE+/n2iGhNnzcAy4FTi71JRMyPiFxE5Orq6kr+AGZm1rlSgn4VcJyk8ZKGALOBA0bPSDoV+DFJyG8uaB8paWj6+mjgDKDwIK6ZmfWxLrtuImKvpGuBx4EaYEFEvCxpHtAcEcuA24ARwMOSADZFxIXABODHkj4g+aNya7vROmZm1scUEZWu4SC5XC6am5srXYaZ2SFD0uqIyBWb5zNjzcwyzkFvZpZxDnozs4xz0JuZZZyD3sws4xz0ZmYZ56A3M8s4B72ZWcY56M3MMs5Bb2aWcQ56M7OMc9CbmWWcg97MLOMc9GZmGeegNzPLOAe9mVnGOejNzDLOQW9mlnEOejOzjCsp6CXNkPSKpPWS5haZP1TSv6TzfympoWDet9L2VyR9uoy1m5lZCboMekk1wJ3A+cCJwBxJJ7Zb7ErgnYj4S+D7wHfTdU8EZgMnATOAH6XbMzOzflLKHv1UYH1EbIiI3cBDwMx2y8wE7ktfLwbOlaS0/aGIeD8iXgPWp9szM7N+MqiEZUYDbxRMtwCndbRMROyV9C5wVNq+st26o4u9iaQmoCmd3C7plRJqK+Zo4O0ertuXXFf3uK7ucV3dk8W66juaUUrQ94uImA/M7+12JDVHRK4MJZWV6+oe19U9rqt7BlpdpXTdtAJjC6bHpG1Fl5E0CPgosKXEdc3MrA+VEvSrgOMkjZc0hOTg6rJ2yywDLktfXwI8FRGRts9OR+WMB44DflWe0s3MrBRddt2kfe7XAo8DNcCCiHhZ0jygOSKWAfcAD0haD2wl+WNAuty/AmuBvcA1EbGvjz5LXq+7f/qI6+oe19U9rqt7BlRdSna8zcwsq3xmrJlZxjnozcwyLjNB39VlGipF0gJJmyW9VOla8iSNlfQLSWslvSzpa5WuKU/SMEm/kvTrtLb/Wuma8iTVSHpe0v+qdC2FJG2U9BtJL0hqrnQ9eZKOlLRY0m8lrZP0iSqo6fj0e8o//izp65WuC0DS36X/5l+StEjSsLJtOwt99OllFV4FziM5KWsVMCci1la0MEDSJ4HtwP0RMbHS9QBIGgWMiog1ko4AVgMXVcn3JaA2IrZLGgysAL4WESu7WLXPSfoGkAM+EhEXVLqePEkbgVxEVNUJQJLuA56NiLvTEXvDI+JPFS5rvzQ3WoHTIuL1CtcymuTf+okRsTMdxPJoRNxbju1nZY++lMs0VEREPEMyEqlqRMSbEbEmfb0NWEcHZyz3t0hsTycHp4+K741IGgP8DXB3pWs5FEj6KPBJkhF5RMTuagr51LnA7ysd8gUGAYen5yINB/5Qrg1nJeiLXaahKoKr2qVXGj0V+GWFS9kv7SJ5AdgMPBkR1VDbHcA3gQ8qXEcxATwhaXV6KZFqMB5oA/5n2t11t6TaShfVzmxgUaWLAIiIVuC/AZuAN4F3I+KJcm0/K0FvPSBpBPAT4OsR8edK15MXEfsiYhLJmdRTJVW0y0vSBcDmiFhdyTo6cWZETCa5wuw1aXdhpQ0CJgP/GBGnAu8B1XTsbAhwIfBwpWsBkDSSpBdiPHAsUCvp0nJtPytB70stdFPa//0TYGFELKl0PcWkP/V/QXKJ60o6A7gw7Qt/CPiUpH+ubEkfSvcGiYjNwFKq4wqxLUBLwa+xxSTBXy3OB9ZExB8rXUjqr4HXIqItIvYAS4C/KtfGsxL0pVymwVLpAc97gHURcXul6ykkqU7Skenrw0kOsP+2kjVFxLciYkxENJD823oqIsq2t9UbkmrTA+qkXSPTgYqP8IqIt4A3JB2fNp1LcoZ8tZhDlXTbpDYBp0sanv7/eS7JsbOyqJqrV/ZGR5dpqHBZAEhaBEwDjpbUAtwUEfdUtirOAL4I/CbtCwe4ISIerVxJ+40C7ktHRBwG/GtEVNVwxipzDLA0yQYGAQ9GxM8rW9J+XwEWpjtfG4ArKlwPsP8P4nnAf6x0LXkR8UtJi4E1JJeLeZ4yXg4hE8MrzcysY1npujEzsw446M3MMs5Bb2aWcQ56M7OMc9CbmWWcg97MLOMc9GZmGff/AbRANO9aKXsyAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "40\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEICAYAAABRSj9aAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAbJklEQVR4nO3dfZRU9Z3n8fdnGpQAJhLtNcjz7BoFiSLUoLNqhHEkmHXEOOYE0ibqaHrjap5PIjGz4y5Zc4zucdzMmDGMMhoHcUYEw+xq1KxBw45EGjQqGAxBxO6Y0IIxIigP+e4f95ZTNFXd1d3VXZfbn9c5daru7z7Ut1r83Fu/+6t7FRGYmVl+/UG9CzAzs77loDczyzkHvZlZzjnozcxyzkFvZpZzDnozs5xz0JuZ5ZyD3qyHJH1S0suS3pL0gKT3d7LsFElrJe1Kn6eUzPuqpOclvSnpJUlf7ZcPYAOGg96sBySdCHwP+BRwDLAL+G6FZQ8DfgD8IzACuAv4QdoOIODT6bzZwNWS5vbpB7ABxUFvhzxJW9Kj4mfTo+s7JB0j6aH0KPlHkkbU+G2bgH+JiCciYifwX4ELJR1RZtkZwCDgloh4JyK+QxLufwIQETdGxLqI2BcRG0l2CqfXuF4bwBz0lhd/DpwDfBD4M+Ah4FqgkeTf+efLrSRprKTfdvL4ZIX3OxH4WXEiIn4J7Enfv9yyz8aB1xt5Nm3vWI+AM4H1nX9cs+oNqncBZjXyNxHxGwBJPwG2RcTT6fRy4OxyK0XEVuDIHrzfcOCNDm1vAOWO6Luz7H8j2TH9Qw9qMivLQW958ZuS17vLTA+v8fvtBN7boe29wJs9XVbS1SR99WdGxDs1qtPMXTc2sKVdNzs7eTRVWHU9cHLJdv4QOBx4scKyJ6XdMkUnUdI9I+kvgPnA2RHR2tvPZVbKR/Q2oKVdNz052l8MPCnpTGAdsABYFhHljuhXAvuBz0u6DfhM2v4YQLoz+RYwMyI296AWs075iN6sByJiPfBZksDfRtLf/l+K89MRP9emy+4BLiDplvkt8BfABWk7wP8AjgLWlHyTuK2fPooNAPKNR8zM8s1H9GZmOddl0EsaI+nHkjZIWi/pC2WWkaTvSNqU/mhlasm8SyT9In1cUusPYGZmneuy60bSSGBkRKxLf/W3lqR/cUPJMh8FPgd8FDgV+F8RcWp67Y8WoABEuu60iHi9Tz6NmZkdpMsj+oh4NSLWpa/fBF4ARnVYbA7w/UisBo5MdxAfAR6NiB1puD9Kci0PMzPrJ90aXilpPHAK8NMOs0YBr5RMt6ZtldrLbbsZaAYYNmzYtBNOOKE7pZmZDWhr1659LSIay82rOuglDQfuB74YEb+rVXFFEbEQWAhQKBSipaWl1m9hZpZbkl6uNK+qUTeSBpOE/OKIWFZmkTZgTMn06LStUruZmfWTakbdCLgDeCEibq6w2Arg0+nom9OANyLiVeBhYJakEellYmelbWZm1k+q6bo5neTmCs9JeiZtuxYYCxARtwEPkoy42URyA4bL0nk7JH0TWJOutyAidtSsejMz61KXQR8Rq0huktDZMgFcVWHeImBRj6ozMwP27t1La2srb7/9dr1LqbshQ4YwevRoBg8eXPU6vqiZmWVea2srRxxxBOPHj+fAi4Dmx8aNyfPxx1deJiLYvn07ra2tTJgwoept+xIIZpZ5b7/9NkcddVRuQ75akjjqqKO6/c3GQW9mh4SBHvJFPfk7OOjNLJdmzEge5qA3M+sX3/rWt+r23g56M8udxYth9Wp4/HEYPz6Z7g/79++vOM9Bb2ZWI4sXQ3MzvJPeXv3ll5Pp3ob9li1bOOGEE2hqamLixIlcdNFF7Nq1i/Hjx3PNNdcwdepU7rvvPpYsWcKHPvQhJk+ezDXXXAPA/Pnz2b17N1OmTKGpKbkN8c0338zkyZOZPHkyt9xyCwCtrVuYOHEin/nMZzjxxBOZNWsWu3fv7l3hkAzXydpj2rRpYWZWtGHDhqqXHTcuAg5+jBvXuxpeeumlAGLVqlUREXHZZZfFTTfdFOPGjYtvf/vbERHR1tYWY8aMiW3btsXevXtj5syZsXz58oiIGDZs2LvbamlpicmTJ8fOnTvjzTffjEmTJsWyZeviRz96KRoaGuLpp5+OiIiPf/zjcffddx9US7m/B9ASFTLVR/Rmlitbt3avvTvGjBnD6aefDsDFF1/MqlWrAPjEJz4BwJo1a5gxYwaNjY0MGjSIpqYmnnjiiYO2s2rVKj72sY8xbNgwhg8fzoUXXkhLy08AmDBhAlOmTAFg2rRpbNmypdd1O+jNLFfGju1ee3d0HNpYnB42bFjvN546/PDD333d0NDAvn37er1NB72Z5cr118PQoQe2DR2atPfW1q1befLJJwG45557OOOMMw6YP336dB5//HFee+019u/fz5IlSzjrrLMAGDx4MHv37gXgzDPP5IEHHmDXrl289dZbLF++nELhzN4XWIGD3sxypakJFi6E4oHxuHHJdHoOtFeOP/54br31ViZOnMjrr7/OlVdeecD8kSNHcsMNNzBz5kxOPvlkpk2bxpw5cwBobm7mpJNOoqmpialTp3LppZcyffp0Tj31VK644gomTTql9wVW0OU9Y+vBNx4xs1IvvPACEydO7NY6xR9LrVxZmxq2bNnCeeedx/PPP1+bDXZQzbVuisr9PSStjYhCueV9UTMzy6VaBXweuOvGzKwK48eP77Oj+b7moDezQ0IWu5nroSd/Bwe9mWXekCFD2L59+4AP+0ivRz9kyJBurec+ejPLvNGjR9Pa2kp7e3u9S+kzv/518vz733e+XPEOU93RZdBLWgScB2yLiMll5n8VKA5cGgRMBBojuV/sFuBNYD+wr9IZYTOzzgwePLhbd1Q6FBVHavbFSeRqum7uBGZXmhkRN0XElIiYAnwdeDwOvAH4zHS+Q97MrA66DPqIeALY0dVyqXnAkl5VZGZmNVWzk7GShpIc+d9f0hzAI5LWSmqu1XuZmVn1anky9s+A/9eh2+aMiGiT9O+ARyX9PP2GcJB0R9AMMLYWVx8yMzOgtsMr59Kh2yYi2tLnbcByYHqllSNiYUQUIqLQ2NhYw7LMzAa2mgS9pPcBZwE/KGkbJumI4mtgFnBo/qzMzOwQVs3wyiXADOBoSa3AdcBggIi4LV3sY8AjEfFWyarHAMvT6zUPAu6JiB/WrnQzM6tGl0EfEfOqWOZOkmGYpW2bgZN7WpiZmdWGL4FgZpZzDnozs5xz0JuZ5ZyD3sws5xz0ZmY556A3M8s5B72ZWc456M3Mcs5Bb2aWcw56M7Occ9CbmeWcg97MLOcc9GZmOeegNzPLOQe9mVnOOejNzHLOQW9mlnMOejOznHPQm5nlXJdBL2mRpG2Snq8wf4akNyQ9kz7+qmTebEkbJW2SNL+WhZuZWXWqOaK/E5jdxTI/iYgp6WMBgKQG4FbgXGASME/SpN4Ua2Zm3ddl0EfEE8COHmx7OrApIjZHxB7gXmBOD7ZjZma9UKs++j+W9DNJD0k6MW0bBbxSskxr2laWpGZJLZJa2tvba1SWmZnVIujXAeMi4mTgb4AHerKRiFgYEYWIKDQ2NtagLDMzgxoEfUT8LiJ2pq8fBAZLOhpoA8aULDo6bTMzs37U66CX9AFJSl9PT7e5HVgDHCdpgqTDgLnAit6+n5mZdc+grhaQtASYARwtqRW4DhgMEBG3ARcBV0raB+wG5kZEAPskXQ08DDQAiyJifZ98CjMzq6jLoI+IeV3M/1vgbyvMexB4sGelmZkNDIsXw+rV8M47MH48XH89NDXVbvv+ZayZWR0tXgzNzUnIA7z8cjK9eHHt3sNBb2ZWR9/4BuzadWDbrl1Je6046M3M6mjr1u6194SD3sysjsaO7V57Tzjozczq6PrrYejQA9uGDk3aa8VBb2ZWR01NsHAhHH54Mj1uXDJdy1E3XQ6vNDOzvtXUBH//98nrlStrv30f0ZuZ5ZyD3sws5xz0ZmY556A3M8s5B72ZWc456M2sT8yYkTys/nIV9P6HZWZ2sFwFvZmZHcxBb2aWcw76AcxdXWYDg4O+HzhQzayeugx6SYskbZP0fIX5TZKelfScpH+VdHLJvC1p+zOSWmpZuJmZVaeaI/o7gdmdzH8JOCsiPgR8E1jYYf7MiJgSEYWelVid4j0XH388uediLW/DZWZ2KKvm5uBPSBrfyfx/LZlcDYyuQV3dUumei1DbS32amR2Kat1HfznwUMl0AI9IWiupubMVJTVLapHU0t7e3q037Y97LpqZHapqFvSSZpIE/TUlzWdExFTgXOAqSR+utH5ELIyIQkQUGhsbu/Xe/XHPxZ5yl5JZtgzEwRE1CXpJJwG3A3MiYnuxPSLa0udtwHJgei3er6P+uOdiT1TqUnLYW975ACdbeh30ksYCy4BPRcSLJe3DJB1RfA3MAsqO3Omt/rjnYk+4S8kGIh/gZE81wyuXAE8Cx0tqlXS5pM9K+my6yF8BRwHf7TCM8hhglaSfAU8B/yciftgHn6Ff7rnYE1nuUjLrKz7AyZ5qRt3M62L+FcAVZdo3AycfvEbf6Ot7LvbE2LHJ0Uy59norfrV+553kq/X119d/x2j54AOc7PEvY/tQVruU/NXa+lJWz5kNZA76PpTVLqWsf7UeiKMi8iSrBzgDWZddN9Y7WexSOhS+Wr/9Nrz4YtfLWfb80R/BggVw7bWwZw8ceyx8+ctJexb+m+7enTxnoZZSu3eD1DfbdtAPQFk+d1AUAcOH17sK66l58+D++5PXS5fWt5aOGhqS56z9+2pogP37+2bbueq6WbkyO0fNWeav1mYDS66C3qqT1XMH1jMXXZQ8zCpx100/yOK3jCyeO8i6YphmrSvCrCs+ordMKY7vf+opmD4dli2rd0Vmhz4HvWVGx/H9bW3wta857DuzbBmsXQtPPukdo1XmoLfMKDe+f/duuOGG+tSTdcuWJTvCPXuSae8YrRIHvWVGpXH8v/pV/9ZxqLjhhn8bE17kHWPnBuo3IAe9ZUalcfzHHtu/dZSTxYCotAP0jrG8gfwNyEFvmVFufP973gPz59ennqKsBkSlHWAWdoxZNJC/ATnoLTM6ju8fNQpuvBEuvLC+dWU1IObPT3aEpbKwYyxaujRbQ1EH8jcgj6O3TCmO79+9G5Yvr3c1iawGRHEH+JWvJN82Ro1KQr7eO8asOvbY5NtYufa8c9APYFn+oZQEO3fWu4rEyJHlQ33kyPrXOGsWTJmSvL777uS53jVl1Ze+BH/5l8kF84qGDEnas/A327/fFzWzAWbIEPjgB+tdReLGG5Px/aVDP4cOTdqzUGOx+yYLtWTZV74CH/gAXH558luNceOydcOdjt1wteSgN+tCMQiyGhBZ/maWNQP10h9VnYyVtEjSNkllb+6txHckbZL0rKSpJfMukfSL9HFJrQq3/MriVUibmuC00+Css2DLluyEvFk1qh11cycwu5P55wLHpY9m4O8AJL0fuA44FZgOXCdpRE+LNTOz7qsq6CPiCWBHJ4vMAb4fidXAkZJGAh8BHo2IHRHxOvAone8wzMysxmo1jn4U8ErJdGvaVqn9IJKaJbVIamlvb69RWWZmlpkfTEXEwogoREShsbGx3uWYmeVGrUbdtAFjSqZHp21twIwO7Str9J5m/SprJ4jNqlWrI/oVwKfT0TenAW9ExKvAw8AsSSPSk7Cz0jYzM+snVR3RS1pCcmR+tKRWkpE0gwEi4jbgQeCjwCZgF3BZOm+HpG8Ca9JNLYiIzk7qmplZjVUV9BExr4v5AVxVYd4iYFH3SzMzs1rIzMlYMzPrGw56M7Occ9CbmeWcL2pmZgPKQBwm6yN6M7Oc8xG9mVkG9OU3DR/Rm5nlnIPezCznHPRmZjnnoDczyzkHvZlZzjnozcxyzkFvZpZzDnozs5xz0JuZ5ZyD3sws5xz0ZmY556A3M8u5qoJe0mxJGyVtkjS/zPy/lvRM+nhR0m9L5u0vmbeihrWbmVkVurx6paQG4FbgHKAVWCNpRURsKC4TEV8qWf5zwCklm9gdEVNqVrGZmXVLNUf004FNEbE5IvYA9wJzOll+HrCkFsWZmVnvVRP0o4BXSqZb07aDSBoHTAAeK2keIqlF0mpJF1R6E0nN6XIt7e3tVZRlZmbVqPXJ2LnA0ojYX9I2LiIKwCeBWyT9+3IrRsTCiChERKGxsbHGZZmZDVzVBH0bMKZkenTaVs5cOnTbRERb+rwZWMmB/fdmZtbHqgn6NcBxkiZIOowkzA8aPSPpBGAE8GRJ2whJh6evjwZOBzZ0XNfMzPpOl6NuImKfpKuBh4EGYFFErJe0AGiJiGLozwXujYgoWX0i8D1JvyfZqdxQOlrHzMz6ng7M5WwoFArR0tJS7zLMzA4Zktam50MP4l/GmpnlnIPezCznHPRmZjnnoDczyzkHvZlZzjnozcxyzkFvZpZzDnozs5xz0JuZ5ZyD3sws5xz0ZmY556A3M8s5B72ZWc456M3Mcs5Bb2aWcw56M7Occ9CbmeWcg97MLOeqCnpJsyVtlLRJ0vwy8y+V1C7pmfRxRcm8SyT9In1cUsvizcysa13eHFxSA3ArcA7QCqyRtKLMTb7/KSKu7rDu+4HrgAIQwNp03ddrUr2ZmXWpmiP66cCmiNgcEXuAe4E5VW7/I8CjEbEjDfdHgdk9K9XMzHqimqAfBbxSMt2atnX055KelbRU0phuroukZkktklra29urKMvMzKpRq5Ox/wKMj4iTSI7a7+ruBiJiYUQUIqLQ2NhYo7LMzKyaoG8DxpRMj07b3hUR2yPinXTydmBateuamVnfqibo1wDHSZog6TBgLrCidAFJI0smzwdeSF8/DMySNELSCGBW2mZmZv2ky1E3EbFP0tUkAd0ALIqI9ZIWAC0RsQL4vKTzgX3ADuDSdN0dkr5JsrMAWBARO/rgc5iZWQWKiHrXcJBCoRAtLS31LsPM7JAhaW1EFMrN8y9jzcxyzkFvZpZzDnozs5xz0JuZ5ZyD3sws5xz0ZmY556A3M8s5B72ZWc456M3Mcs5Bb2aWcw56M7Occ9CbmeWcg97MLOcc9GZmOeegNzPLOQe9mVnOOejNzHLOQW9mlnNVBb2k2ZI2StokaX6Z+V+WtEHSs5L+r6RxJfP2S3omfazouK6ZmfWtLm8OLqkBuBU4B2gF1khaEREbShZ7GihExC5JVwI3Ap9I5+2OiCm1LdvMzKpVzRH9dGBTRGyOiD3AvcCc0gUi4scRsSudXA2Mrm2ZZmbWU9UE/SjglZLp1rStksuBh0qmh0hqkbRa0gXdL9HMzHqjy66b7pB0MVAAzippHhcRbZL+EHhM0nMR8csy6zYDzQBjx46tZVlmZgNaNUf0bcCYkunRadsBJP0p8A3g/Ih4p9geEW3p82ZgJXBKuTeJiIURUYiIQmNjY9UfwMzMOldN0K8BjpM0QdJhwFzggNEzkk4BvkcS8ttK2kdIOjx9fTRwOlB6EtfMzPpYl103EbFP0tXAw0ADsCgi1ktaALRExArgJmA4cJ8kgK0RcT4wEfiepN+T7FRu6DBax8zM+pgiot41HKRQKERLS0u9yzAzO2RIWhsRhXLz/MtYM7Occ9CbmeWcg97MLOcc9GZmOeegNzPLOQe9mVnOOejNzHLOQW9mlnMOejOznHPQm5nlnIPezCznHPRmZjnnoDczyzkHvZlZzjnozcxyzkFvZpZzDnozs5xz0JuZ5ZyD3sws56oKekmzJW2UtEnS/DLzD5f0T+n8n0oaXzLv62n7RkkfqWHtZmZWhS6DXlIDcCtwLjAJmCdpUofFLgdej4j/APw18O103UnAXOBEYDbw3XR7ZmbWT6o5op8ObIqIzRGxB7gXmNNhmTnAXenrpcDZkpS23xsR70TES8CmdHtmZtZPBlWxzCjglZLpVuDUSstExD5JbwBHpe2rO6w7qtybSGoGmtPJnZI2VlFbOUcDr/Vw3b7kurrHdXWP6+qePNY1rtKMaoK+X0TEQmBhb7cjqSUiCjUoqaZcV/e4ru5xXd0z0OqqpuumDRhTMj06bSu7jKRBwPuA7VWua2ZmfaiaoF8DHCdpgqTDSE6uruiwzArgkvT1RcBjERFp+9x0VM4E4DjgqdqUbmZm1eiy6ybtc78aeBhoABZFxHpJC4CWiFgB3AHcLWkTsINkZ0C63D8DG4B9wFURsb+PPktRr7t/+ojr6h7X1T2uq3sGVF1KDrzNzCyv/MtYM7Occ9CbmeVcboK+q8s01IukRZK2SXq+3rUUSRoj6ceSNkhaL+kL9a6pSNIQSU9J+lla23+vd01FkhokPS3pf9e7llKStkh6TtIzklrqXU+RpCMlLZX0c0kvSPrjDNR0fPp3Kj5+J+mL9a4LQNKX0n/zz0taImlIzbadhz769LIKLwLnkPwoaw0wLyI21LUwQNKHgZ3A9yNicr3rAZA0EhgZEeskHQGsBS7IyN9LwLCI2ClpMLAK+EJErO5i1T4n6ctAAXhvRJxX73qKJG0BChGRqR8ASboL+ElE3J6O2BsaEb+tc1nvSnOjDTg1Il6ucy2jSP6tT4qI3ekglgcj4s5abD8vR/TVXKahLiLiCZKRSJkREa9GxLr09ZvAC1T4xXJ/i8TOdHJw+qj70Yik0cB/Am6vdy2HAknvAz5MMiKPiNiTpZBPnQ38st4hX2IQ8J70t0hDgV/VasN5Cfpyl2nIRHBlXXql0VOAn9a5lHelXSTPANuARyMiC7XdAnwN+H2d6ygngEckrU0vJZIFE4B24B/S7q7bJQ2rd1EdzAWW1LsIgIhoA/4nsBV4FXgjIh6p1fbzEvTWA5KGA/cDX4yI39W7nqKI2B8RU0h+ST1dUl27vCSdB2yLiLX1rKMTZ0TEVJIrzF6VdhfW2yBgKvB3EXEK8BaQpXNnhwHnA/fVuxYASSNIeiEmAMcCwyRdXKvt5yXofamFbkr7v+8HFkfEsnrXU076Vf/HJJe4rqfTgfPTvvB7gT+R9I/1LenfpEeDRMQ2YDnZuEJsK9Ba8m1sKUnwZ8W5wLqI+E29C0n9KfBSRLRHxF5gGfAfa7XxvAR9NZdpsFR6wvMO4IWIuLne9ZSS1CjpyPT1e0hOsP+8njVFxNcjYnREjCf5t/VYRNTsaKs3JA1LT6iTdo3MAuo+wisifg28Iun4tOlskl/IZ8U8MtJtk9oKnCZpaPr/59kk585qIjNXr+yNSpdpqHNZAEhaAswAjpbUClwXEXfUtypOBz4FPJf2hQNcGxEP1q+kd40E7kpHRPwB8M8RkanhjBlzDLA8yQYGAfdExA/rW9K7PgcsTg++NgOX1bke4N0d4jnAf653LUUR8VNJS4F1JJeLeZoaXg4hF8Mrzcyssrx03ZiZWQUOejOznHPQm5nlnIPezCznHPRmZjnnoDczyzkHvZlZzv1/mKorYkV4Xj4AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "(1.1121024147048486, 0.05917619688972632)\n" ] } ], "source": [ "array = np.array\n", "# only use eval for trusted files!\n", "proton_0p04=eval(open(\"proton_0p04.json\").read())\n", "proton_0p04={x: proton_0p04[x] for x in proton_0p04 if x >= 400} # keep only thermalized\n", "\n", "proton_0p02=eval(open(\"proton_0p02.json\").read())\n", "proton_0p02={x: proton_0p02[x] for x in proton_0p02 if x >= 400} # keep only thermalized\n", "\n", "\n", "def emp_log(Cf):\n", " E = [np.log(Cf[t].real / Cf[t+1].real) for t in range(2,11)]\n", " return np.array(E)\n", "\n", "def proton(data, nbin, tag, t0):\n", " corrs = bin_data(np.array([data[x][\"proton\"] for x in data]),nbin)\n", " print(len(corrs))\n", " E_p = emp_log(np.mean(corrs,axis=0))\n", " E_p_err = jackknife_variance_estimator(emp_log, corrs)**0.5\n", " \n", " #print(np.mean(corrs,axis=0))\n", " #img = np.mean(corrs,axis=0).imag\n", " #img_err = jackknife_variance_estimator(lambda x: x.imag, corrs)**0.5\n", " #for i in range(len(img)):\n", " # print(i,img[i]/img_err[i])\n", " \n", " t1 = 9\n", " E_p_plateau = plateau(np.mean(corrs,axis=0), t0, t1, E_p_err, emp_log)\n", " E_p_plateau_err = jackknife_variance_estimator(lambda c: plateau(c, t0, t1, E_p_err, emp_log), corrs)**0.5\n", "\n", " fig, ax = plt.subplots()\n", "\n", " plt.ylim(0.0,2.0)\n", " plt.title(tag)\n", "\n", " X=len(E_p)\n", " ax.errorbar(np.arange(X),E_p,E_p_err,marker=\"o\",ls=\"\",color=\"blue\",label=\"proton\")\n", " ax.fill_between(np.arange(t0,t1),[E_p_plateau[0]-E_p_plateau_err[0]]*(t1-t0),\n", " [E_p_plateau[0]+E_p_plateau_err[0]]*(t1-t0),alpha=0.1,color=\"blue\")\n", " \n", " plt.legend()\n", "\n", " plt.show()\n", " \n", " return (E_p_plateau[0], E_p_plateau_err[0])\n", "\n", "\n", "\n", "mproton_0p04 = proton(proton_0p04,1,\"m = 0.04\", 3)\n", "\n", "mproton_0p02 = proton(proton_0p02,1,\"m = 0.02\", 3)\n", "print(mproton_0p02)" ] }, { "cell_type": "code", "execution_count": null, "id": "46792b04", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "id": "6cfbc2a1", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "76c08dd5", "metadata": {}, "source": [ "Finally, let us consider the pion, rho, and proton masses for the ensembles generated so far." ] }, { "cell_type": "code", "execution_count": 183, "id": "4c066df5", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEICAYAAAC3Y/QeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAaMklEQVR4nO3dfZRU9Z3n8fdHHoMoOtAzx8hD41nFxgYUepEdn9o1m1GXyJjxTPQUo2bU3ohmJjOZGcwwJ5k4i0c3exx3dnAyrRNNTIsPccPixhw3owHjRlwaRRQZIiENNGYVkZDRRgT87h9VDUXbTVV33eqq6vt5nVOn6/7urXu/9aP49O17b/2uIgIzM0uX4ypdgJmZDT6Hv5lZCjn8zcxSyOFvZpZCDn8zsxRy+JuZpZDD36xIkuolhaThA3z9X0q6P+m6zAbC4W8VIalD0j5J7+U9/r7Aa5oldQ5WjaXordaIuCMibqxUTWb5BrQHY5aQz0TEPye5QknDI+Jgkus0G4q8529VRdI/SHoib/ouSc9IOh74IfDJvL8UPinpryV9T9J3Jf0auF7SXEkvSPqVpF9K+ntJI/PWGZL+SNJWSe9I+oak43LzjpP0V5K2SXpb0nckjeuj1s9L2iTpX3Pr+k+59mPV+t28118haWOuzlWSGvLmdUj6M0kbJO2V9Kik0Ql3t6WYw9+qzZeBGZKul3QBcANwXUS8D1wGvBkRY3OPN3OvWQB8DzgJaAMOAX8CTAD+HXAJsKjHdq4EmoDZudf/Ya79+tzjYuA0YCzQ1+Got4H5wInA54G/lTS7QK0ASDoDWA58CagDngKezP8lBfw+cCkwFZiZq8ssEQ5/q6QVub3e7sdNEdEF/AFwN/Bd4IsRUeg4/wsRsSIiPoqIfRGxLiLWRMTBiOgA/hG4qMdr7oqIdyNiO3APcE2uPQPcHRFbI+I94CvA1b2d5I2IH0TEzyNrNfC/gQuKfO+fA34QET+KiAPAfwU+Afx23jJ/FxFvRsS7wJPA2UWu26wgH/O3Svrd3o75R8SLkrYCvwk8VsR6duRP5Paq7ya7Zz+G7Od83TFesw34ZO75J3PT+fOGA7/Vc6OSLgO+BpxBdkdqDPBqEfV+bDsR8ZGkHcCpecv8v7znXXk1mpXMe/5WdSTdAowC3gT+Im9WX0PQ9mz/B+BfgNMj4kTgLwH1WGZS3vPJuW2R+zmlx7yDwFs9ahwFPEF2j/23IuIksoduurdTaLjco7YjSbmadhZ4nVkiHP5WVXJ77f8ZWEj28M9fSDo7N/stYHxfJ2DznAD8GnhP0pnAzb0s8+eSTpY0Cfhj4NFc+3LgTyRNlTQWuAN4tJcriEaS/QW1CziY+yvg03nzC9X6GPAfJV0iaQTZcx37gZ8WeG9mifBhH6ukJyUdypv+EdnDHndFxCuQ/WIU8JCkpoj4F0nLga2ShgHT+1jvnwGtZP9qeJlssP/7Hsv8T7KHgsYBDwL/lGv/FtnDK88Bo4GngS/23EBE/KukPyIb4qPIHpNfmTf/mLVGxGZJC4H/nnvP68le+vphH+/JLFHyzVwsbSQF2UNCWypdi1ml+LCPmVkKOfzNzFLIh33MzFLIe/5mZilUsat9JkyYEPX19ZXavJlZTVq3bt07EVFX6noqFv719fW0t7dXavNmZjVJ0rbCSxXmwz5mZink8DczSyGHv5lZClXV8A4HDhygs7OTDz74oNKlVIXRo0czceJERowYUelSzGyIqarw7+zs5IQTTqC+vp7sIIfpFRHs3r2bzs5Opk6dWulyzGyIqarDPh988AHjx49PffADSGL8+PH+K8jMyqKqwh9w8OdxX5hZuVRd+PdXc3P2YWZmxav58B8sN954I6+//nqlyzAzS0RNh39bG6xZA6tXQ319drpc7r//fqZP7+veIWZmtaVmw7+tDVpaYP/+7PS2bdnpUn8BdHR0cOaZZ5LJZGhoaOCqq66iq6uL5ubmw8NRLF++nBkzZtDY2MjixYsPv3bs2LEsWbKEWbNmMW/ePN56662+NmNmVlE1G/5LlkBX19FtXV3Z9lJt3ryZRYsWsWnTJk488UTuvffew/PefPNNFi9ezLPPPsv69etZu3YtK1asAOD9999n3rx5vPLKK1x44YXcd999pRdjZlYGBcNf0rckvS3ptT7mZyRtkPSqpJ9KmpV8mR+3fXv/2vtj0qRJnHfeeQAsXLiQ559//vC8tWvX0tzcTF1dHcOHDyeTyfDcc88BMHLkSObPnw/AnDlz6OjoKL0YM7MyKGbP/0Hg0mPM/wVwUUTMAP6G7I2zy27y5P6190fPSyyLveRyxIgRh5cdNmwYBw8eLL0YM7MyKBj+EfEc8O4x5v80IvbkJtcAExOq7ZiWLoUxY45uGzMm216q7du388ILLwDw8MMPc/755x+eN3fuXFavXs0777zDoUOHWL58ORdddFHpGzUzG0RJH/O/AfhhXzMltUhql9S+a9eukjaUyUBrK4walZ2eMiU7ncmUtFoApk2bxrJly2hoaGDPnj3cfPPNh+edcsop3HnnnVx88cXMmjWLOXPmsGDBgtI3amY2iIq6h6+keuB/RUTjMZa5GLgXOD8idhdaZ1NTU/S8mcumTZtoaGgoWE++7i94rVrVr5f1qaOjg/nz5/Paa72e4hh0A+kTMxu6JK2LiKZS15PIwG6SZgL3A5cVE/xJSir0zczSpOTDPpImA/8D+IOI+FnpJVVWfX191ez1mxk0P9hM84PNlS5jyCm45y9pOdAMTJDUCXwNGAEQEd8EvgqMB+7NXelyMIk/SczMrHwKhn9EXFNg/o3AjYlVZDZEde+9rrp+VUXrqCVtr7axpnMN+w/tp/6eepZespTMjASu6rDa/YavmQ1tba+20fJkC/sPZcdw2bZ3Gy1PttD2ahkH8UqR2g9/j+lsNiQteWYJXQeOHsOl60AXS55JYAwXGwLhX0YdHR00NvZ5dauZldH2vb2P1dJXu/VPbYd/mcd0jgg++uijRNdpZsWZPK73sVr6arf+qd3wL9OYzh0dHUybNo1rr72WxsZG9u3bx0033cRZZ53Fpz/9afbt2wfA+vXrmTdvHjNnzuTKK69kz549BdZsZv2x9JKljBlx9BguY0aMYeklCYzhYjUc/mUc0/mNN95g0aJFbNy4kR07dnDLLbewceNGTjrpJJ544gkArr32Wu666y42bNjAjBkz+PrXv17yds3siMyMDK2faWXUsOwYLlPGTaH1M62+2ichtRv+ZRzTecqUKcybNw+AqVOncvbZZwNHhmneu3cvv/rVrw4P6HbdddcdHtbZrDfdlyyu3raa+nvqfcVKkTIzMsybOI+LplxEx5c6HPwJSmR4h4qYPDl7qKe39hIdf/zxh5+P6h45juwwzd2HfcyK1dcli4DDrAj+XkR51O6efznHdC5g3LhxnHzyyfzkJz8B4KGHHvKwztYnX7Jo1ah29/y7x26+4YbsSd8pU7LBn8SYzkX49re/zRe+8AW6uro47bTTeOCBBwZlu1Z7fMmiVaOihnQuh6SGdE58TOcq4yGda1/9PfVs2/vxQ5RTxk2h40sdg1+Q1bSkhnSu3cM+3VatGrLBb0ODL1m0alT74W9W5XzJolWj2j3mb1ZDMjMy3LfuPsBXr1h18J6/mVkKOfzNzFKo5sPft3gzM+u/mg//anPHHXdUugQzs4JqOvwrNV7KoUOH+pzn8DezWlCz4V+uW7x1dHRw5plnkslkaGho4KqrrqKrq4v6+noWL17M7Nmzefzxx1m+fDkzZsygsbGRxYsXA3Dbbbexb98+zj77bDK5bxrffffdNDY20tjYyD333HN4Gw0NDb0OFW1D16rrV/lKH6saNRv+5RwvZfPmzSxatIhNmzZx4okncu+99wIwfvx4XnrpJS688EIWL17Ms88+y/r161m7di0rVqzgzjvv5BOf+ATr16+nra2NdevW8cADD/Diiy+yZs0a7rvvPl5++WUgO2x0b0NFm5kNhpoN/3KOlzJp0iTOO+88ABYuXMjzzz8PwOc+9zkA1q5dS3NzM3V1dQwfPpxMJtPrkM7PP/88V155Jccffzxjx47ls5/97OHB4HobKtrMbLDUbPiX8xZvknqdzh/quVQ9h4o+ePBgYus2MyukZsO/nOOlbN++nRdeeAGAhx9+mPPPP/+o+XPnzmX16tW88847HDp0iOXLlx8e0nnEiBEcOHAAgAsuuIAVK1bQ1dXF+++/z/e//30uuOCCkuszMytVwfCX9C1Jb0t6rY/5kvR3krZI2iBpdvJlflw5x0uZNm0ay5Yto6GhgT179nDzzTcfNf+UU07hzjvv5OKLL2bWrFnMmTOHBQsWANDS0sLMmTPJZDLMnj2b66+/nrlz53Luuedy4403cs4555Rcn5lZqQoO6SzpQuA94DsR0djL/MuBLwKXA+cC/y0izi204aSGdO7+gldSV1F0dHQwf/58Xnut1991g85DOptZvqSGdC44sFtEPCep/hiLLCD7iyGANZJOknRKRPyy1OKK4UvnzMz6L4lj/qcCO/KmO3NtNam+vr5q9vrNzMplUE/4SmqR1C6pfdeuXb0uU6k7i1Uj94WZlUsS4b8TmJQ3PTHX9jER0RoRTRHRVFdX97H5o0ePZvfu3Q49ssG/e/duRo8eXelSzGwISuJmLiuBWyU9QvaE796BHu+fOHEinZ2d9PVXQdqMHj2aiRMnVroMMxuCCoa/pOVAMzBBUifwNWAEQER8E3iK7JU+W4Au4PMDLWbEiBFMnTp1oC83M7MiFXO1zzUF5gdwS2IVmZlZ2dXsN3zNzGzgHP5mZink8DczSyGHv5lZCjn8zcxSyOFvZpZCDn8zsxRy+JuZpZDD38wshRz+ZmYp5PA3M0shh7+ZWQo5/M3MUsjhb2aWQg5/M7MUcvibmaWQw9/MLIUc/mZmKeTwNzNLIYe/mVkKOfzNzFLI4W9mlkIOfzOrXm1tUF8Pxx2X/dnWVumKhozhlS7AzKxXbW3Q0gJdXdnpbduy0wCZTOXqGiK8529m1WnJkiPB362rK9tuJSsq/CVdKmmzpC2Sbutl/mRJP5b0sqQNki5PvlQzS5Xt2/vXbv1SMPwlDQOWAZcB04FrJE3vsdhfAY9FxDnA1cC9SRdqZikzeXL/2q1fitnznwtsiYitEfEh8AiwoMcyAZyYez4OeDO5Es0slZYuhTFjjm4bMybbbiUrJvxPBXbkTXfm2vL9NbBQUifwFPDF3lYkqUVSu6T2Xbt2DaBcM0uNTAZaW2HKFJCyP1tbfbI3IUmd8L0GeDAiJgKXAw9J+ti6I6I1Ipoioqmuri6hTZvZkJXJQEcHfPRR9qeDPzHFhP9OYFLe9MRcW74bgMcAIuIFYDQwIYkCzcwsecWE/1rgdElTJY0ke0J3ZY9ltgOXAEhqIBv+Pq5jZlalCoZ/RBwEbgWeBjaRvapno6TbJV2RW+zLwE2SXgGWA9dHRJSraLOa1NycfZhVgaK+4RsRT5E9kZvf9tW8568D5yVbmpmZlYu/4Ws2GNraYM0aWL3aY9RYVXD4m5Vb9xg1+/dnp7vHqPEvAKsgh79ZuXmMGqtCDn+zcvMYNVaFHP5m5eYxaqwKOfzNys1j1FgVcviblVv3GDWjRmWnPUaNVQHfyctsMGQyDnurKt7zNzNLIYe/mVkKOfzNzFLI4W9mlkIOfzOzFHL4m5mlkMPfzCyFHP5mZink8DczSyGHv5lZCjn8zcxSyOFvZpZCDn8zsxRy+JuZpZDD38wshRz+ZmYpVFT4S7pU0mZJWyTd1scyvy/pdUkbJT2cbJlmZpakgnfykjQMWAb8B6ATWCtpZUS8nrfM6cBXgPMiYo+k3yxXwWZmVrpi9vznAlsiYmtEfAg8AizoscxNwLKI2AMQEW8nW6aZmSWpmPA/FdiRN92Za8t3BnCGpP8jaY2kS3tbkaQWSe2S2nft2jWwis3MrGRJnfAdDpwONAPXAPdJOqnnQhHRGhFNEdFUV1eX0KbNzKy/ign/ncCkvOmJubZ8ncDKiDgQEb8Afkb2l4GZmVWhYsJ/LXC6pKmSRgJXAyt7LLOC7F4/kiaQPQy0NbkyzcwsSQXDPyIOArcCTwObgMciYqOk2yVdkVvsaWC3pNeBHwN/HhG7y1W0mZmVRhFRkQ03NTVFe3t7RbZtZlarJK2LiKZS1+Nv+JqZpZDD38wshRz+ZmYp5PA3M0shh7+ZWQo5/M3MUsjhb2aWQg5/M7MUcvibWdVqa4P6ejjuuOzPtrZKVzR0FLyZi5lZJbS1QUsLdHVlp7dty04DZDKVq2uo8J6/mVWlJUuOBH+3rq5su5XO4W9mVWn79v61W/84/M2sKk2e3L926x+Hv5lVpaVLYcyYo9vGjMm2W+kc/mZWlTIZaG2FKVNAyv5sbfXJ3qT4ah8zq1qZjMO+XLznb2aWQg5/M7MUcvibmaWQw9/MLIUc/mZmKeTwNzNLIYe/mVkKOfzNzFLI4W9mlkJFhb+kSyVtlrRF0m3HWO73JIWkpuRKNDOzpBUMf0nDgGXAZcB04BpJ03tZ7gTgj4EXky7SzMySVcye/1xgS0RsjYgPgUeABb0s9zfAXcAHCdZnZmZlUEz4nwrsyJvuzLUdJmk2MCkifnCsFUlqkdQuqX3Xrl39LtasVvletFZtSj7hK+k44G7gy4WWjYjWiGiKiKa6urpSN21WE7rvRbttG0QcuRetfwFYJRUT/juBSXnTE3Nt3U4AGoFVkjqAecBKn/Q1y/K9aK0aFRP+a4HTJU2VNBK4GljZPTMi9kbEhIioj4h6YA1wRUS0l6Visxrje9FaNSoY/hFxELgVeBrYBDwWERsl3S7pinIXaFbrfC9aq0ZF3ckrIp4CnurR9tU+lm0uvSyzoWPp0uwx/vxDP74XrVWav+FrVma+F61VI9/D12wQ+F60Vm28529mlkIOfzOzFHL4m5mlkMPfzCyFHP5mZink8DczSyGHv5lZCjn8zcxSyOFvZpZCDn8zsxRy+JuZpZDD38wshRz+ZmYp5PA3M0shh7+ZWQo5/M3MUsjhb2aWQg5/M7MUcvibmaWQw9/MLIUc/mZmKeTwNzNLoaLCX9KlkjZL2iLptl7m/6mk1yVtkPSMpCnJl2pmZkkpGP6ShgHLgMuA6cA1kqb3WOxloCkiZgLfA/5L0oWamVlyitnznwtsiYitEfEh8AiwIH+BiPhxRHTlJtcAE5Mt08zMklRM+J8K7Mib7sy19eUG4Ie9zZDUIqldUvuuXbuKr9LMzBKV6AlfSQuBJuAbvc2PiNaIaIqIprq6uiQ3bWZm/TC8iGV2ApPypifm2o4i6VPAEuCiiNifTHlmZlYOxez5rwVOlzRV0kjgamBl/gKSzgH+EbgiIt5OvkwzM0tSwfCPiIPArcDTwCbgsYjYKOl2SVfkFvsGMBZ4XNJ6SSv7WJ2ZmVWBYg77EBFPAU/1aPtq3vNPJVyXmZmVkb/ha2aWQg5/M7MUcvibmaWQw9/MLIUc/mZmKeTwNzNLIYe/mVkKOfzNzFLI4W9mlkIOfzOzFHL4m5mlkMPfzCyFHP5mZink8DczSyGHv5lZCjn8zcxSyOFvZpZCDn8zsxRy+JuZpZDD38wshRz+ZmYp5PA3M0shh7+ZWQo5/M3MUsjhb2aWQkWFv6RLJW2WtEXSbb3MHyXp0dz8FyXVJ16pmZklpmD4SxoGLAMuA6YD10ia3mOxG4A9EfFvgL8F7kq6UDMzS04xe/5zgS0RsTUiPgQeARb0WGYB8O3c8+8Bl0hScmWamVmShhexzKnAjrzpTuDcvpaJiIOS9gLjgXfyF5LUArTkJvdLem0gRQ9BE+jRVynmvjjCfXGE++KIaUmspJjwT0xEtAKtAJLaI6JpMLdfrdwXR7gvjnBfHOG+OEJSexLrKeawz05gUt70xFxbr8tIGg6MA3YnUaCZmSWvmPBfC5wuaaqkkcDVwMoey6wErss9vwp4NiIiuTLNzCxJBQ/75I7h3wo8DQwDvhURGyXdDrRHxErgn4CHJG0B3iX7C6KQ1hLqHmrcF0e4L45wXxzhvjgikb6Qd9DNzNLH3/A1M0shh7+ZWQqVJfxLGQ5C0ldy7Zsl/U456htMA+0LSfWS9klan3t8c9CLT1gRfXGhpJckHZR0VY9510l6I/e4rudra02JfXEo73PR8+KLmlNEX/yppNclbZD0jKQpefOGzOeixH7o/2ciIhJ9kD0p/HPgNGAk8Aowvccyi4Bv5p5fDTyaez49t/woYGpuPcOSrnGwHiX2RT3wWqXfwyD3RT0wE/gOcFVe+28AW3M/T849P7nS76kSfZGb916l38Mg98XFwJjc85vz/o8Mmc9FKf0w0M9EOfb8SxkOYgHwSETsj4hfAFty66tVHhrjiIJ9EREdEbEB+KjHa38H+FFEvBsRe4AfAZcORtFlUkpfDDXF9MWPI6IrN7mG7HeNYGh9LkrphwEpR/j3NhzEqX0tExEHge7hIIp5bS0ppS8Apkp6WdJqSReUu9gyK+XfNo2fi2MZLald0hpJv5toZYOvv31xA/DDAb62mpXSDzCAz8SgDu9g/fJLYHJE7JY0B1gh6ayI+HWlC7OKmxIROyWdBjwr6dWI+Hmliyo3SQuBJuCiStdSSX30Q78/E+XY8y9lOIhiXltLBtwXuUNfuwEiYh3Z44FnlL3i8inl3zaNn4s+RcTO3M+twCrgnCSLG2RF9YWkTwFLgCsiYn9/XlsjSumHgX0mynDiYjjZEy9TOXLi4qwey9zC0Sc5H8s9P4ujT/hupbZP+JbSF3Xd753sSaCdwG9U+j2Vsy/yln2Qj5/w/QXZk3on556ntS9OBkblnk8A3qDHicFaehT5f+Qcsjs/p/doHzKfixL7YUCfiXK9kcuBn+UKXZJru53sbyuA0cDjZE/o/l/gtLzXLsm9bjNwWaX/USrVF8DvARuB9cBLwGcq/V4GoS/+Ldljne+T/UtwY95r/zDXR1uAz1f6vVSqL4DfBl7NhcOrwA2Vfi+D0Bf/DLyV+7+wHlg5FD8XA+2HgX4mPLyDmVkK+Ru+ZmYp5PA3M0shh7+ZWQo5/M3MUsjhb2aWQg5/M7MUcvibmaXQ/wcfCRKt5A5a6AAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "sqrt_quark_masses = [0.02**0.5, 0.04**0.5]\n", "rho_masses = [mrho_0p02[0], mrho_0p04[0]]\n", "rho_masses_err = [mrho_0p02[1], mrho_0p04[1]]\n", "pi_masses = [mpi_0p02[0], mpi_0p04[0]]\n", "pi_masses_err = [mpi_0p02[1], mpi_0p04[1]]\n", "proton_masses = [mproton_0p02[0], mproton_0p04[0]]\n", "proton_masses_err = [mproton_0p02[1], mproton_0p04[1]]\n", "\n", "def extrapolation():\n", " \n", " fig, ax = plt.subplots()\n", "\n", " plt.ylim(0.0,1.3)\n", " plt.xlim(0,0.25)\n", " \n", " plt.title(\"Extrapolation\")\n", "\n", " ax.errorbar(sqrt_quark_masses,pi_masses,pi_masses_err,marker=\"o\",ls=\"\",color=\"blue\",label=\"pion\")\n", " ax.errorbar(sqrt_quark_masses,rho_masses,rho_masses_err,marker=\"o\",ls=\"\",color=\"red\",label=\"rho\")\n", " ax.errorbar(sqrt_quark_masses,proton_masses,proton_masses_err,marker=\"o\",ls=\"\",color=\"green\",label=\"proton\")\n", " plt.legend()\n", "\n", " plt.show()\n", "\n", "extrapolation()" ] }, { "cell_type": "markdown", "id": "d8a41559", "metadata": {}, "source": [ "To the degree that the $\\rho$ is a stable particle for the given mass configuration, one could then use, e.g., the ratio of $m_\\rho/m_\\pi$ to tune the physical light quark mass $m_l$ and use $m_\\pi / a = m_\\pi^{\\rm experiment}$ to define the lattice spacing. The proton mass at the corresponding mass point is then a prediction of the theory.\n", "\n", "\n", "## $N_f=2+1$\n", "\n", "By introducing a strange quark, new mesonic and baryonic states appear. Particularly important are the\n", "pseudoscalar mesons combining up/down with a strange quark (**kaons**) with a mass of approximately 500 MeV and the vector meson made of a strange quark-anti-quark pair (${\\bf\\phi}$). The latter is a resonance albeit with a narrow width.\n", "\n", "For the baryons a particularly important state is the $J=3/2$, $P=1$ state ${\\bf\\Omega^-}$ made of three strange quarks.\n", "Its mass is known experimentally to very high precision ($1672.45(29)$ MeV), is stable and receives small isospin-breaking and QED corrections. In calculations where a strange quark has been introduced, it is often used to set the lattice spacing.\n", "\n", "## $N_f=2+1+1+1$\n", "\n", "The physics of charm and bottom quarks will be discussed briefly in the next chapter." ] }, { "cell_type": "code", "execution_count": null, "id": "c383b35f", "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 }