{ "cells": [ { "cell_type": "markdown", "id": "ff8465b6", "metadata": {}, "source": [ "| [**Lattice QCD I**](https://homepages.uni-regensburg.de/~lec17310/teaching/wise2122/lqft.html) | - | [**Christoph Lehner**](https://go.ur.de/lehner) |\n", "|:-|--|-:|\n", "\n", "# Chapter 6: gauge theories on a lattice\n", "\n", "In this chapter we explore how we can construct a field theory on a lattice that keeps the gauge symmetry intact." ] }, { "cell_type": "markdown", "id": "ae53a03f", "metadata": {}, "source": [ "## Covariant derivatives in the continuum\n", "\n", "In the last chapter we considered generators $G(x)$ of a gauge symmetry group, which in the case of SU(3) are traceless Hermitian matrices ($G^\\dagger=G$, $\\textrm{Tr} G=0$). We found that we can define covariant derivatives $D_\\mu$ that transformed as\n", "$$\n", "\\begin{align}\n", "D_\\mu \\to [\\mathbb{1}+i G(x)] D_\\mu [\\mathbb{1}-i G(x)]\n", "\\end{align}\n", "$$\n", "under infinitesimal gauge transformations. Finite gauge transformations can then be written as\n", "$$\n", "\\begin{align}\n", "D_\\mu \\to V(x) D_\\mu V(x)^\\dagger\n", "\\end{align}\n", "$$\n", "with\n", "$$\n", "\\begin{align}\n", "V(x) = \\exp\\left( i G(x) \\right) \\,.\n", "\\end{align}\n", "$$\n", "These then allowed for the definition of a gauge invariant kinetic term of a fermionic field $\\psi(x)$\n", "that transforms as\n", "$$\n", "\\begin{align}\n", "\\psi(x) & \\to V(x) \\psi(x) \\,, \\\\\n", "\\bar\\psi(x) & \\to \\bar\\psi(x) V(x)^\\dagger \n", "\\end{align}\n", "$$\n", "using\n", "$$\n", "\\begin{align}\n", "{\\cal L} &= \\bar\\psi i \\gamma^\\mu D_\\mu \\psi \\,.\n", "\\end{align}\n", "$$" ] }, { "cell_type": "markdown", "id": "93df66bd", "metadata": {}, "source": [ "## Discrete derivatives on a lattice\n", "\n", "We now consider fields $\\psi(x)$ that live on a discrete lattice with $x\\in\\Lambda$ and $\\Lambda = \\{ a n \\vert n \\in \\mathbb{Z}^d \\}$ with lattice dimension $d$ and lattice spacing $a$. Then an approximations of the partial derivatives in direction $\\mu$ can be writte, e.g., as\n", "$$\n", "\\begin{align}\n", "a\\overrightarrow{\\partial_\\mu} \\psi(x) &= \\psi(x+a\\hat{\\mu}) - \\psi(x) \\,, \\\\\n", "a\\overleftarrow{\\partial_\\mu} \\psi(x) &= \\psi(x) - \\psi(x-a\\hat{\\mu})\n", "\\end{align}\n", "$$\n", "with unit vectors $\\hat{\\mu}$ in direction $\\mu$. If we imagine $\\psi$ being a sufficiently differentiable field in $\\mathbb{R}^d$\n", "instead, we find\n", "$$\n", "\\begin{align}\n", "\\overrightarrow{\\partial_\\mu} \\psi(x) &= \\frac1a(\\psi(x+a\\hat{\\mu}) - \\psi(x)) = \\partial_\\mu \\psi(x) + O(a) \\,, \\\\\n", "\\overleftarrow{\\partial_\\mu} \\psi(x) &= \\frac1a(\\psi(x) - \\psi(x-a\\hat{\\mu}) = \\partial_\\mu \\psi(x) + O(a)\n", "\\end{align}\n", "$$\n", "and if we take a symmetric combination\n", "$$\n", "\\begin{align}\n", "\\overleftrightarrow{\\partial_\\mu} = \\frac12 (\\overleftarrow{\\partial_\\mu} + \\overrightarrow{\\partial_\\mu})\n", "\\end{align}\n", "$$\n", "we find\n", "$$\n", "\\begin{align}\n", "\\overleftrightarrow{\\partial_\\mu} \\psi(x) &= \\frac1{2a}(\\psi(x+a\\hat{\\mu}) - \\psi(x-a\\hat{\\mu})) = \\partial_\\mu \\psi(x) + O(a^2) \\,.\n", "\\end{align}\n", "$$\n", "\n", "***Homework***: Show this for the symmetric derivative by performing a corresponding Taylor expansion." ] }, { "cell_type": "markdown", "id": "9d6d9e95", "metadata": {}, "source": [ "## Covariant derivatives on a lattice\n", "\n", "We have just seen that on a lattice we need to approximate derivatives with finite differences of fields shifted in discrete steps on the lattice. We can make this more explicit by defining shift operators\n", "$$\n", "\\begin{align}\n", "S^+_\\mu \\psi(x) &= \\psi(x + a\\hat{\\mu}) \\,, \\\\\n", "S^-_\\mu \\psi(x) &= \\psi(x - a\\hat{\\mu})\n", "\\end{align}\n", "$$\n", "such that\n", "$$\n", "\\begin{align}\n", "\\overleftrightarrow{\\partial_\\mu} \\psi(x) = \\frac1{2a}(S^+_\\mu - S^-_\\mu)\\psi(x)\n", "\\end{align}\n", "$$\n", "\n", "If we assume again that under a gauge transformation each field transforms as\n", "$$\n", "\\begin{align}\n", "\\psi(x) \\to V(x) \\psi(x)\n", "\\end{align}\n", "$$\n", "the partial derivatives $\\overleftrightarrow{\\partial_\\mu}\\psi(x)$ again would not have a well-defined behavior under\n", "gauge transformations. What we would need is a covariant shift operation\n", "$$\n", "\\begin{align}\n", "C^+_\\mu \\psi(x) = U_\\mu(x) \\psi(x + a \\hat{\\mu})\n", "\\end{align}\n", "$$\n", "which should obey\n", "$$\n", "\\begin{align}\n", "C^+_\\mu \\psi(x) \\to V(x) C^+_\\mu \\psi(x)\n", "\\end{align}\n", "$$\n", "under a gauge transformation. This can be achieved if the $U_\\mu(x)$ above were to\n", "transform as\n", "$$\n", "\\begin{align}\n", "U_\\mu(x) \\to V(x) U_\\mu(x) V(x + a\\hat{\\mu})^\\dagger\n", "\\end{align}\n", "$$\n", "under gauge transformations. Then we could also define a covariant backwards shift operation\n", "$$\n", "\\begin{align}\n", "C^-_\\mu \\psi(x) = U^\\dagger_\\mu(x - a\\hat{\\mu}) \\psi(x - a \\hat{\\mu})\n", "\\end{align}\n", "$$\n", "that would transform as\n", "$$\n", "\\begin{align}\n", "C^-_\\mu \\psi(x) \\to V(x) C^-_\\mu \\psi(x)\n", "\\end{align}\n", "$$\n", "under gauge transformations.\n", "\n", "Then by replacing all shifts by covariant shifts, we can construct a **covariant derivative**\n", "$$\n", "\\begin{align}\n", "D_\\mu \\psi(x) = \\frac1{2a}(C^+_\\mu - C^-_\\mu)\\psi(x) = \\frac1{2a}\\left(U_\\mu(x) \\psi(x + a \\hat{\\mu}) - U^\\dagger_\\mu(x - a\\hat{\\mu}) \\psi(x - a \\hat{\\mu})\\right)\\,.\n", "\\end{align}\n", "$$\n", "\n", "As operators, we can write\n", "$$\n", "\\begin{align}\n", "C^+_\\mu &= U_\\mu(x) S^+_\\mu \\,, \\\\\n", "C^-_\\mu &= S^-_\\mu U_\\mu^\\dagger(x) \\,.\n", "\\end{align}\n", "$$\n", "We can also consider the $C^\\pm_\\mu$ as **parallel transporters** of the gauge group." ] }, { "cell_type": "markdown", "id": "98fda4cc", "metadata": {}, "source": [ "The parallel transport matrices $U_\\mu(x)$ therefore seem to play the role of transporting a field $\\psi$ between $x$ and $x+a\\hat{\\mu}$. It is therefore only natural to assume geometrically that $U_\\mu(x)$ corresponds to a **link** on the lattice between those two sites.\n", "\n", "If we furthermore make the ansatz\n", "$$\n", "\\begin{align}\n", "U_\\mu(x) &= \n", "\\left(1 + i g \\frac{a}{n} A_\\mu(x)\\right) \\left( 1 + i g \\frac{a}{n} A_\\mu\\left(x+\\frac{a}{n}\\hat{\\mu}\\right)\\right) \\cdots\n", "\\left( 1 + i g \\frac{a}{n} A_\\mu\\left(x+\\frac{a(n-1)}{n}\\hat{\\mu}\\right)\\right) \\\\\n", "&\\stackrel{n\\to\\infty}{=}\n", "{\\cal P} \\exp\\left(i g \\int_{0}^{a} ds A_\\mu(x + s\\hat{\\mu})\\right) = \\exp\\left(i g a A_\\mu\\left(x + \\frac12a\\hat{\\mu}\\right) + O(a^3)\\right)\n", "\\end{align}\n", "$$\n", "with gauge fields $A_\\mu$ living on the links, we find\n", "$$\n", "\\begin{align}\n", "D_\\mu \\psi(x) = (\\partial_\\mu + i g A_\\mu(x)) \\psi(x) + O(a^2) \\,.\n", "\\end{align}\n", "$$\n", "The operator ${\\cal P}$ is the **path-ordering** operator that arises from the limit of multiplying\n", "the individual factors.\n", "\n", "Note that $g$ is the coupling constant that we introduced at the end of the previous chapter.\n", "\n", "***Homework***: Show this equation by explicit Taylor expansion in $a$." ] }, { "cell_type": "markdown", "id": "1b0fd59e", "metadata": {}, "source": [ "## Pure gauge theory\n", "\n", "Consider the smallest possible parallel transport around a $1x1$ square in the $\\mu$-$\\nu$ plane,\n", "the **plaquette**,\n", "$$\n", "\\begin{align}\n", "U_{\\mu\\nu}(x) &= U_\\mu(x) U_\\nu(x+a \\hat{\\mu}) U_\\mu^\\dagger(x + a \\hat{\\nu}) U_\\nu^\\dagger(x)\n", "\\end{align}\n", "$$\n", "transforming to\n", "$$\n", "\\begin{align}\n", "U_{\\mu\\nu}(x) &\\to V(x) U_{\\mu\\nu}(x) V(x)^\\dagger \\,.\n", "\\end{align}\n", "$$\n", "We can then create a gauge-invariant contribution to an action by taking the trace of $U_{\\mu\\nu}$.\n", "Before doing this, however, let us consider what this corresponds to in terms of gauge fields $A_\\mu$.\n", "The plaquette is also the simplest of the **Wilson loops**, i.e., parallel transports around a closed loop.\n", "\n", "If we use the Baker-Campbell-Hausdorff formula\n", "$$\n", "\\begin{align}\n", "\\exp(a A)\\exp(a B) = \\exp\\left(a A + a B + \\frac12 a^2 [A,B] + O(a^3) \\right)\n", "\\end{align}\n", "$$\n", "and combine the terms to a single exponential, we find\n", "$$\n", "\\begin{align}\n", "U_{\\mu\\nu}(x) &= \\exp\\left(i a^2 g F_{\\mu\\nu}(x) + O(a^3)\\right)\n", "\\end{align}\n", "$$\n", "\n", "We can then consider the gauge-invariant trace for a gauge group SU$(N_c)$\n", "$$\n", "\\begin{align}\n", "\\textrm{Tr}\\: U_{\\mu\\nu}(x) = N_c + i a^2 g \\textrm{Tr}\\: F_{\\mu\\nu}(x) - \\frac12 a^4 g^2 \\textrm{Tr}\\: (F_{\\mu\\nu}(x))^2 + O(a^5) \n", "\\end{align}\n", "$$\n", "and\n", "$$\n", "\\begin{align}\n", "\\textrm{Tr}\\: \\left( \\mathbb{1} - U_{\\mu\\nu}(x) \\right) = -i a^2 g \\textrm{Tr}\\: F_{\\mu\\nu}(x) + \\frac12 a^4 g^2 \\textrm{Tr}\\: (F_{\\mu\\nu}(x))^2 + O(a^5) \\,.\n", "\\end{align}\n", "$$\n", "Next we take the real part to create a real-valued contribution to the action and use that\n", "$$\n", "\\begin{align}\n", "F_{\\mu\\nu}^\\dagger = F_{\\mu\\nu}\n", "\\end{align}\n", "$$\n", "such that\n", "$$\n", "\\begin{align}\n", "\\textrm{Re}\\: \\textrm{Tr}\\: \\left( \\mathbb{1} - U_{\\mu\\nu}(x) \\right) = \\frac12 a^4 g^2 \\:\\textrm{Tr}\\: (F_{\\mu\\nu}(x))^2 + O(a^5) \\,,\n", "\\end{align}\n", "$$\n", "where no sum over repeated indices is implied. In order to recover the continuum action at leading order,\n", "we therefore need to divide by $g^2$ and sum over $x\\in \\Lambda$, $\\mu$ and $\\nu$. This yields the **Wilson gauge action**\n", "$$\n", "\\begin{align}\n", "S_W &= \\frac{1}{g^2}\\sum_x \\sum_{\\mu,\\nu=0}^{d-1} \\textrm{Re}\\: \\textrm{Tr}\\: \\left( \\mathbb{1} - U_{\\mu\\nu}(x) \\right) = \\frac12 a^4 \\sum_x \\sum_{\\mu,\\nu=0}^{d-1}\\textrm{Tr}\\: (F_{\\mu\\nu}(x))^2 + O(a^5) \\,,\n", "\\end{align}\n", "$$\n", "where $d$ is the space-time dimensionality of the theory. For gauge group SU$(N_c)$ in the fundamental representation, we customarily define\n", "$$\n", "\\begin{align}\n", "S_W &= \\frac12 \\beta \\sum_x \\sum_{\\mu,\\nu=0}^{d-1} ( 1 - P_{\\mu\\nu}(x)) \n", "= \\beta \\sum_x \\sum_{\\mu<\\nu} ( 1 - P_{\\mu\\nu}(x)) \n", "\\end{align}\n", "$$\n", "with scalar plaquette\n", "$$\n", "\\begin{align}\n", "P_{\\mu\\nu}(x) &= \\frac1{N_c} \\textrm{Re}\\: \\textrm{Tr}\\: U_{\\mu\\nu}(x)\n", "\\end{align}\n", "$$\n", "and\n", "$$\n", "\\begin{align}\n", "\\beta &= \\frac{2 N_c}{g^2} \\,.\n", "\\end{align}\n", "$$\n", "Note that we used\n", "$$\n", "\\begin{align}\n", "U_{\\mu\\nu}(x)^\\dagger &= U_\\nu(x) U_\\mu(x + a \\hat{\\nu})U_\\nu(x+a \\hat{\\mu})^\\dagger U_\\mu(x)^\\dagger = U_{\\nu\\mu}(x)\n", "\\end{align}\n", "$$\n", "to replace the sum over $\\mu,\\nu$ with a sum over $\\mu<\\nu$." ] }, { "cell_type": "markdown", "id": "d236f461", "metadata": {}, "source": [ "## Improved gauge theory\n", "\n", "We can show that also other Wilson loops locally look like the field-strength tensor $F_{\\mu\\nu}$, albeit with different discretization errors. By taking appropriate combinations of Wilson loops, we can create an improved gauge action which reduces the discretization errors to a higher power in $a$. At this point, we could already determine the classically optimal combination of, e.g., the $1x1$ and $2x1$ rectangle terms, however, we will later see that quantum effects can disturb this cancellation. We will revisit this point at a later time." ] }, { "cell_type": "markdown", "id": "cad8d212", "metadata": {}, "source": [ "## Path integral measure" ] }, { "cell_type": "markdown", "id": "1904c84d", "metadata": {}, "source": [ "One can show that for finite groups $G = \\{ g_1, \\ldots, g_n \\}$ that for any element $g_i \\in G$,\n", "the sequence\n", "$g_i g_0,g_i g_1,\\ldots,g_i g_n$\n", "and the sequence\n", "$g_0g_i ,g_1g_i ,\\ldots, g_n g_i$\n", "are just reorderings of the sequence $g_0,\\ldots,g_n$. Therefore averaging a function $f$ over all elements\n", "of a finite group satisfies\n", "$$\n", "\\begin{align}\n", " \\frac1{n} \\sum_g f(g) = \\frac1n \\sum_g f(g_i g) = \\frac1n \\sum_g f(g g_i) \\,.\n", "\\end{align}\n", "$$\n", "For continuous Lie groups that we consider in this chapter as gauge groups\n", "one can define the **Haar integration measure** $d[U]$ that satisfies\n", "$$\n", "\\begin{align}\n", " \\int d[U] f(U) = \\int d[U] f(V U) = \\int d[U] f(U V)\n", "\\end{align}\n", "$$\n", "for group elements $U$ and $V$. Equivalently\n", "$$\n", "\\begin{align}\n", "d[VU] = d[U] = d[UV]\n", "\\end{align}\n", "$$\n", "for each group element $V$. Using such an integration measure, the integral over link variables $U_\\mu(x)$\n", "then also is gauge invariant since\n", "$$\n", "\\begin{align}\n", "d[V(x) U_\\mu(x) V(x + a\\hat{\\mu})^\\dagger] = d[U_\\mu(x)] \\,.\n", "\\end{align}\n", "$$\n", "If the action and integral measure is explicitly invariant under our gauge transformation, then the quantized theory as defined by the path integral is also explicitly gauge invariant. This is crucial since the measure may in principle violate symmetries of the action or in other words, a quantized theory may exhibit lower symmetries than the action itself has.\n", "\n", "In a later discussion we will give a procedure to derive the Haar measure for a given parametrization of a group manifold. For now, we just observe that in a stochastic sampling of the integral all we have to follow is to pick random elements $V$ of the gauge group and to consistently multiply them from the left or right to the existing links $U_\\mu(x)$ to update them. As long as we make sure that any group element can be reached by a product of these random $V$, the construction explicitly follows the Haar integration measure as defined by its invariance properties. For completeness, we note that the Haar measure is only unique up to a constant factor that can be fixed, e.g., by requiring $\\int d[U]=1$." ] }, { "cell_type": "markdown", "id": "23bc8c1a", "metadata": {}, "source": [ "## Markov chains for pure gauge theories" ] }, { "cell_type": "markdown", "id": "2dec2888", "metadata": {}, "source": [ "For pure gauge theories, there are various competitive algorithms to generate a Markov chain:\n", "- The first is a **local metropolis** using a checkerboarding strategy as we have done before for scalar field theory.\n", "- Another is to explicitly parametrize the Haar measure for SU(2) subgroups of the SU$(N_c)$ gauge group and to draw such SU(2) matrices directly from the correct probability distribution. This is the so-called SU(2)-subgroup **heatbath** algorithm. We will discuss it in detail in the next lecture. \n", "- We could also integrate a stochastic differential equation, called a **Langevin** integrator, and show explicitly that the correct stationary distribution is found. We will also discuss this method in the next lecture. \n", "- Finally, one can also use a combination of molecular dynamics evolution with random momentum sampling, the **hybrid monte carlo** method. We will dedicate an entire chapter to this method due to its importance in state-of-the-art calculations.\n", "\n", "In the following, we first demonstrate the ensemble generation with GPT using some of these algorithms.\n", "\n", "First using the local Metropolis:" ] }, { "cell_type": "code", "execution_count": 3, "id": "2741a984", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "GPT : 113.770477 s : Initializing gpt.random(test,vectorized_ranlux24_24_64) took 4.60148e-05 s\n", "GPT : 113.785935 s : Local metropolis 0 has P = 1.0\n", "GPT : 118.091342 s : Local metropolis 50 has P = 0.6232448286480373\n", "GPT : 122.690226 s : Local metropolis 100 has P = 0.558342350853814\n", "GPT : 126.951510 s : Local metropolis 150 has P = 0.5412292811605666\n", "GPT : 131.852604 s : Local metropolis 200 has P = 0.5315725008646647\n", "GPT : 136.840365 s : Local metropolis 250 has P = 0.5241931941774156\n", "GPT : 142.993895 s : Local metropolis 300 has P = 0.5202440818150839\n", "GPT : 147.940617 s : Local metropolis 350 has P = 0.5177053146892124\n", "GPT : 152.788150 s : Local metropolis 400 has P = 0.5131866865687901\n", "GPT : 157.826834 s : Local metropolis 450 has P = 0.5104659663306342\n", "GPT : 162.208752 s : Local metropolis 500 has P = 0.5131155252456665\n", "GPT : 166.561504 s : Local metropolis 550 has P = 0.5085434781180488\n", "GPT : 171.248844 s : Local metropolis 600 has P = 0.5039165417353312\n", "GPT : 175.708371 s : Local metropolis 650 has P = 0.506249725818634\n", "GPT : 180.486332 s : Local metropolis 700 has P = 0.5023225479655796\n", "GPT : 185.004515 s : Local metropolis 750 has P = 0.5003030101458231\n", "GPT : 189.150137 s : Local metropolis 800 has P = 0.4997691445880466\n", "GPT : 193.177941 s : Local metropolis 850 has P = 0.5011237727271186\n", "GPT : 197.519407 s : Local metropolis 900 has P = 0.5038102865219116\n", "GPT : 201.544591 s : Local metropolis 950 has P = 0.5002213915189108\n" ] } ], "source": [ "import gpt as g;\n", "\n", "L = [8, 8, 8, 8]\n", "grid = g.grid(L, g.single)\n", "grid_eo = g.grid(L, g.single, g.redblack)\n", "\n", "rng = g.random(\"test\", \"vectorized_ranlux24_24_64\")\n", "U = g.qcd.gauge.unit(grid)\n", "Nd = len(U)\n", "\n", "# red/black mask\n", "mask_rb = g.complex(grid_eo)\n", "mask_rb[:] = 1\n", "\n", "# full mask\n", "mask = g.complex(grid)\n", "\n", "# simple plaquette action\n", "def staple(U, mu):\n", " st = g.lattice(U[0])\n", " st[:] = 0\n", " Nd = len(U)\n", " for nu in range(Nd):\n", " if mu != nu:\n", " st += g.qcd.gauge.staple(U, mu, nu) / U[0].otype.Nc\n", " return st\n", "\n", "# g.default.push_verbose(\"local_metropolis\", True) # activaties display of acceptance rate\n", "markov = g.algorithms.markov.local_metropolis(rng, step_size=0.5)\n", "beta = 5.5\n", "\n", "plaquette_local_metropolis = []\n", "for it in range(1000):\n", " plaq = g.qcd.gauge.plaquette(U)\n", " plaquette_local_metropolis.append(plaq)\n", " if it % 50 == 0:\n", " g.message(f\"Local metropolis {it} has P = {plaq}\")\n", " for cb in [g.even, g.odd]:\n", " mask[:] = 0\n", " mask_rb.checkerboard(cb)\n", " g.set_checkerboard(mask, mask_rb)\n", "\n", " for mu in range(Nd):\n", " st = g.eval(beta * staple(U, mu))\n", " markov(U[mu], st, mask)" ] }, { "cell_type": "markdown", "id": "30c6fd98", "metadata": {}, "source": [ "Next we consider the Heatbath algorithm." ] }, { "cell_type": "code", "execution_count": 7, "id": "e57ac072", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "GPT : 300.293935 s : SU(2)-subgroup heatbath 0 has P = 1.0\n", "GPT : 313.256864 s : SU(2)-subgroup heatbath 50 has P = 0.5022712548573812\n", "GPT : 326.851409 s : SU(2)-subgroup heatbath 100 has P = 0.4951641890737746\n", "GPT : 341.379575 s : SU(2)-subgroup heatbath 150 has P = 0.49848774406645036\n", "GPT : 356.706945 s : SU(2)-subgroup heatbath 200 has P = 0.5004462268617418\n", "GPT : 373.774078 s : SU(2)-subgroup heatbath 250 has P = 0.49490469694137573\n", "GPT : 389.923751 s : SU(2)-subgroup heatbath 300 has P = 0.49955056111017865\n", "GPT : 405.860566 s : SU(2)-subgroup heatbath 350 has P = 0.49755699104732937\n", "GPT : 420.428378 s : SU(2)-subgroup heatbath 400 has P = 0.499170528517829\n", "GPT : 435.301122 s : SU(2)-subgroup heatbath 450 has P = 0.5017324553595649\n" ] } ], "source": [ "g.default.push_verbose(\"su2_heat_bath\", False)\n", "markov = g.algorithms.markov.su2_heat_bath(rng)\n", "U = g.qcd.gauge.unit(grid)\n", "plaquette_heatbath = []\n", "for it in range(500):\n", " plaq = g.qcd.gauge.plaquette(U)\n", " plaquette_heatbath.append(plaq)\n", " if it % 50 == 0:\n", " g.message(f\"SU(2)-subgroup heatbath {it} has P = {plaq}\")\n", " for cb in [g.even, g.odd]:\n", " mask[:] = 0\n", " mask_rb.checkerboard(cb)\n", " g.set_checkerboard(mask, mask_rb)\n", "\n", " for mu in range(Nd):\n", " st = g.eval(beta * staple(U, mu))\n", " markov(U[mu], st, mask)" ] }, { "cell_type": "markdown", "id": "b1ddf74a", "metadata": {}, "source": [ "Finally, we demonstrate the Langevin algorithm." ] }, { "cell_type": "code", "execution_count": 9, "id": "47ca73bb", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "GPT : 511.601426 s : Langevin_bf(eps=0.02) 0 has P = 0.913337786992391\n", "GPT : 514.370879 s : Langevin_bf(eps=0.02) 50 has P = 0.566335055563185\n", "GPT : 517.300690 s : Langevin_bf(eps=0.02) 100 has P = 0.5438746478822496\n", "GPT : 520.499172 s : Langevin_bf(eps=0.02) 150 has P = 0.5410566197501289\n", "GPT : 523.587544 s : Langevin_bf(eps=0.02) 200 has P = 0.533349891503652\n", "GPT : 527.598893 s : Langevin_bf(eps=0.02) 250 has P = 0.5256600843535529\n", "GPT : 532.062941 s : Langevin_bf(eps=0.02) 300 has P = 0.5235162112447951\n", "GPT : 535.788070 s : Langevin_bf(eps=0.02) 350 has P = 0.5187506278355917\n", "GPT : 541.385091 s : Langevin_bf(eps=0.02) 400 has P = 0.5197881857554117\n", "GPT : 547.443904 s : Langevin_bf(eps=0.02) 450 has P = 0.5210425655047098\n", "GPT : 552.004728 s : Langevin_bf(eps=0.02) 500 has P = 0.5160737898614671\n", "GPT : 555.507547 s : Langevin_bf(eps=0.02) 550 has P = 0.5104968547821045\n", "GPT : 559.059208 s : Langevin_bf(eps=0.02) 600 has P = 0.5076582895384895\n", "GPT : 562.807271 s : Langevin_bf(eps=0.02) 650 has P = 0.5071470008956062\n", "GPT : 566.821152 s : Langevin_bf(eps=0.02) 700 has P = 0.5031868351830376\n", "GPT : 570.149039 s : Langevin_bf(eps=0.02) 750 has P = 0.5077733132574294\n", "GPT : 573.409066 s : Langevin_bf(eps=0.02) 800 has P = 0.498729235596127\n", "GPT : 576.775995 s : Langevin_bf(eps=0.02) 850 has P = 0.5024757583936056\n", "GPT : 579.928553 s : Langevin_bf(eps=0.02) 900 has P = 0.5047208666801453\n", "GPT : 583.567940 s : Langevin_bf(eps=0.02) 950 has P = 0.5034748845630223\n" ] } ], "source": [ "U = g.qcd.gauge.unit(grid)\n", "w = g.qcd.gauge.action.wilson(beta)\n", "l = g.algorithms.markov.langevin_bf(rng, epsilon=0.02)\n", "plaquette_langevin = []\n", "for it in range(1000):\n", " l(U, w)\n", " plaq = g.qcd.gauge.plaquette(U)\n", " plaquette_langevin.append(plaq)\n", " if it % 50 == 0:\n", " g.message(f\"Langevin_bf(eps=0.02) {it} has P = {plaq}\")" ] }, { "cell_type": "code", "execution_count": 16, "id": "97dc4acf", "metadata": { "scrolled": true }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAA4mklEQVR4nO3de3Qcdfn48ffTCw03k0IL1hZ2CxQqbS69UMq3QhbKpUBpQVCoiE1REQSh8BUsfoHdgByreKBfLv5QvkIQ0SIXS0GPSLUbQSm2tRGkQCmwgUDFWJrQW3p9fn/MzmZ2s7vZTTa3yfM6Z0+yc9n5zE7yzGc+n888I6qKMcYY/xrQ0wUwxhjTtSzQG2OMz1mgN8YYn7NAb4wxPmeB3hhjfM4CvTHG+FxOgV5EZojImyKyXkQWpJl/l4jUxV/rRKTJM2+uiLwVf80tYNmNMcbkQNobRy8iA4F1wGlAA7ASmKOqazMs/y1ggqpeKiIHAauAyYACq4FJqrqpcLtgjDEmm1xq9FOA9ar6jqruBBYDs7MsPwf4Vfz3M4DnVfXjeHB/HpjRmQIbY4zJz6AclhkJvO953wAcn25BEQkAo4E/ZVl3ZJr1LgMuA9h///0njR07NodiZfbh5g/5zIGfyTh/9WqYNKlTmzDGmF5l9erV/1HV4enm5RLo83ER8ISq7slnJVX9KfBTgMmTJ+uqVasKXKxkIjBzJkQiXboZY4zpNiJSn2leLk03HwCHed6Pik9L5yJam23yXbfbhMNQXd3TpTDGmO6RS6BfCYwRkdEisg9OMF+aupCIjAWGAi95Jj8HnC4iQ0VkKHB6fFqXi0QjmedlnmWMMb7TbqBX1d3AVTgB+nXg16r6mojcKiKzPIteBCxWzzAeVf0YuA3nZLESuDU+rctV17atskciTrONiPPe/d0CvzHGz9odXtndCtVGL9WChjPvmwj0sl03cbt27aKhoYGWlpaeLooxvU5RURGjRo1i8ODBSdNFZLWqTk63TqE7Y3tUJBpJqslLtVN1D1eGiYQiPVQqk6+GhgYOPPBAgsEg4l5+GWNQVTZu3EhDQwOjR4/OeT1/BfpQJBHQ26vRm96rpaXFgrwxaYgIBx98MI2NjXmt169y3Vgbfd9hQd6Y9Dryv+GrGr1XuDLcZlok0hrUrY3eGNNf+LZGb23ypjMOOOCAPvXZdXV1/O53vyv456aKRqPMnDkTgKVLl7Jw4cIu36bpPN8GetP/9OcmuGyBfvfu3V2yzVmzZrFgQZtktqYX8m2gT3fDlLXR+1tX3+1cV1fH1KlTKSsr47zzzmPTJicJ6/r16zn11FMpLy9n4sSJvP3222zZsoXp06czceJESktLefrpp7N+diwWY+zYsVRVVXH00Udz8cUXs2zZMqZNm8aYMWP429/+BsDWrVu59NJLmTJlChMmTODpp59m586d3HLLLTz22GNUVFTw2GOPEYlEuOSSS5g2bRqXXHIJsViMU045hbKyMqZPn857770HQFVVFZdffjmTJ0/m6KOP5tlnnwWcDvF58+ZRWlrKhAkTWL58eZsy19TUcNVVVwHw+OOPM378eMrLyznppJMK9p2bAlHVXvWaNGmSFgIRss9HNRwuyKZMga1du7ZD65H9kOdl//33bzOttLRUo9GoqqrefPPNes0116iq6pQpU/Spp55SVdXt27fr1q1bddeuXdrc3Kyqqo2NjXrkkUfq3r17M372u+++qwMHDtRXXnlF9+zZoxMnTtR58+bp3r17dcmSJTp79mxVVb3xxhv1kUceUVXVTZs26ZgxY3TLli360EMP6ZVXXpn4vHA4rBMnTtRt27apqurMmTO1pqZGVVV/9rOfJT5v7ty5esYZZ+iePXt03bp1OnLkSN2+fbv+6Ec/0nnz5qmq6uuvv66HHXaYbt++XZcvX65nn322qmrSNsePH68NDQ2Jcpmule5/BFilGeKqb2v0ubB8N31fd12lNTc309TURGVlJQBz587lz3/+M5s3b+aDDz7gvPPOA5ybWfbbbz9Ule9+97uUlZVx6qmn8sEHH/DRRx9l3cbo0aMpLS1lwIABjBs3junTpyMilJaWEovFAPjDH/7AwoULqaioIBQK0dLSkqidp5o1axb77rsvAC+99BJf+tKXALjkkkt48cUXE8t98YtfZMCAAYwZM4YjjjiCN954gxdffJEvf/nLAIwdO5ZAIMC6desyln3atGlUVVXxwAMPsGdPXjkNTTfwVaCPRCNItSRulHJ/T9eME247KMf0QZGIM3rKHUHl/t7TzXGPPvoojY2NrF69mrq6Og499NB27/QdMmRI4vcBAwYk3g8YMCDRzq6qPPnkk9TV1VFXV8d7773HZz/72bSft//+++dU1tTheh0Zvnf//ffzve99j/fff59JkyaxcePGvD/DdB1/BfpQBA1r4kYp93fvCBy3BujW5q2d3uSiuLiYoUOH8sILLwDwyCOPUFlZyYEHHsioUaNYsmQJADt27GDbtm00NzdzyCGHMHjwYJYvX059fcYMsnk544wzuOeee9D4mW3NmjUAHHjggWzevDnjev/1X//F4sWLAeckdOKJJybmPf744+zdu5e3336bd955h2OOOYYTTzyRRx99FIB169bx3nvvccwxx2T8/Lfffpvjjz+eW2+9leHDh/P+++9nXNZ0P9+Oo8/ExtL7VyGv0rZt28aoUaMS76+77joefvhhLr/8crZt28YRRxzBQw89BDhB/xvf+Aa33HILgwcP5vHHH+fiiy/mnHPOobS0lMmTJ9PZh+m4br75ZubPn09ZWRl79+5l9OjRPPvss5x88smJJp0bb7yxzXr33HMP8+bN44477mD48OGJsgMcfvjhTJkyhU8++YT777+foqIivvnNb3LFFVdQWlrKoEGDqKmpSbriSHX99dfz1ltvoapMnz6d8vLyguyvKQzfJjWLRCPtjqW3QN87vf766xmbI0xhVVVVMXPmTC644IKeLorJQ7r/kWxJzXzVdJMLG2JpjOlvfFujzyWpmdXoeyer0RuTndXorWpujDFJfBXoI9EIItVZh1da040xpr/xX9NNvD0m16abcNiCfG9jTTfGZNc/m27SVdPd6e2wu2ONMX7nn0CfcntkuLL9qrrdHWuM6Q/8EejTiIQiaVMfgN0d61sFPIC3334748aNo6ysjIqKCl5++WUAgsEg//nPfxLLefOzAyxZsoRbb70VgDvvvJNjjz02kTHSvTu2sbGRGTNmpN1uLBZj/PjxKbsV4Uc/+lGH9iMajfLXv/418b6qqoonnngi5/VjsRi//OUvE++9GStN3+G/QO+ppnsfFO7VW/OjmE4qUDvcSy+9xLPPPsvf//53XnnlFZYtW8Zhhx2W07o//OEP+eY3vwnAhAkTWLVqFa+88goXXHABN9xwAwDDhw9nxIgR/OUvfylIebNJDfT5Sg30pm/yX6C3aG06acOGDQwbNixxy/+wYcP4zGc+0+5669atY8iQIQwbNgyAk08+mf322w+AqVOn0tDQkFj23HPPTeSSycfbb7/NjBkzmDRpEieeeCJvvPEGAM888wzHH388EyZM4NRTT+Wjjz4iFotx//33c9ddd1FRUZHI07Ns2bI2+edjsRgnnngiEydOZOLEiYmTw4IFC3jhhReoqKjgrrvuAuDDDz9kxowZjBkzJnHyMr1cpvzFPfXqbD768PKwEqHNK7w83LpM2K3DJ78sP33vkFc++i44mJs3b9by8nIdM2aMXnHFFYkc9KqqgUBAGxsbE++9+dkffPBBve6669J+5pVXXqm33XZb4n1DQ4OOHz++zXLvvvuuFhUVaXl5eeJ16KGH6h133KGqqqeccoquW7dOVVVXrFihJ598sqqqfvzxx4l89w888ECiHOFwOLGuaub881u3btXt27erquq6devU/T/07p+qk4N+9OjR2tTUpNu3b9fDDz9c33vvvZy+V1M4+eaj911Ss0ioNcdNpiGWltjMR7rgYB5wwAGsXr2aF154geXLl3PhhReycOFCqqqq0qbwdadt2LCB4cOHt5n/i1/8glWrVlFbW5uYdsghh/Dhhx+m3f6RRx5JXV1d4n0kvn9btmzhr3/9K1/4whcS83bs2AFAQ0MDF154IRs2bGDnzp2MHj064/6lyz8/evRorrrqKurq6hg4cGDW3PPTp0+nuLgYgGOPPZb6+vqcm7ZMz/Bf000O7KYp056BAwcSCoWorq7m3nvv5cknnwTg4IMPTjxCEODjjz9ONNXsu+++bXLOL1u2jNtvv52lS5cmZX9saWlJPBQkV3v37qWkpCSRi76uro7XX38dgG9961tcddVVvPrqq/zkJz/Jmvs+Xf75u+66i0MPPZR//OMfrFq1ip07d2Zc37sfAwcO7LJn0prC6beB3jpjfahA42XffPNN3nrrrcT7uro6AoEAAKFQiEceeQSAPXv28Itf/IKTTz4ZgM9+9rOsX78+sd6aNWv4xje+wdKlSznkkEOStrFu3bo2o2va86lPfYrRo0fz+OOPA06z6z/+8Q/AeQLWyJEjAXj44YcT66TLU58u/3xzczMjRoxgwIABPPLII4mnRLWX5970DTkFehGZISJvish6EUn72HcR+aKIrBWR10Tkl57pe0SkLv5aWqiCG9NGgc7UW7ZsYe7cuYmhkWvXrk00n9x8882sX7+e8vJyJkyYwFFHHZV45N5JJ53EmjVrEg8Fuf7669myZQtf+MIXqKioYNasWYltLF++nLPPPjvvsj366KP87Gc/o7y8nHHjxiUeOh6JRPjCF77ApEmTElcYAOeccw6/+c1vkjpj3fzzZ555ZlL++Ycffpjy8nLeeOONxNOpysrKGDhwIOXl5YnOWNP3tJsCQUQGAuuA04AGYCUwR1XXepYZA/waOEVVN4nIIar67/i8Lap6QK4F6mwKhEg0knZYZbgynDY/vbXR9z59OQXCNddcwznnnMOpp56adbmTTjqJp59+mqFDh3ZTyYyfdEUKhCnAelV9R1V3AouB2SnLfB24T1U3AbhBvifk8zhBa6M3hfbd736Xbdu2ZV2msbGR6667zoK86Ta5BPqRgPcBkA3xaV5HA0eLyF9EZIWIeG/7KxKRVfHp53auuIWR2kYfDlsbvSmMQw89NKmJJp3hw4dz7rnndk+BjKFwnbGDgDFACJgDPCAiJfF5gfjlxJeARSJyZOrKInJZ/GSwqrGxsUBFyp0lNjPG+Fkugf4DwDtIdlR8mlcDsFRVd6nquzht+mMAVPWD+M93gCgwIXUDqvpTVZ2sqpPTjUPuSpbYzBjjd7kE+pXAGBEZLSL7ABcBqaNnluDU5hGRYThNOe+IyFARGeKZPg1YSxeKRCOJB45A+oePgCU2M8b0H+0GelXdDVwFPAe8DvxaVV8TkVtFxG2MfA7YKCJrgeXA9aq6EfgssEpE/hGfvtA7WqcrpHbGhivDbTpjwcbSm+wOOCDngWLd4pZbbmHZsmV5rbNmzRq++tWvAk7WyeHDh1NRUcG4ceO44IILEp3GkUiEkSNHUlFRQUVFBQsWpB1BnbPe9t11NVXl6quv5qijjqKsrIy///3vaZdbvXo1paWlHHXUUVx99dVJw3DHjh1LWVkZ5513Hk1NTQC8+uqrVFVVFa6QvenV2Vw3Xm6em3aXa38R043yynXj4c1n1Fn7779/wT6rp1xwwQVaV1enqk6OmiuvvDIxb86cOfrggw+qatt8OJ3lh+8uH7/97W91xowZunfvXn3ppZd0ypQpaZc77rjj9KWXXtK9e/fqjBkz9He/+52qqj733HO6a9cuVVW94YYb9IYbbkisM336dK2vr2/zWfnmuvH1nbHhyswN8DbE0n8ypaUulHQZIsGpEV966aWEQiGOOOII7r777sQ6t912G8cccwyf+9znmDNnTiKvfLoslM3NzQQCAfbu3QvA1q1bOeyww9i1a1dSHvlgMEg4HGbixImUlpYmMlh6bd68mVdeeYXy8vI283bv3s3WrVvzGt559913J24gu+iiiwDnxrJ58+ZRWlpKWVlZIk0EwP/8z/9QXl7O1KlTE99TY2Mj559/PscddxzHHXdcIk1zJBJh7ty5nHjiiQQCAZ566iluuOEGSktLmTFjBrt27QKcGnFlZSWTJk3ijDPOYMOGDVnLvHXrVi699FKmTJnChAkTEjeX1dTUMHv2bEKhEGPGjKE63n67detWzj77bMrLyxk/fjyPPfZYTt/N008/zVe+8hVEhKlTp9LU1NSmbBs2bOCTTz5h6tSpiAhf+cpXWLJkCQCnn346gwY5acdSs5yec845LF68OKdyZJXpDNBTr0LU6HPJYJl8Juz0Jk0BdbRGn8vVW67S1UqzZYg84YQTtKWlRRsbG/Wggw7SnTt36t/+9jctLy/X7du36yeffKJHHXVUu1koZ82apX/6059UVXXx4sX61a9+VVWdrJOPP/64qjoZNO+++25VVb3vvvsSy3j96U9/0s9//vOJ9w899JAOGzZMy8vL9ZBDDtHPfe5zunv37kT5P/OZzySyZf7+979v83kjRozQlpYWVVXdtGmTqjq1z2uuuSbp+1FVBXTp0qWqqnr99dcnsnbOmTNHX3jhBVVVra+v17Fjxya2P23aNN25c6fW1dXpvvvum6jtnnvuufqb3/xGd+7cqSeccIL++9//Tnw38+bNa1NOrxtvvFEfeeSRRJnHjBmjW7Zs0Yceekg//elP63/+8x/dtm2bjhs3TleuXKlPPPGEfu1rX0us39TUpKqq8+fPT8om6r6+//3vq6rq2WefndgvVefYrly5MqksK1eu1OnTpyfe//nPf07KCuqaOXNmosyqqi+++KLOnDmzzXL9Pnsl5J7B0jus0q3Z28PC+5bUO6HdTvhMd0J3RrYMkWeffTZDhgxhyJAhHHLIIXz00Uf85S9/Yfbs2RQVFVFUVMQ555wDZM9CeeGFF/LYY49x8skns3jx4sRDTFJ9/vOfB2DSpEk89dRTbeany6R54YUXcu+996KqXHnlldxxxx2J9vhrr72Wb3/72xn3vaysjIsvvphzzz03cQ/AsmXLkmqb7hXCPvvsk3jq1qRJk3j++ecTy69d29pF98knn7BlyxYAzjzzTAYPHkxpaSl79uxJPIGrtLSUWCzGm2++yT//+U9OO+00wMkzNGLEiIzlBfjDH/7A0qVLE1dRLS0tvPfeewCcdtppHHzwwYnv8sUXX+Sss87iv//7v/nOd77DzJkzOfHEEwG6LfXD7bffzqBBg7j44osT07JlOc2HLwN9LlKz24KlQuiLcjmpF8q3vvUtrrvuOmbNmkU0Gk3kv4H8Mjp6s1CmmjVrFt/97nf5+OOPWb16Naecckraz3C3l2lb6TJpukSEc845h3vuuSfnjtff/va3/PnPf+aZZ57h9ttv59VXX8247ODBgxMZMr3l27t3LytWrKCoqCjj/gwYMCBp/QEDBrB7925UlXHjxvHSSy/lVF5wWiuefPJJjjnmmKTpL7/8ctoMnkcffTR///vf+d3vfsdNN93E9OnTueWWW7j22mtZvnx5m8+/6KKLWLBgASNHjuT991vvKW1oaEgkmHONHDkyqUkmdZmamhqeffZZ/vjHPyaVrSNZTtPxdRu9MYWUKUNkJtOmTeOZZ56hpaWFLVu2JJ7mlC0L5QEHHMBxxx3HNddcw8yZMxk4cGCHypqaSTPViy++yJFHtrl3Ma29e/fy/vvvc/LJJ/ODH/yA5uZmtmzZwmmnncZ9992XWM6bvjmd008/nXvuuSfxPt2JLpNjjjmGxsbGRKDftWsXr732GgD33nsv9957b5t1zjjjDO65557E6JY1a9Yk5j3//PN8/PHHbN++nSVLljBt2jQ+/PBD9ttvP7785S9z/fXXJ0bP3HXXXUmpod2Xe5KcNWsWP//5z1FVVqxYQXFxcZurjREjRvCpT32KFStWoKr8/Oc/Z/ZsJ5PM73//e374wx+ydOnSxBPJXB3JcpqOLwN9vmPpXdYh27dl63zP17Zt2xg1alTideedd2bMEJnJcccdx6xZsygrK+PMM8+ktLQ08cCOTFkowWli+cUvfsGFF17Y4fKPHTuW5ubmpBTDjz32GBUVFZSVlbFmzRpuvvnmnD5rz549fPnLX6a0tJQJEyZw9dVXU1JSwk033cSmTZsYP3485eXlaWu9XnfffTerVq2irKyMY489lvvvvz/n/dlnn3144okn+M53vkN5eTkVFRWJxx2+8cYbiWYYr5tvvpldu3ZRVlbGuHHjkvZ3ypQpnH/++ZSVlXH++eczefJkXn31VaZMmUJFRQXV1dXcdNNNOZXtrLPO4ogjjuCoo47i61//Oj/+8Y8T8yoqKhK///jHP+ZrX/saRx11FEceeSRnnnkmAFdddRWbN2/mtNNOo6KigssvvzyxTkeznLaRqfG+p16FHF6pmnsHnXXI9h4d7YztjTZv3qyqqlu3btVJkybp6tWru23bd955pz7wwAPdtr2ecvbZZ+uOHTtyXj51qGlv1dLSoscff3xi6KWXDa/Mgw2xNF3tsssuo6KigokTJ3L++eczceLEbtv2FVdckdR34FfPPvss++yzT08Xo+Dee+89Fi5cmBh62Rnt5qPvbp3NR58q1w46ERtx01v05Xz0xnSHrshH3+elts0npkeSa/TV1Vaj7y16WwXEmN6iI/8bvgz0qZ2x1bXVGTtjVZOHVarlu+lxRUVFbNy40YK9MSlUlY0bN6YdoppNv2i6ATI236TeOOWyZpyes2vXLhoaGjKOAzemPysqKmLUqFEMHjw4aXq2phvfBvp8nx3rBvxe9nUYY0xO+mWgd7kBP9cO2V72dRhjTE76dWdse/lObIilMcbvfB/o22MPCjfG+J1vA32+aRBsiKUxxq9830YP+d00BdZOb4zpe/p1G30uLLmZMcbPfBvoc226McYYv+s3TTe5PnHIhlgaY/qiftl0k2saBLAhlsYYf/PtowRTHzEHmdMgGGOMn/m2Rg+ttXpXtiGWltzMGONX/g70oQga1sQj5jSsaFgz5rqxkTfGGD/yd6CP1+jd5GY28sYY0x/1i1E3kP/IG7DRN8aYvqNfjrqBzo28AWu6Mcb4Q06BXkRmiMibIrJeRBZkWOaLIrJWRF4TkV96ps8Vkbfir7mFKngu3DZ672ibbG304XDbz6iutmBvjOnb2m26EZGBwDrgNKABWAnMUdW1nmXGAL8GTlHVTSJyiKr+W0QOAlYBkwEFVgOTVHVTpu11VT76VNmacazpxhjT13S26WYKsF5V31HVncBiYHbKMl8H7nMDuKr+Oz79DOB5Vf04Pu95YEZHdqI7pavZG2NMX5VLoB8JvO953xCf5nU0cLSI/EVEVojIjDzWRUQuE5FVIrKqsbEx99LnILX5Jpchlu4zZG2IpTHGDwrVGTsIGAOEgDnAAyJSkuvKqvpTVZ2sqpOHDx/euZJYVDbGmCS5BPoPgMM870fFp3k1AEtVdZeqvovTpj8mx3ULqzq5PT6fLJZ2h6wxxo9yCfQrgTEiMlpE9gEuApamLLMEpzaPiAzDacp5B3gOOF1EhorIUOD0+LRuk8/IG7A7ZI0x/tNuoFfV3cBVOAH6deDXqvqaiNwqIrPiiz0HbBSRtcBy4HpV3aiqHwO34ZwsVgK3xqcVVp7pJ+3OWGNMf+K/O2PTJJQP1YSora9ts2hloJJoVTTrR4XDVps3xvR+/fbOWFcoGMp5uj0s3BjjN/4L9GkGwUdCkUQGS6/q2mprxjHG+J7/An0nq96WCsEY4zf+C/QFEI3mN90YY3oz/3XG5iDXRwtGIq3D8nvZ12SMMUn6fWcstL1xCtq/ecp775WNpzfG9FW+fTh4qnQPCwd7YLgxxv/6TY0e2j4sHLBHCxpjfK9fBfpM0g2ztNE3xhi/6FeBPtN4+ozLR9p2wtqdssaYvqZfjrqB3NrpQyGobZs5gcpKG2ppjOldbNSNRz7t9NGoE9RT1dY6JwFjjOkL+l+gD0WoDLSN3tW11YRqQjl/Tm2tNeEYY/qGfjO8sqPcJhpvM4413Rhj+hJ/1ujbqWpHq6Jpa/W19bVpa/WpbfW1tc7NU9Z8Y4zpC/wZ6FMeJ5hOPqmLo1Fn9I23vd5q9caYvsKfgT4H+bbVW63eGNNX+SfQ5/k4QXCacNINrUz31Kl0tfpw2Gr1xpjez1+BXrX1Dif39yyBPt1QS0g/3NI9j3hr9e7Tp6xWb4zpzfwT6DsgU/MNOE04JQtLWpeNpB9Tb4wxvZ0/A326JDUZRKuiGdMiNO9oTl42mv6jbUy9MaY382egzzPqZqvZe2v12SxalNcmjTGm2/gz0BdQ847mnFIZNzdbrd4Y0zv5M9B3IOJmuonKtWjFosRHZ2qrr6nJe7PGGNPl/Jm9UqRTD3lNNxIHIFAcIDY/BkAwCPX16dcPBCAW6/DmjTEmb5a9Mk+Zavb1zfWJm6liscw1+2DQqfnbsEtjTG/gn0DfgRumMolWRQkUB9LO8+bDiUad2nubZWqdMfaWztgY0xvkFOhFZIaIvCki60VkQZr5VSLSKCJ18dfXPPP2eKYvLWThk3TghqlsYvNjGWv2tfW1lCwsIRKNUFWV/XPcYG8dtcaYntJuG72IDATWAacBDcBKYI6qrvUsUwVMVtWr0qy/RVUPyLVAvaGNvs3HZWizBwir5pJDDXBq/1VVFvSNMYXX2Tb6KcB6VX1HVXcCi4HZhSxgbxZcFMy+QCiCavomnFT19Tkl1jTGmILKJdCPBN73vG+IT0t1voi8IiJPiMhhnulFIrJKRFaIyLnpNiAil8WXWdXY2Jhz4ZMUsI3eq6qiKuv86tpqpFoIRkJQFaJ4dvvbKynpVJGMMSYvheqMfQYIqmoZ8DzwsGdeIH458SVgkYgcmbqyqv5UVSer6uThw4cXqEiFEQlF0LBmfIC4q7a+FoK1NE+ohqpQ1mWbmwtyDjLGmJzkEug/ALw19FHxaQmqulFVd8Tf/h8wyTPvg/jPd4AoMKET5c2swJ2xnRKspfiaEIG5GbYdikAowsIVEYLBtrPtBGCMKaRcnhm7EhgjIqNxAvxFOLXzBBEZoaob4m9nAa/Hpw8FtqnqDhEZBkwDfliowmcVChU8WXy4Mkw0FqXuX3VtEp6lah5aS/PQWorDNTTvaIKFTc6MqhAEnVzHO4D6HRAMOgHfvcnKvRHLAr4xphByujNWRM4CFgEDgQdV9XYRuRVYpapLReT7OAF+N/AxcIWqviEi/wX8BNiLc/WwSFV/lm1bBRl1E4k4vZ5ddNdvJBqhujbPXtVYJcRCEEqznjsvJZ+OPa7QGJOrbKNu/JkCAQo+xDKd4KIg9c0Z8iB0RCRe3lAkEfSLi6GiwgK+MSa7/pMCIRRKP/Kmi25PdfPeFEwo4jTteGr9zc3pb7pK91xbY4xJx1+B3n2wa2qHbBdWh8OV4US6hHBlOGsGzHaFqhPt90kjd0KRRFoFgGAklJSKwRhjsrGmmwIJ1YQSDxUP1TiBuCCiYecEEA1DRY0zraS1uagyUEm0KkokGiESihAKOZ26lj3TmP6l/7bRh8M9NnQlEo0QjUWJNcVoamlKGqXj1vpzGb1DUyApsKdTPKSY5h3NVGqY2uoIYKmSjelv+k+gd0fbpOrBgO9y8+W4NXBX0feK2LFnR4a1OiBWCZ+ug5YSKtfECEWcmr4xxt/6T2ds6k1T0CuCPLS233uDPMCCz7VJBto5wVooaoaSemonBKmurW4/X48xxtf8FeihNeeNq7q6V+QbiIQibYK8O13D2rlO3EziTT71zfXtPvPWGONf/mq68fI24/SyfcxEqiUR8AvWmesxcFcxg1bNZ8ECIBrp6XOfMaaAsjXd5JICoW8JhZyB515uDb+XNONkEq4MJ9rTvXffVgYqE4F/4OYAB6yvcpKn5WnP4Gb2nFBNdS3QXAPRKgCisSihWJQaQlSFQtTU1RT+HgFjTI/xb40eWgN8L9vHXLl33mpYCdWECAVDiRNBsCpC/cOReIK0AiS5Txnd415ZhIIhorFo0tBRbzncYZ3ZeIeeGmO6Rv8ZdZPKDfR9OGlMtiCZaJ2KSOt4e5zhllv+U8KeAwuXnqF4SDElRSWJlA+B4gDBkiC19bVpUziHakLEmmLE5seQaiFQHLCrBGO6UP8N9N5mnF62n4USiQChCDVVEYJVzg1TkZDT/l7dHGx3DH4hVAYqiTXFgNa0EO5w0kBxIOnkUFVR1e4VQCQaoaauJmnZkoUllBSVJJbJdtLI5SrDGL/pf4G+F4+n727BqoiT9rgQzTs5GLI9wO5BTewZnPlGMO8VgBuU3VFB0Vg0qSPaTS+Rmjwu3We4pFqS+jtc3pQR3qYoa1YyfmCB3qsPN+N0lAhO7pySWLfU8NszcHOAmw6IEQ06qSICDyn18zI/gD2dykAloWAIcB7nqGFNnCzcTmwNK8FFQWLzY21SS3tPIOlOCsb0Nf0v0Hv18Q7ZQohEoKYGqqpgEUGa64PZg34OaRcKqgDbc9NA5DvPFa4MA06zl1vzT+2Aducb0xtZoId+2WyTifvwLZkXgpIYxe9W0TK2hh371jspFGqiTgdvP+TtUygeUsz8qfOTrgSs9m96q/4Z6NONp4d+2XSTiYjzYJOmJud9JBqBaMRp9QpFIBh1ZsRCMHURtJRAU7A1lXJ73Jp6d18hdDG3U9kr3XBTN7Gd9QGY7tA/A72XNd+kFYmkv8iJRJxzYbrzZOtC8e/UM6wzKaBHw85TsuYHoa4q/87glmInZ4/n/cAhLeyRAiaAKyD3XgfvcFN35JH7rOGMw2TthGAKoP8Geht90ylZA757o1ZEIRRhyBDY8VzECexFTa0PQyf+OMT5EWolx2DvbT5qcjpNWRQDnPTLQNbOW2kOoMXdewXhvXvZHW6aOlIoU7OPe0JIdz+CMbnqv4He1Q/G03elYNDpyE09ZxbPjtD8dARwvlb33BmNOl95dXVyS1kkAtWSIUDHKp0mIfdKAJKenduGtw8htWkooq3zU+fFKgnPDbHw9zXskKbkq4Zu4AZztxYfCoYSfQDFQ4qp+HQFsaYYVRVV7V4FWF+B8bJAD8mdsmA1+g5wn15VX996URQKOfPSdXukaxoKVkWoH13NkO0BdrxcBcEogfIY9dfGsgf2VFWh1t9ros6VRFOw9X0o4jyRa1EMFpTAvyqcvobUzy9UCok85DIKyOUOEQXnJjH3jmM3NUYqOwH0X/030GfqkIXu7ZTN1BjeR3V2d9xgFAw6VwvRaGszkXviqKlx5q1YATu6ulneG+ybAk6fAjgd0EXNrf0QPdCp7B0FpGFNNPO487x3D7t9BG66Cfe+glwDf2dPEnaS6Vn9N9BD9mDfXW313fz8Wr/J1NpTUN4+B6+qUOsVQjTSozeeZboScB9o4z0JePsMwpVhFkVrKCkhYwoK94aycGU4YwrrTIHcvc/AXT8aiybm1f2rjqYFTfntaIFFopFEJSK1/MFFQYIlwaREfR29WzrfE2uh9e9AD9nvlO2OYJ8a6LMNd+lNNf9eUh63CKmH0O2YbWqC5mbnfX2a+Ftc7MxvVz5NRxFx+hWg19xxnKshGyopKoL5JVHnEKfcNUysknAg2ubQu8npqiqqErmIaupq2nQ6p3KvMHoq3YT3BKhhJRh0miAHfTuYlPivWAOUlCTfLe2e3LxlT7cf3u+wpzrVLdBnC/TQNcE+24if6urW3kvvdntbzb+XlcftIwCnc9j71bn/vMGgE/grKpz37nJuM1HqxV1xceuy6eZ7hcNOk5KTOyiSfFJwM4gGo06nstsEFIw6z/DN1unrLtvNfQVEw1TOjaZ/yE2sksDyqPNrrLVZqKPcKww362khA362JqOShSVJV0FDXgqz47mI82/YzqXikO0Bduxbn2g+c9NsVNdWJ66i3AR83pOde1XjZm/tLhbo2wv0UNhg7wbw1Kqo+12LJAd8lzu9o+UoZA3c/c562d9HZ3n/FNJ103hPJt6rg0CgdbqbUsLVJvC7zT1tNp4hsLjNRQtKnBNCb7nBLBofuNAFJ6CAVkJJ5wOhG3i9/RXpgm8hePtLILl5LBN3qG26/Uz0VaVpPuoIC/SQva0+EHCqftA2UHYkeLo14VwalzO1K4TDThRyeyrTlctbvkikMDXwQtx70JVNUwX4jBVFIVqmhgi100zjBv3Uq4d0y7nDSbN/YCQ5aEbDrSODvMu4fQHQenXgCfwDdxVnzQ5aUCogXRcjvAG/KFJCS6SpzTJpR28tClJVUZXc5AQUrwl36OlrXc29AnCDutunknry6EyzT6cDvYjMAP4XGAj8n6ouTJlfBdwBfBCfdK+q/l983lzgpvj076nqw9m21WWB3tVe8HUDWmrwTA2s7u9e3nUyNRjnyq3xV1Y6J6jUNn43qrjLeZcPhdKXN5NQKGWwe4ZA787PJhh09jtdddn9Lt3tuUNtMo1+8h4Hd7sdPZl59zHdndIFOIG4NX3vYS8uhvnznU3HYtBwaog9b4ecmbn0B8wPOicC753I0Uhr7d8nijVAs9RTvCZMSV0EqkIE48erthYCdwWhrop/FUXZPayu+050BeTuI7S9OnC5J4SO6FSgF5GBwDrgNKABWAnMUdW1nmWqgMmqelXKugcBq4DJgAKrgUmquinT9ro80OfajJMaZL1Pq0p381Uun9sZ3jb9XK4U3MDsbX5Jd4Jyy+1eQaS76kn9rNTg6/6eun5qsPc2WWX6HtOVy/2ZeuKNRJIDeCbeZijvcUptNstUjtTvq50TgnuDmbt4uvmpdYAhQ9oZRhqKtK39Q3JzkHsScK8c3JvQoPUO42xNQnsHwoA9WQrRDu/2OsP9nKZA6/72syR7HUme19lAfwIQUdUz4u9vBFDV73uWqSJ9oJ8DhFT1G/H3PwGiqvqrTNvr8kAPyXf+dIa3OaOrA32+UgNptgCXD/dKxXvyy9ZM5S1HJqlXTLmWK93Vjst7Imhv+6mf473SSD0BFKiDOl3HcoeGkXqbg7xDQ90moFAEglGG/CrqnEiyXQnkEqhTcxA1BZJvVHObnGLxn9GIc1XS0T6HXE5QORi4OZD+0ZrxFB5tOtE70i9RwL6Vngj0FwAzVPVr8feXAMd7g3o80H8faMSp/V+rqu+LyLeBIlX9Xny5m4HtqvqjlG1cBlwGcPjhh0+q72wAzlW6qlV/4DaftBcAs/H2LbQ3fnHIECgqym2MY85jIT3S9R+I5HaSae9z3asPSB/4Cyhd/SOR2yfbn2k8uLqjZFKXdS+sEhdDUWdceSL3UKwycddwYK5z53Ja7omkKuQsn+6+g0zc4aiFqPHno6UYVswnXBlhUV0kqf2+Up0rIPeE21QRoaQu4qT7EGk9iZXEkrO2ZtqPiDrfjTsv9aSYp3ybcboj0B8MbFHVHSLyDeBCVT0l10Dv1S01eihsDbyzAcUUhrepqDvusuqiezCCwdbfvSN9MrWsuWMJ3KK463tTVbif4S2uVDvBt/LdaGI77k/v2HNnRmXbkUTekUZxlZUZ7mb2XmFA4UfyZKhNF9+liTTcQCKdxL/+BS0LYxk/ruTcCBVNkeTv202z4d5A5+6DOzrJO+oqfuJMJObLs6bfK5tuUpYfCHysqsW9tunG1duaW0zHubXu7jjh9uBINe9FaGqAT7esG7jTcZ8/kG59bzoF6qpoWhJp90LLOwTVvZ8h3TrFxdB8rXQutYS7TqySIRtCTubUeDu+NAfQNVUA6PI0O5eHrCHCDcTZOtWrQpkf5OPNzOo5cXR05E22QI+qZn0Bg4B3gNHAPsA/gHEpy4zw/H4esCL++0HAu8DQ+Otd4KBs25s0aZJ2O+dft2te4XDyz0Cg9Xc/voqLe74MuR6Tjr4qK7v/b9QjHHa+5i7fzvJw0vvKSufPN90hz/SVZPrqwsvDzmeFws68+QHnZ4T2XwuKNRBQHfKNytbPC6tWhsNa+VBloqzhsBaE++dSWZnfn0nSdxUKO/sYCrf+nrpOBKWqMunPNB/AKtUMcTzTjKSF4Cyctve3gf+JT7sVmBX//fvAa/GTwHJgrGfdS4H18de89rbVI4E+3yOYT0BRjX/NmnzkvMHG/d1d1js9W1AKBJLL7q6f7r+xN7268kTXXfseCHTpn2Rv5gbRysr2z3nhcOshSRe4vEE0HFatfKgyEdCLw4HW4B4KayBcqeHl4R45z7pld/fH/desrGyt27gnQvenauZl077ckx6Zv69sOh3ou/PVI4FeNTkoFyJopAb1TNt053mrJt7AnTotdZ53fmdOWIUKkO5fdGeXcV+97QqhUNXEfibb15Y6L7w8nLiaCNwV0MBdgZ6+iGpXPuVr71+toyzQ5yPdN546LVuQKsQ1o/cE4PJeN2eqGnl/epdzy5Wt3KlNTKn7nmvATZVuXiFq9Pme1NpbPtcTnTEFkO3PraMs0OfDDbLpmlm8te9MwaI7ypeL1MDq/sxWS3XLn63JKPWKI/XaO105sv0Fp36uOy31GKSeaF3tBfB028pULu93ke1kZLV6UwCBQHJTT3Fx5/60LNAXQqYj0Fv/+bNdS6ZrHmpvWa/U/W3vujzdFUq6z/KePLItl1q21NeQIenL7S1H6slbte13lnoyM6YXs0DflXpjkM9FPuVOFzALJZ8TRiYdvcJqbz+861ugN71ctkDff7JXdpVccq30Rvkk8SpAwq8u570Ttiv+pvvCd2D6NUtTbPwvNduoMf1MtkA/oLsLY0yXcGvbbqZNY0yCBXrjL9a8YkwbFuiNMcbnLNAbY4zPWaA3xhifs0BvjDE+Z4HeGGN8zgK9Mcb4nAV6Y4zxOQv0xhjjcxbojTHG5yzQG2OMz1mgN8YYn7NAb4wxPmeB3hhjfM4CvTHG+JwFemOM8TkL9MYY43MW6I0xxucs0BtjjM9ZoDfGGJ/LKdCLyAwReVNE1ovIgizLnS8iKiKT4++DIrJdROrir/sLVXBjjDG5GdTeAiIyELgPOA1oAFaKyFJVXZuy3IHANcDLKR/xtqpWFKa4xhhj8pVLjX4KsF5V31HVncBiYHaa5W4DfgC0FLB8xhhjOimXQD8SeN/zviE+LUFEJgKHqepv06w/WkTWiEitiJzY8aIaY4zpiHabbtojIgOAO4GqNLM3AIer6kYRmQQsEZFxqvpJymdcBlwGcPjhh3e2SMYYYzxyqdF/ABzmeT8qPs11IDAeiIpIDJgKLBWRyaq6Q1U3AqjqauBt4OjUDajqT1V1sqpOHj58eMf2xBhjTFq5BPqVwBgRGS0i+wAXAUvdmararKrDVDWoqkFgBTBLVVeJyPB4Zy4icgQwBnin4HthjDEmo3abblR1t4hcBTwHDAQeVNXXRORWYJWqLs2y+knArSKyC9gLXK6qHxei4MYYY3IjqtrTZUgyefJkXbVqVU8Xwxhj+hQRWa2qk9PNsztjjTHG5yzQG2OMz1mgN8YYn7NAb4wxPmeB3hhjfM4CvTHG+JwFemOM8TkL9MYY43MW6I0xxucs0BtjjM9ZoDfGGJ+zQG+MMT5ngd4YY3zOAr0xxvicBXpjjPE5C/TGGONzFuiNMcbnLNAbY4zPWaA3xhifs0BvjDE+Z4HeGGN8zgK9Mcb4nAV6Y4zxOQv0xhjjcxbojTHG5yzQG2OMz1mgN8YYn7NAb4wxPpdToBeRGSLypoisF5EFWZY7X0RURCZ7pt0YX+9NETmjEIU2xhiTu0HtLSAiA4H7gNOABmCliCxV1bUpyx0IXAO87Jl2LHARMA74DLBMRI5W1T2F2wVjjDHZ5FKjnwKsV9V3VHUnsBiYnWa524AfAC2eabOBxaq6Q1XfBdbHP88YY0w3abdGD4wE3ve8bwCO9y4gIhOBw1T1tyJyfcq6K1LWHZm6ARG5DLgs/naLiLyZQ7kyGQb8pxPr90W2z/7X3/YXbJ/zFcg0I5dAn5WIDADuBKo6+hmq+lPgp50tS7w8q1R1cvtL+ofts//1t/0F2+dCyiXQfwAc5nk/Kj7NdSAwHoiKCMCngaUiMiuHdY0xxnSxXNroVwJjRGS0iOyD07m61J2pqs2qOkxVg6oaxGmqmaWqq+LLXSQiQ0RkNDAG+FvB98IYY0xG7dboVXW3iFwFPAcMBB5U1ddE5FZglaouzbLuayLya2AtsBu4shtG3BSkCaiPsX32v/62v2D7XDCiql3xucYYY3oJuzPWGGN8zgK9Mcb4nG8Cfa5pGvoaETlMRJaLyFoReU1ErolPP0hEnheRt+I/h8ani4jcHf8eXonf49AnichAEVkjIs/G348WkZfj+/ZYfHAA8c7+x+LTXxaRYI8WvINEpEREnhCRN0TkdRE5we/HWUSujf9d/1NEfiUiRX47ziLyoIj8W0T+6ZmW93EVkbnx5d8Skbn5lMEXgd6TpuFM4FhgTjz9gh/sBv5bVY8FpgJXxvdtAfBHVR0D/DH+HpzvYEz8dRnw/7q/yAVzDfC65/0PgLtU9ShgE/DV+PSvApvi0++KL9cX/S/we1UdC5Tj7Ltvj7OIjASuBiar6nicwR4X4b/jXAPMSJmW13EVkYOAMM7NqlOAsHtyyImq9vkXcALwnOf9jcCNPV2uLtrXp3HyDr0JjIhPGwG8Gf/9J8Acz/KJ5frSC+eeiz8CpwDPAoJzx+Cg1GOOMyLshPjvg+LLSU/vQ577Wwy8m1puPx9nWu+6Pyh+3J4FzvDjcQaCwD87elyBOcBPPNOTlmvv5YsaPenTNLRJtdDXxS9VJ+AkjjtUVTfEZ/0LODT+u1++i0XADcDe+PuDgSZV3R1/792vxD7H5zfHl+9LRgONwEPx5qr/E5H98fFxVtUPgB8B7wEbcI7bavx9nF35HtdOHW+/BHrfE5EDgCeB+ar6iXeeOqd434yTFZGZwL9VdXVPl6UbDQImAv9PVScAW2m9nAd8eZyH4iQ+HI2T3XZ/2jZx+F53HFe/BHpfp1oQkcE4Qf5RVX0qPvkjERkRnz8C+Hd8uh++i2nALBGJ4WRLPQWn/bpERNyb/Lz7ldjn+PxiYGN3FrgAGoAGVXXTfD+BE/j9fJxPBd5V1UZV3QU8hXPs/XycXfke104db78E+qxpGvoyERHgZ8DrqnqnZ9ZSwO15n4vTdu9O/0q8934q0Oy5ROwTVPVGVR2lTkqNi4A/qerFwHLggvhiqfvsfhcXxJfvUzVfVf0X8L6IHBOfNB3njnLfHmecJpupIrJf/O/c3WffHmePfI/rc8DpIjI0fiV0enxabnq6k6KAnR1nAeuAt4H/6enyFHC/PodzWfcKUBd/nYXTNvlH4C1gGXBQfHnBGYH0NvAqzoiGHt+PTux/CHg2/vsROLmS1gOPA0Pi04vi79fH5x/R0+Xu4L5WAKvix3oJMNTvxxmoBt4A/gk8Agzx23EGfoXTB7EL58rtqx05rsCl8X1fD8zLpwyWAsEYY3zOL003xhhjMrBAb4wxPmeB3hhjfM4CvTHG+JwFemOM8TkL9MYY43MW6I0xxuf+P4ZC7LsmCtdSAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "from matplotlib import pyplot as plt\n", "fig, ax = plt.subplots()\n", "\n", "plt.ylim([0.4,0.7])\n", "ax.plot(range(len(plaquette_local_metropolis)), plaquette_local_metropolis, marker='+', ls='', c='blue', label=\"Local metropolis\")\n", "ax.plot(range(len(plaquette_heatbath)), plaquette_heatbath, marker='+', ls='', c='red', label=\"SU(2) Heatbath\")\n", "ax.plot(range(len(plaquette_langevin)), plaquette_langevin, marker='+', ls='', c='green', label=\"Langevin (BF scheme, eps=0.02)\")\n", "\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "3a4f16e1", "metadata": {}, "source": [ "**Homework**: Measure the average 1x1 and 2x1 Wilson loops for $\\beta=5.5$. Write your own implementation based on the **g.cshift** function and compare to **g.qcd.gauge.rectangle**. See [this test](https://github.com/lehner/gpt/blob/master/tests/qcd/gauge.py) for an explanation of the latter function." ] }, { "cell_type": "markdown", "id": "86931375", "metadata": {}, "source": [ "## Haar integration for SU(2)\n", "\n", "Before explaining in detail how the heat bath algorithm works, we first discuss the Haar integration for SU(2) in some detail.\n", "\n", "### Haar measure for explicit parametrization\n", "If we parametrize a general 2x2 matrix\n", "$$\n", "\\begin{align}\n", "U = \\begin{pmatrix} a & b \\\\ c & d \\end{pmatrix}\n", "\\end{align}\n", "$$\n", "its inverse (if it exists) can be written as\n", "$$\n", "\\begin{align}\n", "U^{-1} = \\frac1{\\det(U)} \\begin{pmatrix} d & -b \\\\ -c & a \\end{pmatrix}\n", "\\end{align}\n", "$$\n", "such that for $U \\in \\textrm{SU}(2)$ we find\n", "$$\n", "\\begin{align}\n", "U^{-1} = \\begin{pmatrix} d & -b \\\\ -c & a \\end{pmatrix} = U^\\dagger = \\begin{pmatrix} a^* & c^* \\\\ b^* & d^* \\end{pmatrix}\n", "\\end{align}\n", "$$\n", "and therefore\n", "$$\n", "\\begin{align}\n", "d &= a^* \\,, & c &= -b^* \\,.\n", "\\end{align}\n", "$$\n", "A general SU(2) matrix can thus be written as\n", "$$\n", "\\begin{align}\n", "U = \\begin{pmatrix} a & b \\\\ -b^* & a^* \\end{pmatrix}\n", "\\end{align}\n", "$$\n", "with\n", "$$\n", "\\begin{align}\n", "\\vert a \\vert^2 + \\vert b \\vert^2 = 1 \\,.\n", "\\end{align}\n", "$$\n", "can therefore also be expressed as\n", "$$\n", "\\begin{align}\n", "U = \\Sigma \\cdot x\n", "\\end{align}\n", "$$\n", "with $a = x_0 + i x_3$, $b = x_2 + i x_1$ and\n", "$$\n", "\\begin{align}\n", "\\Sigma_0 &= \\mathbb{1} \\,, \\\\\n", "\\Sigma_i &= i \\sigma_i\n", "\\end{align}\n", "$$\n", "and Pauli matrices\n", "$$\n", "\\begin{align}\n", "\\sigma_1 &= \\begin{pmatrix} 0 & 1 \\\\ 1 & 0 \\end{pmatrix} \\,, &\n", "\\sigma_2 &= \\begin{pmatrix} 0 & -i \\\\ i & 0 \\end{pmatrix} \\,, &\n", "\\sigma_3 &= \\begin{pmatrix} 1 & 0 \\\\ 0 & -1 \\end{pmatrix} \\,.\n", "\\end{align}\n", "$$\n", "such that\n", "$$\n", "\\begin{align}\n", "\\det(U) = \\vert x \\vert^2 = 1 \\,.\n", "\\end{align}\n", "$$\n", "This parametrization allows for a straightforward definition of a Haar integration measure.\n", "It is straightforward to show that for\n", "$$\n", "\\begin{align}\n", "V &= \\Sigma \\cdot y \\,, \\\\\n", "U V &= \\Sigma \\cdot z \\,,\n", "\\end{align}\n", "$$\n", "we have\n", "$$\n", "\\begin{align}\n", "\\det\\left( \\frac{\\partial z_\\mu}{\\partial x_\\nu} \\right) = \\det(V)^2 = 1\n", "\\end{align}\n", "$$\n", "and therefore\n", "$$\n", "\\begin{align}\n", "d[U] = dx_0 dx_1 dx_2 dx_3 \\delta(\\vert x \\vert - 1)\n", "\\end{align}\n", "$$\n", "is a Haar integration measure. For concreteness, we can transform this to spherical coordinates\n", "with\n", "$$\n", "\\begin{align}\n", " x_0 &= \\cos(\\varphi_0) \\,, \\\\\n", " x_1 &= \\sin(\\varphi_0)\\cos(\\varphi_1) \\,, \\\\\n", " x_2 &= \\sin(\\varphi_0)\\sin(\\varphi_1)\\cos(\\varphi_2) \\,, \\\\\n", " x_3 &= \\sin(\\varphi_0)\\sin(\\varphi_1)\\sin(\\varphi_2) \\,,\n", "\\end{align}\n", "$$\n", "and write\n", "$$\n", "\\begin{align}\n", " \\int d[U] f(U) = \\int_0^{\\pi} d\\varphi_0 \\int_0^{\\pi} d\\varphi_1 \\int_0^{2\\pi} d\\varphi_2 \\sin(\\varphi_0)^2 \\sin(\\varphi_1) f(U) \\,.\n", "\\end{align}\n", "$$\n", "Therefore we obtain the standard normalization\n", "$$\n", "\\begin{align}\n", "\\int d[U] = 1\n", "\\end{align}\n", "$$\n", "for\n", "$$\n", "\\begin{align}\n", " \\int d[U] f(U) = \\frac1{2\\pi^2}\\int_0^{\\pi} d\\varphi_0 \\int_0^{\\pi} d\\varphi_1 \\int_0^{2\\pi} d\\varphi_2 \\sin(\\varphi_0)^2 \\sin(\\varphi_1) f(U)\n", "\\end{align}\n", "$$\n", "or\n", "$$\n", "\\begin{align}\n", "d[U] = \\frac1{2\\pi^2} dx_0 dx_1 dx_2 dx_3 \\delta(\\vert x \\vert - 1) \\,.\n", "\\end{align}\n", "$$\n", "\n", "One can show that the Haar integration satisfies\n", "$$\n", "\\begin{align}\n", "\\int_{\\textrm{SU}(2)} d[U] U_{ab} U^\\dagger_{cd} = \\frac12 \\delta_{ad} \\delta_{bc} \\,.\n", "\\end{align}\n", "$$\n", "It is straightforward to show this using the explicit parametrization that we just obtained,\n", "however, we can also check this against our simple stochastic sampling implementation discussed\n", "in the previous subsection.\n", "\n", "### Stochastic integration using Haar measure" ] }, { "cell_type": "code", "execution_count": 3, "id": "a1d1fa04", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\\int dU U_00 U^\\dagger_00 = (0.4963044822891241, 0.002874731877928811)\n", "\\int dU U_00 U^\\dagger_01 = (0.007828394800862332, 0.0028740171419067277)\n", "\\int dU U_00 U^\\dagger_10 = (-0.001844956947146266, 0.0028647142161262475)\n", "\\int dU U_00 U^\\dagger_11 = (0.0026751189098447098, 0.0040610681617341)\n", "\\int dU U_01 U^\\dagger_00 = (0.0004471683358340047, 0.0029073253833936242)\n", "\\int dU U_01 U^\\dagger_01 = (0.0034909122072718794, 0.004106166013781911)\n", "\\int dU U_01 U^\\dagger_10 = (0.4978861456473543, 0.0028930697294849385)\n", "\\int dU U_01 U^\\dagger_11 = (0.0026898438724813297, 0.0028914647379603485)\n", "\\int dU U_10 U^\\dagger_00 = (0.000477677892407116, 0.002863568675674644)\n", "\\int dU U_10 U^\\dagger_01 = (0.49465169552991295, 0.0029002791419250867)\n", "\\int dU U_10 U^\\dagger_10 = (0.0025135099129404356, 0.00402159882946753)\n", "\\int dU U_10 U^\\dagger_11 = (0.008233227681888599, 0.0028798633019951943)\n", "\\int dU U_11 U^\\dagger_00 = (0.004582479741889893, 0.004133579983541516)\n", "\\int dU U_11 U^\\dagger_01 = (0.0020417697565920964, 0.0028761547175931634)\n", "\\int dU U_11 U^\\dagger_10 = (-0.004925965656714123, 0.0028903740046162135)\n", "\\int dU U_11 U^\\dagger_11 = (0.4987628291665141, 0.002880979866428887)\n" ] } ], "source": [ "import gpt as g\n", "import numpy as np\n", "import scipy.linalg\n", "\n", "pauli = [np.matrix([[0,1],[1,0]]), np.matrix([[0,-1j],[1j,0]]), np.matrix([[1,0],[0,-1]])]\n", "U = [np.matrix(scipy.linalg.expm(2j*sum([np.random.normal()*p for p in pauli]))) for i in range(40)]\n", "for u in U:\n", " assert np.linalg.norm(u * u.H - np.identity(2)) < 1e-14\n", " \n", "def haar_integrate(f, U, N, NskipToMakeIndependent):\n", " u0 = U[np.random.randint(0,len(U))]\n", " val = 0.0\n", " val2 = 0.0\n", " n = 0.0\n", " for i in range(N):\n", " v = f(u0).real\n", " val += v\n", " val2 += v**2\n", " n += 1\n", " for j in range(NskipToMakeIndependent):\n", " u0 = U[np.random.randint(0,len(U))] * u0\n", " # just use biased estimator for error\n", " return (val/n, (val2/n - (val/n)**2)**0.5 / n**0.5 )\n", "\n", "for a in range(2):\n", " for b in range(2):\n", " for c in range(2):\n", " for d in range(2):\n", " print(f\"\\int dU U_{a}{b} U^\\dagger_{c}{d} =\",\n", " haar_integrate(lambda u: u[a,b]*u.H[c,d], U, 10000, 2))" ] }, { "cell_type": "markdown", "id": "2395cc2f", "metadata": {}, "source": [ "### Haar measure using invariant length elements\n", "\n", "Before moving on to the heat bath, we observe that there is also another straightforward way to determine the Haar integration measure by studying invariant length elements\n", "$$\n", "\\begin{align}\n", "ds^2 = \\textrm{Tr}\\: \\left[ (dU U^{\\dagger}) (dU U^{\\dagger})^\\dagger \\right]\n", "\\end{align}\n", "$$\n", "that by construction is invariant under $U \\to V U$ and $U \\to U V$. For the angular parametrization\n", "of SU(2), this yields\n", "$$\n", "\\begin{align}\n", "ds^2 &= 2 d\\varphi_0^2 + 2 d\\varphi_1^2 \\sin(\\varphi_0)^2 + 2 d\\varphi_2^2 \\sin(\\varphi_0)^2 \\sin(\\varphi_1)^2 \n", "= \\begin{pmatrix} d\\varphi_0 & d\\varphi_1 & d\\varphi_2 \\end{pmatrix} G \\begin{pmatrix} d\\varphi_0 \\\\ d\\varphi_1 \\\\ d\\varphi_2 \\end{pmatrix}\n", "\\end{align}\n", "$$\n", "with\n", "$$\n", "\\begin{align}\n", "G = \\begin{pmatrix} 2 & & \\\\ & 2\\sin(\\varphi_0)^2 & \\\\ & &2 \\sin(\\varphi_0)^2 \\sin(\\varphi_1)^2 \\end{pmatrix} \\,.\n", "\\end{align}\n", "$$\n", "Therefore there is a different basis of 1-forms $d\\beta_\\mu$\n", "for which the length element is simply\n", "$$\n", "\\begin{align}\n", "ds^2 = d\\beta_0^2 + d\\beta_1^2 + d\\beta_2^2\n", "\\end{align}\n", "$$\n", "corresponding to the integration measure\n", "$$\n", "\\begin{align}\n", "d[U] = d\\beta_0 d\\beta_1 d\\beta_2\n", "\\end{align}\n", "$$\n", "and\n", "$$\n", "\\begin{align}\n", "\\begin{pmatrix} d\\beta_0 \\\\ d\\beta_1 \\\\ d\\beta_2 \\end{pmatrix} = G^{1/2} \\begin{pmatrix} d\\varphi_0 \\\\ d\\varphi_1 \\\\ d\\varphi_2 \\end{pmatrix}\n", "\\end{align}\n", "$$\n", "and thus\n", "$$\n", "\\begin{align}\n", "d[U] = d\\varphi_0 d\\varphi_1 d\\varphi_2 \\sqrt{\\det(G)} = \\sqrt{8} d\\varphi_0 d\\varphi_1 d\\varphi_2 \\sin(\\varphi_0)^2 \\sin(\\varphi_1)\n", "\\end{align}\n", "$$\n", "which agrees with our previous result up to the constant factor that only changes the normalization of the integral.\n", "Since we could re-define the invariant length-element $ds^2$ by a constant prefactor, this procedure gives us the properly normalized Haar measure only after normalizing the integral to $\\int d[U]=1$." ] }, { "cell_type": "markdown", "id": "dc39e9de", "metadata": {}, "source": [ "## Heat bath algorithm\n", "\n", "Let us now consider the heat bath algorithm to perform the integral over a single gauge link $U \\in$ SU(2). In order to extend this to integrals over SU(3) or more generally SU($N_c$), we can then consider the SU(2) subgroups of SU($N_c$) and re-write the integral over SU($N_c$) as a sufficient number of integrals over SU(2) subgroups. This is what the heat bath studied numerically in the last lecture does.\n", "\n", "With $U \\in$ SU(2), we need to integrate over\n", "$$\n", "\\begin{align}\n", "dP(U) = d[U] \\exp\\left(\\frac{\\beta}{2} \\textrm{Re}\\: \\textrm{Tr}\\: U A \\right)\n", "\\end{align}\n", "$$\n", "with staple matrix $A$, which itself is a sum over SU(2) matrices. As we have seen above, each SU(2) matrix can\n", "be written of the form $\\Sigma \\cdot x$ such that sums of SU(2) matrices always have to be proportional to a SU(2) matrix again. Therefore, we can make the ansatz\n", "$$\n", "\\begin{align}\n", "A = a V\n", "\\end{align}\n", "$$\n", "with\n", "$$\n", "\\begin{align}\n", "a = \\sqrt{\\det(A)} \\,.\n", "\\end{align}\n", "$$\n", "It can also be shown that $\\det(A) \\geq 0$. If $\\det(A)=0$, then $a=0$ and $A=0$ such that we\n", "can simply integrate over the Haar measure. If $a\\neq 0$, we proceed as follows.\n", "\n", "We define a matrix $X=U V$ in SU(2) and use the Haar invariance $d[U]=d[X]$ to find\n", "$$\n", "\\begin{align}\n", "dP(U) = d[X] \\exp\\left(\\frac{\\beta a}{2} \\textrm{Re}\\: \\textrm{Tr}\\: X \\right)\\,.\n", "\\end{align}\n", "$$\n", "The strategy is now to generate matrices $X$ following this distribution and to use $U=X V^\\dagger=X A^\\dagger / a$ to recover the desired gauge link $U$.\n", "We write\n", "$$\n", "\\begin{align}\n", "X = \\Sigma \\cdot x\n", "\\end{align}\n", "$$\n", "with\n", "$$\n", "\\begin{align}\n", "d[X] = \\frac1{2\\pi^2} dx_0 dx_1 dx_2 dx_3 \\delta(\\vert x \\vert - 1)\n", "\\end{align}\n", "$$\n", "and\n", "$$\n", "\\begin{align}\n", "\\textrm{Tr}\\: X = 2 x_0\n", "\\end{align}\n", "$$\n", "such that\n", "$$\n", "\\begin{align}\n", "dP(U) = \\frac1{2\\pi^2} dx_0 dx_1 dx_2 dx_3 \\delta(\\vert x \\vert - 1) \\exp\\left(\\beta a x_0 \\right)\\,.\n", "\\end{align}\n", "$$\n", "We note that only $x_0$ remains in the exponential and therefore it is useful to separate integral by writing\n", "$$\n", "\\begin{align}\n", "dP(U) &= \\frac1{2\\pi^2} dx_0 dx_1 dx_2 dx_3 \\delta(\\sqrt{x_0^2 + \\vert \\vec{x}\\vert^2} - 1) \\exp\\left(\\beta a x_0 \\right) \\\\\n", "&=\\frac1{2\\pi^2} dx_0 dx_1 dx_2 dx_3 \\frac{\\Theta(1-x_0^2)}{\\vert \\vec{x} \\vert} \\delta(\\vert \\vec{x} \\vert - \\sqrt{1 - x_0^2}) \\exp\\left(\\beta a x_0 \\right)\\,.\n", "\\end{align}\n", "$$\n", "with $\\vec{x}^T = (x_1,x_2,x_3)$. We now write $\\vec{x}$ in spherical coordinates $r, \\varphi,\\theta$ such that\n", "$$\n", "\\begin{align}\n", "dP(U) &=\\frac1{2\\pi^2} dx_0 d\\theta d\\varphi dr r^2 \\sin(\\theta) \\frac{\\Theta(1-x_0^2)}{r} \\delta(r - \\sqrt{1 - x_0^2}) \\exp\\left(\\beta a x_0 \\right) \\\\\n", "&=\\frac1{2\\pi^2} dx_0 d\\theta d\\varphi dr \\sin(\\theta) \\Theta(1-x_0^2) \\sqrt{1 - x_0^2}\\delta(r - \\sqrt{1 - x_0^2}) \\exp\\left(\\beta a x_0 \\right)\\,.\n", "\\end{align}\n", "$$\n", "\n", "We now need to sample $x_0$ from the distribution\n", "$$\n", "\\begin{align}\n", "dp(x_0) = dx_0 \\Theta(1-x_0^2) \\sqrt{1 - x_0^2}\\exp\\left(\\beta a x_0 \\right)\n", "\\end{align}\n", "$$\n", "which can be performed by first defining\n", "$$\n", "\\begin{align}\n", "x_0 = 1 - 2 \\lambda^2\n", "\\end{align}\n", "$$\n", "such that we need to sample $\\lambda \\in [0,1]$ from\n", "$$\n", "\\begin{align}\n", "dp(\\lambda) \\propto d\\lambda \\lambda^2 \\sqrt{1 - \\lambda^2}\\exp\\left(-2\\beta a \\lambda^2 \\right)\n", "\\end{align}\n", "$$\n", "this in turn can be generated by sampling four random numbers $r_1,r_2,r_3,r_4$ uniformly distribute in $]0,1]$\n", "and defining\n", "$$\n", "\\begin{align}\n", "\\lambda^2 = -\\frac1{2a\\beta}(\\ln(r_1) + \\cos(2\\pi r_2)^2\\ln(r_3))\n", "\\end{align}\n", "$$\n", "which we accept only if $r_4^2 \\leq 1 - \\lambda^2$.\n", "\n", "The value $r$ is then fixed by $r=\\sqrt{1 -x_0^2}$. The remaining random variables $z=\\cos(\\theta)$ and $\\varphi$\n", "can then be sampled from a uniform distribution $z\\in[-1,1]$ and $\\varphi\\in[0,2\\pi[$.\n", "\n", "The implementation of this algorithm in GPT can be studied [here](https://github.com/lehner/gpt/blob/master/lib/gpt/algorithms/markov/su2_heat_bath.py).\n", "\n", "## The Langevin algorithm\n", "\n", "The Langevin algorithm demonstrated above generates a Markov chain by integrating a stochastic differential equation. We will illustrate the algorithm for the simple case of a continuous random variable $x$ following the distribution\n", "$$\n", "\\begin{align}\n", "p(x) = e^{-S(x)} \\,,\n", "\\end{align}\n", "$$\n", "however, the generalization of the algorithm to integrals over\n", "gauge groups is straightforward albeit tedious.\n", "\n", "We define a discrete stochastic process by\n", "$$\n", "\\begin{align}\n", " x_{n+1} = x_{n} + \\Delta_n(x_n)\n", "\\end{align}\n", "$$\n", "with\n", "$$\n", "\\begin{align}\n", " \\Delta_n(x) &= -\\varepsilon S'(x) + \\sqrt{\\varepsilon} \\eta_n\n", "\\end{align}\n", "$$\n", "with $\\eta$ satisfying\n", "$$\n", "\\begin{align}\n", " \\langle \\eta^{2N+1}_n \\rangle &= 0 \\,, & \\langle \\eta_n \\eta_m \\rangle &= 2 \\delta_{nm}\n", "\\end{align}\n", "$$\n", "for all integer $N$.\n", "The details of distribution of $\\eta$ beyond the first two moments do not matter\n", "at this point.\n", "\n", "### Fokker-Planck equation\n", "For now, let us assume there is a stationary probability distibution of $x_n$ for sufficiently large $n$.\n", "Then the probability $P(n,x)$ to find the value $x$ at step $n$ follows from\n", "the **Fokker-Planck** equation\n", "$$\n", "\\begin{align}\n", " P(n+1,x_{n+1}) &= \\int dx_n P(n,x_n \\cap n+1,x_{n+1}) \\\\\n", "&= \\int dx_n P(n+1,x_{n+1} \\vert n,x_n) P(n,x_n)\n", "\\end{align}\n", "$$\n", "with\n", "$$\n", "\\begin{align}\n", "P(n+1,x_{n+1} \\vert n,x_n) &= \\left\\langle \\delta\\left( x_{n+1} - x_{n} - \\Delta_n(x_n) \\right) \\right\\rangle \\\\\n", "&=\\delta(x_{n+1} - x_{n}) \\notag\\\\\n", "&\\quad + \\varepsilon \\left[ S'(x_n)\\delta'(x_{n+1} - x_{n}) + \\delta''(x_{n+1} - x_{n})\\right] \\notag\\\\&\\quad + O(\\varepsilon^2) \\,.\n", "\\end{align}\n", "$$\n", "We can define the operation of a transition matrix $T$ on a function $f$ as\n", "$$\n", "\\begin{align}\n", " (T \\circ f)(x) &\\equiv \\int dy P(n+1,x \\vert n,y) f(y) \\\\\n", " &=\\biggl\\{ 1 + \\varepsilon \\partial_{x} \\left[ S'(x) + \\partial_x\\right] + O(\\varepsilon^2) \\biggr\\} f(x) \\,.\n", "\\end{align}\n", "$$\n", "Therefore\n", "$$\n", "\\begin{align}\n", " P(n+1,x_{n+1}) &= (T \\circ P(n))(x_{n+1})\n", "\\end{align}\n", "$$\n", "and more generally\n", "$$\n", "\\begin{align}\n", " P(n+N,x_{n+N}) &= (T^N \\circ P(n))(x_{n+N})\n", "\\end{align}\n", "$$\n", "with $T^N \\circ f = T^{N-1} \\circ (T \\circ f)$ and $T^1 \\circ f = f$.\n", "\n", "### Stationary probability distribution\n", "\n", "We now need to study the right-eigensystem of $T$. A stationary probability distribution $\\overline{P}$ requires\n", "the presence of an eigenvalue 1 and corresponds to its eigenvector. All other eigenvalues need to be smaller\n", "than one such that the stationary distribution is approached in the limit $n\\to\\infty$.\n", "\n", "The stationary distribution exists if $T \\overline{P} = \\overline{P}$ or equivalently\n", "$$\n", "\\begin{align}\n", "0 = \\partial_{x} \\left[ \\overline{P} S' + \\overline{P}'\\right] + O(\\varepsilon)\n", "\\end{align}\n", "$$\n", "or\n", "$$\n", "\\begin{align}\n", "\\overline{P}'(x) = -\\overline{P}(x) S'(x) + c_1 + O(\\varepsilon)\n", "\\end{align}\n", "$$\n", "has a solution for a constant $c_1$. This is solved by\n", "$$\n", "\\begin{align}\n", " \\overline{P}(x) &= c_0 e^{-S(x)} + c_1 e^{-S(x)} \\int_1^x dy \\: e^{S(y)}\n", "\\end{align}\n", "$$\n", "for all $c_0$. If we now insist that $\\overline{P}$ and its derivative vanish for sufficiently large $x$, only $c_1=0$ is allowed.\n", "On compact spaces one similarly imposes the compactification on $\\overline{P}$ and its derivatives to determine additional\n", "integration constants.\n", "The overall constant is determined by probability conservation.\n", "\n", "Therefore we have shown that if a stationary probability distribution exists, it follows the desired\n", "$$\n", "\\begin{align}\n", " \\overline{P}(x) &= e^{-S(x)}\\,.\n", "\\end{align}\n", "$$\n", "We next show that the system indeed converges to a stationary probability distribution up to the $O(\\varepsilon)$ corrections which we have previously neglected. One can construct higher-order stochastic differential equations, which eliminate additional powers of $\\varepsilon$.\n", "\n", "### Convergence\n", "\n", "A similarity transform of the operator\n", "$$\n", "\\begin{align}\n", "\\Delta T &= \\partial_{x} \\left[ S'(x) + \\partial_x\\right]\n", "\\end{align}\n", "$$\n", "by\n", "$$\n", "\\begin{align}\n", " \\tilde{\\Delta} T &= e^{S(x)/2} \\Delta T e^{-S(x)/2} \\\\\n", " &=e^{S(x)/2} \\partial_{x} e^{-S(x)/2} \\left[ \\frac12 S'(x) + \\partial_x\\right] \\\\\n", " &=-\\left[ \\frac12 S'(x) -\\partial_{x} \\right] \\left[ \\frac12 S'(x) + \\partial_x\\right] \\\\ &= - Q^\\dagger Q\n", "\\end{align}\n", "$$\n", "with\n", "$$\n", "\\begin{align}\n", " Q &= \\frac12 S'(x) + \\partial_x\n", "\\end{align}\n", "$$\n", "shows that the self-adjoint operator $\\tilde{\\Delta} T$ is real and semi-negative definite.\n", "Therefore, the stochastic process is driven towards the stationary solution $\\overline{P}$.\n", "\n", "The implementation of this algorithm in GPT can be studied [here](https://github.com/lehner/gpt/blob/master/lib/gpt/algorithms/markov/langevin.py).\n" ] }, { "cell_type": "markdown", "id": "404fa097", "metadata": {}, "source": [ "***Homework***: Calculate $\\langle x^2 \\rangle$ for $x$ following a normal distribution with 0 mean and variance 1 by implementing a Langevin process." ] }, { "cell_type": "code", "execution_count": null, "id": "8c11723f", "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 }