Improved Hamiltonian for Minkowski Yang-Mills Theory

TL;DR

This paper introduces an improved Hamiltonian for Minkowski Yang-Mills theory that reduces lattice artifacts, enabling more accurate investigation of infrared field dynamics, Chern-Simons number response, and Lyapunov exponents.

Contribution

The paper develops a higher-order corrected Hamiltonian for classical Minkowski Yang-Mills theory to minimize lattice spacing effects in simulations.

Findings

The improved Hamiltonian reduces lattice artifacts to within 10% for key quantities.

The maximal Lyapunov exponent's small $a$ limit differs by about 5%, indicating lattice effects influence its value.

Chern-Simons number response shows potential dependence on lattice regulation, but data is inconclusive.

Abstract

I develop an improved Hamiltonian for classical, Minkowski Yang-Mills theory, which evolves infrared fields with corrections from lattice spacing $a$ beginning at $O(a^4)$. I use it to investigate the response of Chern-Simons number to a chemical potential, and to compute the maximal Lyapunov exponent. Both quantities have small $a$ limits, in both cases within $10\% $ of the limit found using the unimproved (Kogut Susskind) Hamiltonian. For the maximal Lyapunov exponent the limits differ by about $5 \% $, significant at about $5 \sigma$, indicating that while a small $a$ limit exists, its value is corrupted by lattice artefacts. For the response of Chern-Simons number the statistics are not good enough to resolve $ 5 \% $ differences, but it seems possible in analogy with the Lyapunov exponent that the final answer depends on the lattice regulation.

Tables

Table 1. Table 1: Stepsize dependence of thermodynamic properties in the thermalization algorithm for the improved Hamiltonian lattice, for a 16 3 superscript 16 3 16^{3} lattice at lattice inverse temperature β L = 5 subscript 𝛽 𝐿 5 \beta_{L}=5 .

quantity	at $\Delta t=0.05$	at $\Delta t=0.20$	at $\Delta t=0.30$
$H_{U\mathrm{improved}}{\beta }_L/\left(3N^3\right)$	$1.031\pm .0016$	$1.1081\pm .0015$	$1.2278\pm .0016$
$1\times 1$ Wilson loop	$\mathrm{\hspace{0.33em}.84865}\pm .00024$	$\mathrm{\hspace{0.33em}.83774}\pm .00021$	$\mathrm{\hspace{0.33em}.82087}\pm .00025$
$1\times 2$ Wilson loop	$.72376\pm .00035$	$.70763\pm .00034$	$.68274\pm .00038$
$1\times 3$ Wilson loop	$.6201\pm .0005$	$.6005\pm .0005$	$.5703\pm .0006$
$1\times 4$ Wilson loop	$.5316\pm .0007$	$.5104\pm .0007$	$.4768\pm .0007$
$1\times 5$ Wilson loop	$.4560\pm .0008$	$.4338\pm .0008$	$.3988\pm .0008$
$2\times 2$ Wilson loop	$.5412\pm .0010$	$.5222\pm .0010$	$.4909\pm .0010$
$2\times 3$ Wilson loop	$.4121\pm .0009$	$.3925\pm .0008$	$.3592\pm .0009$
$2\times 4$ Wilson loop	$.3149\pm .0010$	$.2967\pm .0009$	$.2644\pm .0010$
$2\times 5$ Wilson loop	$.2410\pm .0010$	$.2246\pm .0010$	$.1949\pm .0010$
$3\times 3$ Wilson loop	$.2843\pm .0015$	$.2660\pm .0015$	$.2344\pm .0015$
$3\times 4$ Wilson loop	$.1975\pm .0011$	$.1824\pm .0010$	$.1548\pm .0010$
$3\times 5$ Wilson loop	$.1377\pm .0011$	$.1255\pm .0011$	$.1018\pm .0010$
$4\times 4$ Wilson loop	$.1251\pm .0015$	$.1143\pm .0015$	$.0924\pm .0015$
$4\times 5$ Wilson loop	$.0848\pm .0011$	$.0722\pm .0011$	$.0546\pm .0010$
$5\times 5$ Wilson loop	$.0467\pm .0015$	$.0424\pm .0015$	$.0286\pm .0015$

Table 2. Table 2: Rate of motion of N C ​ S subscript 𝑁 𝐶 𝑆 N_{CS} under a chemical potential at several lattice spacings.

${\beta }_L$	$N$	number of runs	${\kappa }_{\mu }$
3	16	3	$.181\pm .012$
4	16	6	$.310\pm .017$
5	16	16	$.389\pm .021$
6	16	27	$.413\pm .024$
8	24	16	$.487\pm .033$
10	30	18	$.537\pm .036$

Table 3. Table 3: Lyapunov exponent for improved and Kogut-Susskind Hamiltonian at several lattice spacings. The improved data is new, the KS data is taken from [ 7 ] . Each Hamiltonian system shows a fine lattice limit, but the value of the limit differs by about 6 % percent 6 6\% between the Hamiltonian systems.

${\beta }_L$	$N$	improved ${\lambda }_{max}$	KS ${\lambda }_{max}$	improved ${\beta }_L{\lambda }_{max}$	KS ${\beta }_L{\lambda }_{max}$
2	12		$.659\pm .002$		$1.318\pm .004$
2.5	12	$.4932\pm .005$	$.586\pm .002$	$1.233\pm .013$	$1.465\pm .005$
3	12	$.399\pm .0034$	$.497\pm .004$	$1.197\pm .010$	$1.491\pm .012$
4	12	$.2863\pm .0021$	$.350\pm .003$	$1.145\pm .008$	$1.400\pm .012$
5	20	$.2280\pm .0018$	$.260\pm .002$	$1.140\pm .009$	$1.300\pm .010$
6	20	$.1908\pm .0010$		$1.145\pm .006$
7	20	$.1630\pm .0027$		$1.141\pm .019$
7.5	20		$.163\pm .001$		$1.222\pm .008$
8.75	20		$.137\pm .001$		$1.199\pm .009$
10	20		$.120\pm .001$		$1.20\pm .01$

Equations

$$L=\int d^{3}x\left(\frac{{\dot{\varphi }}^2}{2}-\frac{\overrightarrow{(\nabla }\varphi {)}^2}{2}-\frac{m^2}{2}{\varphi }^2-\frac{\lambda }{4}{\varphi }^4\right),$$

$$\frac{d\varphi }{dt}$$

$$\frac{d\pi }{dt}$$

$$\int d^{3}x(\vec{\nabla }\varphi {)}^2=a^3\sum _x\sum _{\hat{i}}\frac{(\varphi (x+a\hat{i})-\varphi (x){)}^2}{a^2}$$

$${\nabla }^2\varphi (x)\to \sum _{\hat{i}}\frac{\varphi (x-a\hat{i})+\varphi (x+a\hat{i})-2\varphi (x)}{a^2}$$

$$\sum _{\hat{i}}\frac{\varphi (x-a\hat{i})+\varphi (x+a\hat{i})-2\varphi (x)}{a^2}=\sum _i{\nabla }_i^2\varphi +\frac{a^2}{12}\sum _i{\nabla }_i^4\varphi +\frac{a^4}{360}\sum _i{\nabla }_i^6\varphi +\dots$$

$${\nabla }_{\mathrm{Latt}}^2$$

$$=$$

$$\frac{-1}{2}\left(\varphi (x)-\varphi (x-a\hat{i})\right)\left(\frac{-1}{12}\varphi (x+a\hat{i})+\frac{5}{4}\varphi (x)-\frac{5}{4}\varphi (x-a\hat{i})+\frac{1}{12}\varphi (x-2a\hat{i})\right),$$

$${\varphi }_x(x+a\hat{i})=\mathcal{P}\exp \left({\int }_{x+a\hat{i}}^x\frac{i{\tau }^a}{2}{A}_a^j(y)d{y}_j\right)\varphi (x+a\hat{i}),$$

$${\varphi }_{t+\Delta t}(x,t)=\mathcal{P}\exp \left({\int }_t^{t+\Delta t}i{\tau }^a{A}_0^a\right)\varphi (x,t)$$

$$U_j(x)U_i(x+j)U_j^{†}(x+i)U_i^{†}(x)\equiv U_{\square }(x)$$

$$\frac{1}{2}\mathrm{Tr}-i{\tau }^a{U}_{\square }(x),$$

$$-L_U=\sum _{\square }1-\frac{1}{2}\mathrm{Tr}{U}_{\square }$$

$$\frac{d{U}_i(x)}{dt}=i{\tau }^a{\dot{U}}_i^a(x)U_i(x)$$

$${\dot{U}}_i^a(x)=A_0(x)-U_i(x)A_0(x+i)U_i^{†}(x)+E_i^a(x).$$

$$L_E=\sum _{x,i}\frac{E_i^a(x)E_i^a(x)}{2}.$$

$$\frac{\partial L}{\partial A_0^a(x)}=0=\sum _j\left(E_j^a(x)-E_j^b(x-i)\frac{1}{2}\mathrm{Tr}-i{\tau }^a{U}_j^{†}(x-j)i{\tau }^b{U}_j(x-j)\right),$$

$$\frac{d{E}_i^a(x)}{dt}=-D{F}_i^a(x),D{F}_i^a(x)\equiv \sum _{4{\square }_i(x)}\frac{1}{2}\mathrm{Tr}-i{\tau }^a{U}_{\square },$$

$$E_i^a(x,t+\Delta t/2)$$

$$U_i(x,t+\Delta t)$$

$$I=\left(\sum _{{\square }_s}-\frac{1}{(\Delta t{)}^2}\sum _{{\square }_t}\right)\left[1-\frac{1}{2}\mathrm{Tr}{U}_{\square }\right],$$

$$U_i(x,t+\Delta t)$$

$$\frac{1}{2}\mathrm{Tr}-i{\tau }^a{E}_i(x,t+\Delta t/2)$$

$$\mathrm{Energy}(t)=\sum _{\square }1-\frac{1}{2}\mathrm{Tr}{U}_{\square }(t)+\sum _{x,i}\frac{(E_i^a(x,t+\Delta /2){)}^2+(E_i^a(x,t-\Delta /2){)}^2}{4}$$

$$\mathrm{Energy}(t)$$

$$\frac{d{N}_{CS}}{dt}=\int d^{3}x\frac{E_i^a(x)B_i^a(x)}{8{\pi }^2}.$$

$$\frac{d{N}_{CS}}{dt}=\sum _{x,i}\frac{E_i^a(x)B_i^a(x)}{2{\pi }^2},B_i^a(x)\equiv \frac{1}{8}\sum _{8\square }\frac{1}{2}\mathrm{Tr}-i{\tau }^a{U}_{\square }$$

$$2{\pi }^2\Delta N_{CS}=\Delta t\sum _{x,i}\frac{E_i^a(x,t+\Delta t/2)+E_i^a(x,t-\Delta t/2)}{2}{B}_i^a(x,t),$$

$$\frac{d{U}_i(x)}{dt}$$

$$\frac{d{E}_i^a(x)}{dt}$$

$$\sum _i\left(B_i^a(x)-B_i^a(x-i)\right)\ne 0$$

$$2{\pi }^2\frac{d{N}_{CS}}{dt}=E^c\cdot B+E^{\ast }\cdot B.$$

Full text

Preprint PUP-TH-1618

e-Print Archive: hep-lat/9605001

Improved Hamiltonian for Minkowski Yang-Mills Theory

Guy D. Moore111e-mail: guymoore@puhep1.princeton.edu

Princeton University

Joseph Henry Laboratories, PO Box 708

Princeton, NJ 08544, USA

Abstract

I develop an improved Hamiltonian for classical, Minkowski Yang-Mills theory, which evolves infrared fields with corrections from lattice spacing $a$ beginning at $O(a^{4})$. I use it to investigate the response of Chern-Simons number to a chemical potential, and to compute the maximal Lyapunov exponent. Both quantities have small $a$ limits, in both cases within $10\%$ of the limit found using the unimproved (Kogut Susskind) Hamiltonian. For the maximal Lyapunov exponent the limits differ by about $5\%$, significant at about $5\sigma$, indicating that while a small $a$ limit exists, its value is corrupted by lattice artefacts. For the response of Chern-Simons number the statistics are not good enough to resolve $5\%$ differences, but it seems possible in analogy with the Lyapunov exponent that the final answer depends on the lattice regulation.

1 Introduction

Several interesting and outstanding questions in high energy physics and cosmology, such as the production of quark-gluon plasma in relativistic heavy ion collisions and the conditions in the early universe under which baryogenesis may have occurred, involve finite temperature field theories where there are interactions between bosons, and some of the bosons have very small masses in comparison to the temperature, $m\ll T$. The infrared physics of such systems can be strongly coupled; in perturbation theory, the contribution of loops involving very soft ($O(g^{2}T)$, $g$ a generic coupling strength characterizing the theory) momenta can be order 1, even when $g$ is small enough that the vacuum theory is weakly coupled. This strongly coupled infrared physics is often very interesting; in the electroweak theory it may determine the strength of the electroweak phase transition, and it sets the rate of Chern-Simons number ($N_{CS}$) diffusion, which determines the rate of baryon number violation. The former problem is thermodynamical, and there are Euclidean tools which can be brought to bear on the nonperturbative physics involved[1, 2]. However, the latter question is dynamical, and cannot be addressed by Euclidean methods. There is also interest in a more general understanding of the infrared dynamics of Yang-Mills and Yang-Mills Higgs theory at finite temperature, for instance to understand its implications in heavy ion collisions. These questions cannot be addressed by Euclidean methods, which are only appropriate for answering thermodynamical questions. They are also very hard to address analytically in any reliable way, and the prospects for numerical evaluation of Minkowski path integrals are more or less hopeless.

Several years ago, Grigoriev and Rubakov made an interesting proposal for a new numerical approach to the understanding of the infrared problem in finite temperature bosonic systems[3]. The key observation is that, if the coupling constant $g$ is small, then the nonperturbative physics arises because occupation numbers are large in the infrared; but this is precisely the regime where a classical field approximation becomes reasonable. The physics is nonperturbative but in an essentially classical way. Therefore, it is reasonable to hope that, by analyzing the analogous classical field theory at finite temperature, one can understand at leading order the physics of the quantum system of interest.

Based on this idea, there has been a revival of interest in classical, thermal Yang-Mills theory[4, 5, 6, 7, 8, 9]. The work is based on lattice implementation of classical Yang-Mills theory, using the Hamiltonian formulism worked out by Kogut and Susskind years ago[10]. This work has yielded interesting results; it appears that we now know the $N_{CS}$ diffusion rate in the symmetric electroweak phase to within around $5\%$ accuracy[8], and there is numerical evidence for the interesting conjecture that the maximal Lyapunov exponent of Yang-Mills theory is twice the damping rate of a plasmon at rest[6].

However, there have been some frustrations in the application of this technology; for instance, in [9] I find that, to determine the rate of $N_{CS}$ diffusion in SU(3) gauge theory, the lattice spacing must be so fine that the calculation is computationally prohibitive. Also, there is a more fundamental challenge to the whole approach, raised by Bodecker et. al. [11], who consider the “hard thermal loops,” the one loop contributions of ultraviolet modes to the dispersion relations and interactions of the infrared modes. On the lattice, the strength of the hard thermal loops grows linearly with inverse lattice spacing, and as is often the case with linearly divergent quantities, the functional form of the hard thermal loops may depend on the form of the cutoff; they posess lattice artefacts which make their $\vec{q}/q_{0}$ dependence different than in the continuum, quantum theory. Although the results for $N_{CS}$ diffusion in SU(2) apparently have a fine lattice limit[8], which must be independent of the overall amplitude of the hard thermal loops, the value of the fine lattice limit may depend on the functional form of the hard thermal loops, and hence on the specifics of the lattice cutoff. I briefly explored this possibility in [9], where I found some evidence for a weak dependence on the form of the cutoff; however the case was far from clear, and it is an open and interesting question whether the $N_{CS}$ diffusion rate depends on the presence of hard thermal loops and on their functional form, and whether the convergence to a fine lattice limit represents the disappearance of nonrenormalizeable operators or a suppression of effects which die as a power of $m_{D}^{2}/m^{2}_{\rm nonperturbative}$.

One or both of these problems might be solved by finding an “improved” Hamiltonian for the classical, Minkowski Yang-Mills theory. By “improved” I mean a Hamiltonian such that the interaction of infrared modes would be the same as in the continuum, plus finite $a$ effects starting at $a^{4}$, rather than at $a^{2}$, as they do in the simplest lattice implementation. That is, for infrared fields the lattice action can be expanded in operator dimension, and an “improved” action is one in which, in addition to the desired dimension 4 $F_{\mu\nu}F^{\mu\nu}$ term, there are no dimension 6 operators, but only those starting at dimension 8. It is well known that in evolving scalar partial differential equations, and in Euclidean lattice gauge theory, such improvement greatly reduces the numerical demands by allowing accurate results on much coarser lattices[12]. Implementing such an improved Hamiltonian for the systems under consideration should either lighten the numerical demands, or verify that hard thermal loops are indeed important, either because the lattice fineness requirement arises from the need to make the Debye mass large, or because the results depend on the functional form of the hard thermal loops, and hence on the form of the regulation.

The basic idea of an improved lattice Hamiltonian is quite simple, and is easiest to present in the case of scalar field theory. Consider the evolution of the continuum system derived from the Lagrangian

[TABLE]

with $\phi$ a single component, real scalar field. The Hamilton’s equations of motion (defining $\pi$ as the conjugate momentum of $\phi$) are

[TABLE]

To implement these equations on a lattice, we discretize the fields $\phi$ and $\pi$ so they each take on a single value at each lattice point. The only difficulty lies in choosing a lattice definition of $\nabla^{2}$. If we make the replacement

[TABLE]

in the Lagrangian, then the appropriate substitution in the equation of motion is

[TABLE]

which is the standard lattice Laplacian.

Now consider the quality of this implementation when evolving some smooth field configuration. If the characteristic wavelength of the field configuation of interest is on order the lattice length, it is inevitable that the lattice implementation will poorly reproduce the continuum dynamics; so it is only interesting to test the quality of the lattice implementation for fields which are smooth on the lattice scale, in which case we can study the quality of the lattice Laplacian by expanding $\phi$ in a Taylor series about the point $x$; it is straightforward to show that

[TABLE]

where I write in terms of a positive 3 dimensional measure, as I will for the remainder of the paper. All but the first of the derivative terms are undesirable. The second term will produce lattice artefacts in the evolution of a mode of wave vector $k$ which will be on order $a^{2}k^{2}$; if $a\ll 1/k$, which is necessary for the lattice implementation to at all resemble the desired continuum system, then the quality of the implementation would be improved if we could write a new lattice Laplacian which did not produce an $O(a^{2})$ term.

This is quite easy; just include contributions from next nearest neighbors as well as nearest neighbors, and choose the coefficients so the coefficients of $\nabla^{2}$ and $\nabla_{i}^{4}$ are as desired. The correct choice is

[TABLE]

which is free of $O(a^{2})$ contamination. We can get this form if we replace the lattice Lagrangian term $-(\phi(x)-\phi(x-a\hat{(}i)))^{2}/2$ with

[TABLE]

which means we can still expect to have a conserved energy, which must however be defined in terms of the above gradient term.

This system is more complicated to evolve, perhaps doubling the number of computations in an update. What we have gained, though, is that the lattice size necessary to model the desired system with a given accuracy is greatly reduced; the number of lattice points needed is order the square root of the number previously required. This greatly reduces the numerical demands of the evolution.

The implementation of an improved lattice action is not quite as simple in nonabelian gauge theory, but the techniques involved have been worked out in the Euclidean case some time ago[13]. The purpose of this paper is to develop and implement the correct improved Minkowski Hamiltonian, and to use it to study the cutoff dependence of the system.

The paper is organized as follows. Sec. 2 will review the Kogut-Susskind implementation of Minkowski, classical lattice Yang-Mills theory, the lattice definition of $dN_{CS}/dt$, the addition of a chemical potential for $N_{CS}$, and the thermalization of the lattice system. Everything in the section has appeared in previous literature; I include it for didactic reasons, to review necessary tools and ideas for the remainder of the paper. In Sec. 3 I will construct a lattice Hamiltonian in continuum time which is free of dimension 6 lattice artefacts and find a discrete time algorithm for its evolution. The kinetic (electric field) part of the Hamiltonian is not diagonal in the obvious basis for electric fields, so the modifications to the definition of $dN_{CS}/dt$, chemical potential, and the thermalization algorithm are nontrivial. With the Hamiltonian, chemical potential, and thermalization fully implemented, I analyze the motion of $N_{CS}$ under a chemical potential in SU(2) gauge theory, presenting numerical results in Sec. 4. I also compute the maximal Lyapunov exponent for the system. Finally, Sec. 5 concludes.

2 Kogut-Susskind Hamiltonian and $N_{CS}$

In this section I will review the construction of a lattice, Minkowski model for classical Yang-Mills theory, developed by Kogut and Susskind [10]. I present a leapfrog algorithm for its update and a definition for $\Delta N_{CS}/\Delta t$ due to Ambjorn et. al. [4, 8], and a way to apply a chemical potential for $N_{CS}$ developed in [9]. Krasnitz has developed a thermalization algorithm for the system, based on a set of Langevin equations [7]; but I will present a much simpler thermalization algorithm developed in [9], because it is this algorithm I will extend to the improved model in Sec. 3.

2.1 Degrees of Freedom and Hamiltonian

To implement Yang-Mills theory on a lattice, Kogut and Susskind followed Wilson’s ideas [14], taking very literally the nature of the guage field as an affine connection for SU(2) objects in space. If there were a scalar field taking values in the fundamental representation of SU(2) at each point on a lattice, the gauge field should provide a rule for parallel transporting the field from lattice site to lattice site. In the continuum the rule is that the value of $\phi(x+a\hat{i})$, parallel transported to the point $x$, is

[TABLE]

with the path ordering placing points near $x$ to the left. (Here and throughout I will use nonrelativistic notation and a positive metric on the space indicies.) The lattice analog is to have a matrix $U\in$SU(2) which lives on the link between the lattice sites, so the parallel transported object is just $U_{i}(x)\phi(x+a\hat{i})$. The left index of $U$ lies at the basepoint of the link and the right index lies at its endpoint, so to perform a gauge transform in which $\phi(x)$ is rotated by $G\in$SU(2) to $G\phi$, the link $U_{i}(x)$ is rotated to $GU_{i}(x)$ and the link $U_{i}(x-a\hat{i})$ becomes $U_{i}(x-a\hat{i})G^{\dagger}$; so transporting $G\phi(x)$ back to $x-a\hat{i}$ gives $U_{i}(x-a\hat{i})G^{\dagger}G\phi(x)$, which is unchanged, as it should be, by the guage transform at $x$. To forward transport $\phi(x)$ to $x+a\hat{i}$ we multiply by $U_{i}^{\dagger}(x)$. Note that my convention is to label the link by its basepoint, so $U_{i}(x)$ runs from $x$ to $x+a\hat{i}$; also from now on I will drop $a$, measuring everything in lattice units, exept in those cases where it is useful to include it to count dimensions. I will also drop hats where the meaning is obvious, ie $x+i$ means $x+a\hat{i}$.

For the time being I will consider a spatial lattice but with continuous time. The connection in the time direction then remains a Lie algebra element $A_{0}^{a}$ living at each lattice point, with parallel transport of a field $\phi(x,t)$ to time $t+\Delta t$ given by

[TABLE]

where the path ordering places earlier times to the right. $A_{0}$ transforms under a gauge change as $A_{0}\rightarrow G(x,t)A_{0}(x,t)G^{\dagger}(x,t)+\partial_{0}G(x,t)$.

The field strength is a curvature tensor and it reflects the failure of an object, on transport around a closed curve, to return to its original value. The smallest nontrivial closed curve is the square plaquette, so $F_{ij}\neq 0$ in the continuum case corresponds to the product of the links around a plaquette,

[TABLE]

not being the identity. A reasonable definition of the magnetic field strength is the Lie algebra element

[TABLE]

which, although it transforms under gauge change as an adjoint object at the point $x$, is really associated with the plaquette lying in the positive $i$ and $j$ directions from $x$. The continuum Lagrangian depends on the $F_{ij}$ as $\int F_{ij}F_{ij}/4$, which is reproduced at leading order of a weak field expansion by

[TABLE]

where the sum is over all plaquettes in the lattice. (The coefficient does not agree with the continuum one, but I will account for this when I thermalize the system.)

Next consider the electric field part of the Lagrangian. It is convenient to write the time derivative of a link variable as

[TABLE]

and to expand the Lie algebra element $\dot{U}_{i}^{a}(x)$ as

[TABLE]

The field $E_{i}^{a}$ transforms under gauge transformations as an adjoint object and can be recognized in the weak field limit as the electric field; alternately we can notice that it is the small $\Delta t$ limit of a plaquette in the $i$ and time direction, traced against the Lie algebra and divided by $\Delta t$. The electric part of the Lagrangian should be

[TABLE]

Now we can derive Lagrange equations of motion. Note first that the time derivative of $A_{0}$ never appears, so $A_{0}$ enters the Lagrangian like a Lagrange multiplier. Varying with respect to $A_{0}$ gives

[TABLE]

where the second term is just the $E$ field on the link leading into $x$, parallel transported so its indicies reside at the point $x$. This equation is a constraint, called the Gauss constraint because it corresponds to the continuum condition that $D_{j}E_{j}=0$, which is Gauss’ Law. If we initialize the evolution by choosing values for the $U$ and $E$, then the Gauss constraint is a restriction on the allowed initial conditions.

The equations of motion derived by varying with respect to a link are

[TABLE]

where $4\Box_{i}(x)$ refers to the 4 plaquettes which contain the link leaving $x$ in the $i$ direction, with orientation chosen so that $U_{i}(x)$ and not $U^{\dagger}_{i}(x)$ appears, and beginning and ending at the point $x$. I use the notation $DF_{i}^{a}(x)$ because it is compact and remeniscient of the equivalent continuum quantity, $(D_{j}F_{ji})^{a}(x)$. Note that Eq. (19) automatically preserves the Gauss constraint if the initial conditions satisfy it, because the time derivative of Eq. (18) will, according to Eq. (19), consist of the sum of traces over plaquettes, and every plaquette which contributes to one $E$ contributes to a neighbor with opposite sign. This exact satisfaction is necessary because the evolution would not make any sense if it excited the modes constrained to be zero; in effect we keep track of more $E$ fields than are actually dynamical degrees of freedom, and we must make sure that only the dynamical degrees of freedom actually get excited.

Note that Eq. (19) does not actually uniquely specify the evolution of the system. $E_{i}^{a}(x)$ does not give us $\dot{U}_{i}^{a}(x)$, and the equations say nothing about how to evolve the value of $A_{0}$. This is because the evolution should not in fact be uniquely specified, yet; we have the freedom to make arbitrary time dependent gauge transformations. To write a numerical update algorithm we have to decide how to handle this freedom; the simplest and easiest choice is to pick the gauge such that $A_{0}^{a}$ is everywhere zero. The system can then be considered as a constrained Hamiltonian system with coordinates $U$ and conjugate momenta $E$. The Hamiltonian is just $L_{E}-L_{U}\equiv H_{E}+H_{U}$.

2.2 Numerical evolution of system

One way to evolve this system is to use a leapfrog algorithm in which the $U$ fields are defined at integer multiples of $\Delta t$, and the $E$ are taken to be Lie algebra elements living at half time steps; the rule is

[TABLE]

This choice still identically preserves the Gauss constraint [7], whether or not the parallel transport involved in the definition of the Gauss constraint is performed with the past or future $U$. The reason is that $E$ commutes with its own exponent, so the parallel transport is the same for the updated $U$ as for the $U$ before update. This algorithm forms the basis of the work in [7, 8, 9].

An alternate evolution algorithm developed in [4], which is closer to what I will need for the improved Hamiltonian system, is to begin with a discrete time lattice, so the $A_{0}$ fields are also connection matricies and the $E$ fields are plaquette matricies for plaquettes with one space and one time step. The action is

[TABLE]

where $\Box_{s}$ are plaquettes entirely in space and $\Box_{t}$ are plaquettes in one space and one time direction. The equations of motion in the temporal gauge follow immediately;

[TABLE]

The only differences between the two update algorithms are that matrix exponentiation is replaced by inverting $(1/2){\rm Tr}-i\tau^{a}$, and the definition of energy is changed from $E^{a}E^{a}/2$ to $(1-\sqrt{(\Delta t)^{2}E^{a}E^{a}})/(\Delta t)^{2}$. These changes are quadratic in $|E^{a}|\Delta t$, which is no larger than the natural level of error of the leapfrog algorithm, so the difference is unimportant.

For finite $\Delta t$ there is no identically preserved quantity which can be understood as a system energy; but a close reasonable approximation is the quantity

[TABLE]

for the first update algorithm and

[TABLE]

for the second. In each algorithm, the energy is preserved in the mean, with fluctuations suppressed by $(\Delta t)^{2}$ and further washed out because they generically add incoherently between different degrees of freedom. The central value is absolutely stable, because the update algorithm has two important properties: (1) it bases the evolution of the system on the minimum necessary amount of state information; and (2) it is exactly time reversal invariant (to evolve the system backwards, flip the sign of $E$ in the first version or take Hermitian conjugate in the second). It is therefore impossible that the energy should tend to rise exponentially, since under time reversal it would tend to shrink exponentially, but the algorithm is the same. further, by good fortune the roundoff errors which are inevitable in a numerical implementation add incoherently and do not tend to increase the energy. The stability of the system energy is necessary for long Hamiltonian evolutions and will also allow checks that the chemical potential method is behaving properly.

2.3 Magnetic field and $N_{CS}$

To track the evolution of $N_{CS}$ in the realtime evolution of the Yang-Mills fields described in the last section it is necessary to develop a definition of $dN_{CS}/dt$. (It is impossible to define $N_{CS}$ directly on a lattice, see [9]). In the continuum the relation is

[TABLE]

We have already found lattice analogs of $E_{i}$ and $F_{jk}\sim B_{i}$. $E_{i}$ is associated with a link, and $F_{jk}$ is associated with a plaquette in the $j,k$ directions (cyclic permutation implied), and the simplest definition of $dN_{CS}/dt$ which is invariant under the cubic point group is

[TABLE]

where $8\Box$ means a sum over the 8 plaquettes which are orthogonal to the link $x,i$ and touch one of its endpoints. The plaquettes should be traversed according to the righthand rule, starting and ending at the point along the link, see Figure 1; and parallel transport of the plaquettes which touch the endpoint of the link back to the basepoint is implied. In the discrete time, leapfrog algorithm evolution of the system the appropriate choice for $\Delta N_{CS}$ from one update is

[TABLE]

where $E_{i}^{a}(x)$ should be replaced by $(1/2\Delta t){\rm Tr}-i\tau^{a}E_{i}(x)$ for the algorithm used by Ambjorn et. al.

The coefficient is $2\pi^{2}$ rather than $8\pi^{2}$, as it is in the continuum, because the continuum equations are based on fields defined in terms of $\tau^{a}/2$, and throughout I use $\tau^{a}$ in the discrete system. Similarly the $g^{2}$, which can appear in the continuum equations depending on the normalization of $A_{\mu}$, does not appear here, but will appear in the definition of lattice temperature instead.

The definition of $\Delta N_{CS}$ above corresponds to that used in [8, 9]. The definition used in [4] was simpler, but not invariant under the cubic point group; so although it does correspond to $N_{CS}$ in the continuum limit it is contaminated by nonrenormalizeable operators already at dimension 5, whereas the symmetric definition is only contaminated at dimension 6.

2.4 Chemical potential for $N_{CS}$

A particularly efficient way to study the rate of topology change in the symmetric electroweak phase is to apply a chemical potential $\mu$ for $N_{CS}$ and to measure $\dot{N}_{CS}T/\mu$. To do so it is necessary to develop a definition for the chemical potential, since there is no good lattice analog to the continuum definition of $N_{CS}$. The most naive solution is that a chemical potential modifies the dynamics to

[TABLE]

which look just like the corresponding continuum equations. However, this does not work, because the $B$ term in the second equation does not preserve the Gauss constraint, ie

[TABLE]

(where $B_{i}^{a}(x-i)$ should be parallel transported as an adjoint object along the link leading to $x$). The corresponding continuum expression, $D_{i}B_{i}$, is zero, by the Bianchi identity; but the proof of the Bianchi identity depends on continuity of the $A$ fields, and the notion of continuity is lost on the lattice. Eq. (32) is zero in abelian gauge theory. For SU(2) in the low temperature or small $a$ limit, the violations are higher order in temperature than $B$, but nevertheless the failure of Eq. (32) to equal zero means that the proposed implementation will spoil the Gauss constraints, exciting unphysical modes. As I discuss in [9], this makes the results unreliable.

The solution is to take very literally the constraints. What they mean is that we are not really living in the phase space of all $E$, $U$, but on the constraint subspace, where the $E$ satisfy Gauss’ law. The only actual dynamical fields are the linear combinations of $E$ which are orthogonal to the Gauss constraints; if we could, we would express the system and its evolution in terms only of these.

Consider the most obvious basis of electric fields, where a basis element is a specification of a lattice point $x$, a direction $i$, and a Lie algebra element $a$. I will write this basis as $\{E_{\alpha}\}$. This basis has the important property that, in terms of the inner product defined by $2H_{E}$, it is orthonormal. I write the Gauss constraints in this basis as $G_{\alpha}=c_{\alpha\beta}E_{\beta}=0$, where note that the first index of $c$ only runs over the $3N^{3}$ Gauss constraints ($N$ the number of lattice points on a side of a cubic lattice), while the second runs over the $9N^{3}$ independent electric field basis elements. (Sums on repeated greek indicies are always implied.) Note also that some of the coefficients $c_{\alpha\beta}$ are $U$ dependent, because the $E$ fields on the links leading into a point must be parallel transported to the point. Now assuming that all Gauss constraints are independent222in the abelian theory with periodic boundary conditions this is in fact not the case; there is precisely one linear dependence between the Gauss constraints. But for generic $U$ the Gauss constraints in the nonabelian theory are independent. we can find a different orthonormal basis in which $c_{\alpha\beta}$ is nonzero only for $\beta\leq 3N^{3}$, so the basis can be partitioned into two subsets, the $E^{c}$, which are an orthonormal basis for the Gauss constraints, and the $E^{*}$, which are the dynamical degrees of freedom. Now the definition of the discrete $dN_{CS}/dt$ can be written as

[TABLE]

Since we are only interested in the situation in which all $E^{c}$ are exactly zero, we may as well define $dN_{CS}/dt$ as

[TABLE]

and hence the action of the chemical potential should be to modify the $E^{*}$ in proportion to the parallel component of $B$, but to leave the $E^{c}$ alone. In other words, if we write the orthogonal transformation between the $E^{c},E^{*}$ basis and the $E_{\alpha}$ basis as

[TABLE]

where orthogonality ensures

[TABLE]

then the update rule should be

[TABLE]

or equivalently

[TABLE]

This rule identically preserves the Gauss constraint while reflecting the meaning we want for a chemical potential for $N_{CS}$, that it changes electric fields in proportion to the parallel magnetic field component. Particularly important, it will satisfy the property that the change in the system energy

[TABLE]

is $-\mu N_{CS}$, as it should be. This allows a check of the algorithm by comparing the accumulated energy to $N_{CS}$. (Note that this is where it is essential that the algorithm without chemical potential should preserve system energy, in the mean, to high precision.)

The problem with this update rule lies in finding the orthogonal transformation to the basis which splits into constraints and dynamical fields, which requires diagonalizing the matrix $c_{\alpha\beta}$. Even if the diagonalization procedure could be done efficiently, using it via Eq. (37) or (38) would take $O(N^{6})$ computations, which is prohibitive. However, notice that Eq. (38) shows that the correct algorithm is to perform the naive implementation, but then to orthogonally project back to the constraint surface. This can be done approximately, but very accurately, by an algorithm with $O(N^{3})$ computations per update.

The idea is the following [9]; the most general way of moving orthogonally to the constraint surface is to change $E_{\alpha}$ by

[TABLE]

with $\kappa$ completely arbitrary. (It is essential that the motion be strictly orthogonal to the constraint surface so the dynamical fields are not changed.) The right choice of $\kappa$ is the one which will change $E$ so as to satisfy all Gauss constraints. A reasonable first guess, which is guaranteed to reduce the sum of squares of Gauss constraint violations, is $\kappa_{\alpha}=-\gamma c_{\alpha\beta}E_{\beta}$, $0<\gamma<1/6$, which is a step in a dissipative or quenching algorithm for the Gauss constraints. (The condition $\gamma<1/6$ is demanded by stability.) A crude approximate projection algorithm is to do one application of this quench at each time step. A more refined algorithm is to guess that the correction should be almost the same as in the previous timestep, and to first apply $\kappa_{\beta}=m\kappa_{\beta,{\rm previous\;step}}$, where $m<1$ ensures algorithm stability and $\kappa_{\beta,{\rm previous\;step}}$ is the sum of all corrections applied at the previous timestep. Then the dissipative step is applied twice, with different stepsizes $\gamma_{1}$ and $\gamma_{2}$. I find that using the choices $m=1-1.6\Delta t$, $\gamma_{1}=5/24$, $\gamma_{2}=5/48$ is stable and highly efficient. (Note that when two applications of the quench algorithm are performed in a row, the stability criterion changes, see [9].) Applying this approximate projection algorithm requires about as many floating point operations as determining the $DF_{\alpha}$, so it is not a great strain on algorithm efficiency. It is also close enough to an exact projection that the residual failure has virtually no effect on the results [9].

2.5 Thermalization

Finally I discuss the thermalization of the lattice Yang-Mills system developed in this section. What a thermalization algorithm needs to do is to draw a set of $E$, $U$ from the constraint submanifold of the phase space of $E$, $U$, choosing with weight $\exp-\beta_{L}H$ times the natural measure on the constraint surface. This definition of $\beta_{L}$ corresponds to $\beta_{L}=4\beta_{\rm continuum}/(g^{2}a)$, with $a$ the lattice spacing in the same physical units as $\beta_{\rm continuum}$. (The factor of $g^{2}$ is the usual Yang-Mills wave function normalization and the factor of 4 arises because I work in terms of $\tau^{a}$ rather than $\tau^{a}/2$, which are used in the definition of the continuum fields.) It is possible to perform the thermalization for a general constraint manifold, using a set of Langevin equations [7]; but by using the special properties that the subspace of $E$ fields at a fixed value of $U$ is a vector space, that the constraints are linear, and that the Hamiltonain is separable into a part which depends only on $U$ and a part which is quadratic in $E$ at fixed $U$, it is possible to find a simple and efficient thermalization algorithm [9].

The idea is as follows; while the $U$ dependence of the Hamiltonian, $H_{U}$, is quite nontrivial, all of its complexities appear in the equations of motion; and the equations of motion preserve the thermal probability distribution, while mixing excitation between the $E$ and $U$ fields. Since we must implement the equations of motion anyway, it is best to use them to thermalize the system. All we need is a thermalization algorithm for the $E$ fields at fixed $U$; we then repeatedly draw the $E$ fields from the fixed $U$ thermal distribution, evolve the system to mix the thermalization between $E$ and $U$ fields, and throw away the $E$ fields and repeat. The idea is identical in spirit to the “molecular dynamics” Monte-Carlo algorithm of Euclidean lattice gauge theory. The only complication here is that the $E$ fields must satisfy constraints.

The thermalization of the $E$ fields is possible because the $E$ fields at fixed $U$ form a vector space, the constraints are linear, and $H_{E}$ is quadratic in $E$. In terms of the orthogonal basis which partitions into $E^{*}$ and $E^{c}$, the Hamiltonian is just $H_{E}=E^{*}_{\alpha}E^{*}_{\alpha}/2+E^{c}_{\alpha}E^{c}_{\alpha}/2$; we just need to choose each $E^{*}_{\alpha}$ from the Gaussian distribution

[TABLE]

and set each $E^{c}$ to zero. The easiest way to do this, given the difficulty of finding the basis which partitions into $E^{*}$ and $E^{c}$, is to set both equal to the Gaussian distribution, and then orthogonally project to the constraint surface. To get the initial Gaussian distribution, I simply draw the $E_{\alpha}$ from the same Gaussian distribution; because the desired basis is gotten by orthogonal transformation from this natural basis, the resulting distribution will also be Gaussian in the basis of $E^{*}$ and $E^{c}$. The orthogonal projection is performed by using the quenching step of the algorithm of the last section repeatedly until the projection is complete. I find that order 100 applications reduce the Gauss constraint violation to the level of roundoff error (for single precision numerics).

In practice I find this algorithm to be very efficient. For instance, the total energy of the system approaches the thermal average value by a factor of $1/2$ on each $E$ randomization, if the Hamiltonian evolutions are taken to be of length $\simeq\beta_{L}$; within 7 or 8 applications the energy has reached the level of random fluctuation expected in a system of finite size, and the values of Wilson loops equal their thermal averages to within expected thermal fluctuation. Two or three applications are enough to re-randomize a thermalized distribution.

This supplies all necessary equipment for the numerical study of classical 3-D Yang-Mills theory under the minimal action.

3 Improved Hamiltonian

In this section I construct an improved Hamiltonian for the classical lattice Yang-Mills theory. I also present an improved definition of $dN_{CS}/dt$ and discuss how to implement a chemical potential for this quantity, and how to thermalize the system.

3.1 a false start

The simplest guess for an improved Hamiltonian takes inspiration from the example of scalar theory presented in the introduction and replaces the gradient term, $H_{U}$, with an improved term,

[TABLE]

which is known from the study of Euclidean lattice gauge theory [13]. Here $\sum_{\sqsubset\!\sqsupset}$ refers to a sum over all $1\times 2$ rectangular plaquettes, of which there are $6N^{3}$. The evolution of the system is the same as in the previous section, except that the definition of $DF_{i}^{a}(x)$ changes to

[TABLE]

where $12\sqsubset\!\sqsupset$ refers to the 12 rectangular plaquettes containing the link $U_{i}(x)$, traversed so as to begin and end at $x$ and to contain $U_{i}(x)$ and not $U^{\dagger}_{i}(x)$.

To see that this proposal fails, it is useful to consider a sinusoidal excitation in a single Lie algebra element, and to compute the dispersion relation. Consider the infinite volume limit of the lattice and start from the initial conditions

[TABLE]

where $\epsilon$ is infinitesimal. (The coordinate in the argument of sin is the center of the link.) The energy per site of this arrangement is $\sum_{i}(E_{i})^{2}\epsilon^{2}/2$. An infinitesimal time $\delta t$ later some of this energy has appeared as magnetic energy; the fraction is $\delta t^{2}\omega^{2}$, and $\omega$ gives the dispersion relation. At time $a\delta t$, the $U$ fields are

[TABLE]

and a lengthy but straightforward calculation gives the energy from square plaquettes as

[TABLE]

(where I introduce the shortand notation $s_{i}\equiv 2\sin(ak_{i}/2)$), the energy from rectangular plaquettes as

[TABLE]

and the Gauss constraint as

[TABLE]

by adding the square of the Gauss constraint, which is zero, to the energy in magnetic fields, I readily find that the dispersion relation is

[TABLE]

The contribution to $\omega^{2}$ from the first line is $k_{i}^{2}$ plus zero times terms of form $a^{2}k^{4}$, plus $O(a^{4}k^{6})$, and is therefore the correct dispersion relation, up to corrections which vanish as $a^{4}$ in comparison to the leading term. The other terms, which vanish when $E_{1}s_{1}=E_{2}s_{2}=E_{3}s_{3}=0$, spoil the dispersion relations by introducing polarization dependent $O(a^{2}k^{4})$ contibutions, so the evolution is not “improved”.

An even simpler way of seeing that the Hamiltonian is not “improved” is to consider the system’s thermodynamics, which are determined by integration of the fields over the weight $\exp-\beta_{L}H$. Here it is convenient to follow [8] and introduce Lagrange multipliers $A_{0}^{a}(x)$ for the Gauss constraints, so the Hamiltonian is

[TABLE]

where again parallel transport of $E^{a}_{i}(x-i)$ to the point $x$ is to be performed along the link leading into $x$. In terms of the notation of Sec. 2.4 I can write this as

[TABLE]

If we are only interested in the values of gauge invariant operators constructed from the link matricies then we can perform the integral over the $E$ fields, which is now Gaussian, and gives

[TABLE]

Note that $A_{0\beta}c_{\beta\alpha}E_{\alpha}$ is (performing the sum on $\beta$ and expressing the sum on $\alpha$) $\sum_{x,i}E_{i}^{a}(x)(A_{0}^{a}(x)-A_{0}^{a}(x+i))$, parallel transport of the forward $A_{0}$ implied, as usual. So, the $A_{0}$ expression above can be untangled to read

[TABLE]

which is the lattice action of the finite temperature quantum theory in the approximation of dimensional reduction, with zero bare Debye mass. It has the improved plaquette action, but the unimproved scalar ($A_{0}$) gradient term, as a comparison with the scalar theory presented in the introduction will show. Hence the Hamiltonian has not been “improved,” even at the level of thermodynamics, and this failure has nothing to do with the choice for an improved magnetic term.

3.2 Hamiltonian and Gauss constraints

The problem with the approach above is that, while I have removed nonrenormalizeable operators such as $a^{2}\sum_{ij}(D_{i}F_{ij})^{a}(D_{i}F_{ij})^{a}$, I have left in the nonrenormalizeable operator $a^{2}\sum_{i}(D_{i}E_{i})^{a}(D_{i}E_{i})^{a}$. An easy way to see why is to recall how the minimal Hamiltonian system could be interpreted as arising from an action on a lattice in both space and time, consisting of all $1\times 1$ plaquettes, in space and time. The $E^{2}$ term in $H$ is just the sum over the $1\times 1$ plaquettes with time links. In the previous section, I included $1\times 2$ plaquettes only if each direction were in space, but continued to use a kinetic $E^{2}$ term derived only from the square plaquettes. Instead I should include rectangular plaquettes here too. The result is that the continuum time Hamiltonian should be

[TABLE]

where the first term comes from square plaquettes, the second comes from plaquettes which stretch 2 units in time, and the third comes from plaquettes which stretch 2 units in space, and I have already taken the small time step limit. As usual, $E_{i}^{a}(x+i)$ must be parallel transported along the straight path. The second term should be absorbed in the first, changing its coefficient to $4/3$. I have retained my earlier definition of $E_{i}^{a}(x)$, so in the $A_{0}^{a}=0$ gauge the update rule for the $U$ is

[TABLE]

The Hamiltonian is conveniently summarized by defining the operator $M$,

[TABLE]

(parallel transports again implied), in terms of which the Hamiltonian is

[TABLE]

The Gauss constraint is derived as before, by varying $A_{0}(x)$; it is the sum of all plaquettes which traverse the time direction and terminate at $x$, with appropriate signs:

[TABLE]

or in somewhat more compact notation,

[TABLE]

Now I will address the question of whether this Hamiltonian is improved in the sense of replicating the continuum evolution for smooth fields without contamination from dimension 6 operators. There are two general types of dimension 6, nonrenormalizeable operators which can appear in cubic point group invariant lattice implementations of Yang-Mills theory; there are those with two powers of the field strength tensor, such as $\sum_{ij}(D_{i}F_{ij})^{a}(D_{i}F_{ij})^{a}$, and terms which contain three powers of the field strength, like $f_{abc}F_{ij}^{a}F_{jk}^{b}F_{ki}^{c}$. In the latter term, the indicies of the $F$ must differ, so that they represent curvatures in different planes, since $f_{abc}F_{ij}^{a}F_{ij}^{b}$ automatically vanishes. Since all the plaquettes used in Eq. (59) are planar, terms involving three powers of the field strength will not appear at dimension 6. I can therefore restrict attention to the other sort of term, which appear already in the Abelian theory and which modify the dispersion relations of weak disturbances which vary sinusoidally in space. It is sufficient, then, to check whether the dispersion relation for $\omega^{2}$ is free of $k^{4}a^{2}$ terms in the small $k$ expansion. It is again sufficient to start a sinusoidal excitation with $U=I$ everywhere and see how quickly the energy is transferred into the magnetic term of the Hamiltonian. For

[TABLE]

the electric energy is $\epsilon^{2}/2$ times

[TABLE]

with the first term from the single $E$ term in the Hamiltonian and the second term from the two $E$ term in the Hamiltonian. Again I use the notation $s_{i}=2\sin k_{i}/2$.

The new Gauss constraint is

[TABLE]

and the magnetic energy at time $a\delta t$ is still

[TABLE]

as in the previous subsection. Squaring the Gauss constraint, which equals zero, and adding it to the magnetic energy term, I find that the oscillation frequency is

[TABLE]

The first term here is

[TABLE]

which is the usual dispersion relation for an improved, scalar Laplacian. It contributes no $a^{2}k^{4}$ terms to $\omega^{2}$, as can be readily verified by expanding the cosines as power series. The second term is polarization dependent, but has 6 powers of $2\sin ak/2$, so it only starts to contribute at order $a^{4}k^{6}$, and the dispersion relation is therefore free of the unwanted $a^{2}k^{4}$ terms, and the Hamiltonian is free of dimension 6 nonrenormalizeable operators. Also note that the new Gauss constraint demands that $E$ be transverse to a higher order in $k$ than the unimproved Gauss constraint does.

I can also check the thermodynamics as before. Introducing $A_{0}(x)$ as a Lagrange multiplier for the constraint, the thermodynamics is governed by the probability distribution $\exp-\beta_{L}H$, with

[TABLE]

which, on performing the Gaussian integration over $E$, becomes

[TABLE]

The $A_{0}$ term untangles to

[TABLE]

parallel transports again implicit. The derivative term for $A_{0}$ is precisely the improved derivative term for the scalar theory discussed in the introduction, made covariant. Hence the Hamiltonian does give improved thermodynamics. This also points out a relation between the improvement scheme for the plaquette action and for the discrete Laplacian.

3.3 Numerical evolution of the improved system

Rather than derive continuum time equations of motion from the Hamiltonian, I will skip to finding a discrete time evolution algorithm. Naively, I should base this on the action

[TABLE]

which includes plaquettes which go two steps in time. These plaquettes prevent a nonrenormalizeable operator of form $(\Delta t)^{2}(D_{0}F_{0i})^{a}(D_{0}F_{0i})^{a}$ and should make the evolution fourth order in stepsize. This is unnecessary if we take $\Delta t<<1$, and is also undesirable since the inclusion of plaquettes which stretch two steps in time makes the numerical evolution unstable. The reason is that the update rule derived by varying one link then determines $U$ based on the 4 previous time slices, rather than just 2 as the minimal algorithm did; since the Yang-Mills equations are second order, this is too much state information, which opens the possibility (which is realized) of an unstable growing mode. Instead I choose $\Delta t$ small enough to make those plaquettes unnecessary, and choose

[TABLE]

As before, I will define $E_{i}(x,t+\Delta t/2)$ as the square plaquette, so in the gauge $U_{0}(x)=I$,

[TABLE]

The Gauss law is that the sum of Lie algebra traces of all plaquettes containing a time link be zero, which in the same gauge is

[TABLE]

or

[TABLE]

Implicitly parallel transporting using the connection at time $t$, this is

[TABLE]

whereas if the parallel transports are performed using the connection at time $t+\Delta t$, this is

[TABLE]

The equations of motion, derived by varying one spatial link, are also somewhat more complicated. In terms of the definition of $DF_{i}^{a}(x)$ appearing in Eq. (43), and implicitly parallel transporting objects using the connections at time $t$ and the straight line paths, the update rule is

[TABLE]

Writing this equation in the convenient shorthand

[TABLE]

one can see that the Gauss constraint is

[TABLE]

and that the difference between the value of this quantity at times $t+\Delta t/2$ and $t-\Delta t/2$ will be

[TABLE]

which vanishes identically, because each plaquette which appears in the evaluation of one term appears with opposite sign in the evaluation of another term, as before. Hence the Gauss constraint will be preserved identically by the evolution.

Eq. (79) does not give the Lie algebra content of the electric fields, $(1/2){\rm Tr}-i\tau^{a}E_{i}(x)$, directly, but gives the result of a nondiagonal, nonlinear transformation on them; to implement the update it is necessary to invert this transformation. The inversion problem is not too bad (for $\Delta t\ll 1$), because the expression on the lefthand side of Eq. (79) is close to diagonal in the $E$ fields and the inversion can be performed iteratively, ie

[TABLE]

can be solved by making a guess for $E$, substituting into the righthand side, and getting a better guess for $E$. I have written the formula so that the righthand side depends almost only on $E$ at the points 1 space ahead and behind of the point of interest, so the convergence rate is doubled by alternately updating even and odd lattice points.

I have implemented this algorithm numerically. At single precision, choosing the initial guess for $E_{i}(x,t+\Delta t/2)$ to be $E_{i}(x,t-\Delta t/2)$, it is only necessary to iterate over the lattice 4 times before the remaining correction comes on order roundoff error. If the inversion were conducted completely to infinite precision then the algorithm would not gain or lose energy in the mean, because the update rule is time reversal invariant and only depends on the minimal amount of state information. For a single precision implementation, because the iteration converges from above, the roundoff errors systematically boost the energy, so that it rises by about 2 parts in $10^{8}$ per time step; but this level is small enough to be of no concern if the total length of a simulation is kept under $t\simeq 5000$. By going to double precision and making the timestep very small to reduce the fluctuations in system energy I was able to verify that the algorithm does preserve energy in the mean. I also observe that violation of the Gauss constraint is only generated by roundoff error, as in the evolution with the unimproved action. As a final check I verified that the implementation produces the expected dispersion relations when the initial state is a weak sinusoidal excitation in one Lie algebra direction.

3.4 Magnetic field and $N_{CS}$

The next step in the program is to find the appropriate definition of $dN_{CS}/dt$ for the improved Hamiltonian. The approach is similar to the unimproved case, where $\int(dN_{CS}/dt)dt$ is defined as

[TABLE]

where the first sum is over all points in the 4 dimensional grid, the second is an average over the 4 plaquettes in the $jk$ plane ($ijk$ a cyclic permutation) which touch the point, and the last is a sum over the 4 plaquettes which proceed one step forward or backward in the time direction and a step forward or backward in the $i$ direction from the point. That is, it is an average over the electric fields leading into and out of the point, taken a half timestep before and after the spacetime point. In the previous section I rearranged this sum into a sum over the electric fields of a contribution coming from each end of the link on which the electric field resides, which is how the sum over 8 square plaquettes arose.

To modify this definition so as to remove dimension 6 contamination, I need to replace the sum over square plaquettes with a sum on square and $1\times 2$ rectangular plaquettes. Again, because the plaquettes are planar and the improvement need only correct dimension 6 errors (dimension 4 in the electric or magnetic field), terms cubic in the field strength will not arise and I need only consider weak sinusoidal fields in one Lie algebra direction. In a small $k$ expansion, the choice of spatial plaquettes should not contain any $k^{3}$ term.

The contribution from the 4 square plaquettes around a point $x$ and in the $1,2$ plane for $U_{i}(x)=1+\epsilon i\tau^{a}\sin(k\cdot x+k_{i}/2)$ is

[TABLE]

where $c_{i}=\cos k_{i}/2$. The contribution from the 4 rectangular plaquettes which have midpoints at $x$ is twice this, and the choice of whether to use square plaquettes or midpoints of rectangular ones is ours (though it may become forced when improvement is carried up to higher order). The contribution from the 8 rectangular plaquettes which have a corner at $x$ is

[TABLE]

The correct combination should produce no $k^{3}$ term. The requirement is solved by using $5/3$ the sum over square plaquettes (or $5/6$ the sum over rectangular ones with midpoints at $x$) plus $-1/6$ the sum over rectangular plaquettes with corners at $x$. For the electric field it is not necessary to use rectangular plaquettes which stretch two units in the time direction, because the field variation in this direction is suppressed by a factor of $(\Delta t)^{2}$, and the right combination at the point $x,t$ (making parallel transports using the links at time $t$) is

[TABLE]

This differs only by terms of order $(\Delta t)^{2}$ from the same expression with the trace taken over each $E$ in the two $E$ terms, rather than over the product. This is because the order $\Delta t$ correction arises from a commutator term which enters with opposite sign in the $t-\Delta t/2$ and $t+\Delta t/2$ contributions. Since there are already inevitable $O((\Delta t)^{2})$ errors in the evolution of the system, I am free to choose the more convenient definition without worrying about the order $(\Delta t)^{2}$ change it will make in the results. Choosing to trace each $E$ separately, I can rearrange the definition to be

[TABLE]

where all definitions should be clear from earlier usage and parallel transports of $F_{jk}^{a}$ are implied along the straight line paths to the point $x$.

This gives a usable definition for the evolution of $N_{CS}$ which is free of dimension 6 contamination.

3.5 Chemical potential for $N_{CS}$

Now I should apply a chemical potential for $N_{CS}$ in order to study its motion due to infrared physics. In analogy with the unimproved Hamiltonian case, the chemical potential term should modify the evolution of the $E$ fields in a way which makes the overall energy change by $-\mu N_{CS}$. In the unimproved case this was accomplished by making

[TABLE]

so the change in the electric energy was

[TABLE]

Then I modified the definition so that the linear combination of $E$ fields representing a Gauss constraint violation were not excited after all, by orthogonally projecting to the constraint surface. This did not change the energy because these modes had started out zero anyway.

In the improved case, it is easiest to begin by thinking about the continuum time system. The change in $N_{CS}$ is

[TABLE]

but the shift in the energy is $\delta E_{\alpha}M_{\alpha\beta}E_{\beta}$, so $\delta E$ should be

[TABLE]

before worrying about preserving the Gauss constraint. Note that the evolution equation for $E$ already includes a term like $-M^{-1}_{\alpha\beta}DF_{\beta}dt$, so the magnetic field appears in the same place as $DF_{i}^{a}(x)$. This describes a possible implementation in the discrete case as well, namely that Eq. (79) is modified by replacing $-DF_{i}^{a}(x)$ by $-DF_{i}^{a}(x)-\mu B_{i}^{a}(x)/2\pi^{2}$.

Another way of arriving at this result is to note that the Hamiltonian defines a metric on the space of $E$, and that in the improved Hamiltonian the natural basis $E_{\alpha}$ is no longer orthogonal. Rather, the orthogonal basis is $(M^{1/2})_{\alpha\beta}E_{\beta}$. In terms of this basis, the definition of $N_{CS}$ is

[TABLE]

and so the change in the electric fields should be

[TABLE]

which leads to Eq. (94).

Now although the new definition of the magnetic field should be free of the leading order of operators responsible for it not satisfying the Bianchi identity, so the satisfaction of the Gauss constraint should be better in the infrared, the ultraviolet violation of the Gauss law is just as bad as in the unimproved case, and it is again necessary to change the definition of the action of a chemical potential to prevent the violation of the Gauss constraint. In analogy with the unimproved case, the correct modification is to orthogonally project back to the constraint surface. Again, orthogonal should mean with respect to the metric defined by $H_{E}$, so the most general motion orthogonal to the constraint surface is

[TABLE]

For small $\kappa$, the former changes the system energy by

[TABLE]

which is suppressed by a power of the adherence to the Gauss constraint, whereas the latter changes the system energy by

[TABLE]

which is nonzero even if the Gauss constraint is satisfied; so the modification to $E$ is not orthogonal to the constraint surface.

A nice way to understand this conclusion is to remember that $M^{1/2}E$ are orthogonal; the Gauss constraint is

[TABLE]

so the most general correction to the electric field which is orthogonal to the constraint surface should be

[TABLE]

which is the same as Eq. (97).

I can now find an approximate projection algorithm in strict analogy with the unimproved Hamiltonian case. A step in a dissipative or quench algorithm for the orthogonal projection is

[TABLE]

where stability now demands $\gamma<1/8$. This step is orthogonal to the constraint surface, and repeated application is guaranteed to be a projection; a single application at each step is enough to keep the departure from the constraint surface under control. A better algorithm, again in analogy with the unimproved system, is to make an initial correction with $\kappa$ equal to $m<1$ times the total $\kappa$ of the previous step, and then to apply the dissipative algorithm twice, with different stepsizes.

Note that Eq. (97) is of exactly the same form as in the unimproved case. To understand why, recall the origin of the Gauss constraint; it is the Lagrange equation for the $A_{0}$ field (or the time direction link matrix), and when it is not satisfied it means we have inadvertently not set $A_{0}$ to zero; since the electric field only corresponds to the time derivative of the link matrix when $A_{0}$ is zero, we should restore $A_{0}$ to zero, and Eq. (101) is precisely how the system changes as we modify $A_{0}$ by $\kappa$.

This latter observation explains how to project to the constraint surface in the discrete time case. The modification of the electric fields should take the form of a time dependent gauge transformation at each point, with magnitude chosen to restore the Gauss constraint everywhere. In other words the correct orthogonal projection to the constraint surface is a choice of a group element $g(x)$ at each point $x$ which modifies the $E$ matricies through

[TABLE]

The $g(x)$ should be chosen such that

[TABLE]

(parallel transports at time $t$ implied) is zero. Again the job of finding the right $g$ is intractable, but a reasonable guess (the basis of the dissipative algorithm) is

[TABLE]

and a better guess is to use $m$ times the Lie algebra element used in the previous time step, and then to apply the above dissipative step twice with different coefficients $\gamma$. The stability condition is now $-1/3<(1-\omega\gamma_{1})(1-\omega\gamma_{2})<1$ for all $0<\omega<16$, and I find that a good choice is $m\sim 1-1.6\Delta t$, $\gamma_{1}=5/64$, $\gamma_{2}=5/32$. Using these values, the Gauss constraint is as nearly satisfied as in the unimproved Hamiltonian, and the relation $\Delta H=-\mu\Delta N_{CS}$ is satisfied to $<1\%$ (after energy gain due to roundoff errors are taken into account; I checked this in a double precision implementation, where there is no energy gain from roundoff error), which should assure that the systematic errors in the determination of the linear response to a chemical potential, due to not using an exact orthogonal projection, are of the same size [9].

This provides a satisfactory implementation of chemical potential in the improved Hamiltonian system.

3.6 Thermalization

Next I extend the thermalization algorithm developed for the unimproved Hamiltonian to the improved one. I will begin by showing how to thermalize the system in the continuum time case. The technique does not strictly work for finite $\Delta t$, but the failures are of order $(\Delta t)^{2}$, and since the time evolution and the definition of $N_{CS}$ already introduce systematic errors of at least this size I will ignore this problem and find a discrete implementation which is correct up to these stepsize ambiguities.

In the continuum time case, I want to thermalize the system with Hamiltonian given in Eq. (59), subject to the Gauss constraints given in Eq. (61). The subspace of momenta at fixed $U$ is an inner product space with $2H_{E}$ the inner product, and the Gauss constraints are linear, so I can implement the thermalization discussed in Subsec. 2.5, modified only to account for the new inner product; in other words I need only find a thermalization algorithm for the fixed $U$ electric field part of the Hamiltonian and then use the equations of motion to stir the thermalization between electric and magnetic modes.

If we had an orthonormal basis which partitioned into Gauss constraints and dynamical degrees of freedom, we would choose the dynamical degrees of freedom to be Gaussian distributed and the constraints to be zero, and the electric fields would be thermalized. Equivalently I could Gaussian distribute both of them and then orthogonally project to the constraint surface. The continuum time orthogonal projection algorithm is given already in the previous section, and a Gaussian distribution in any orthonormal basis is Gaussian in any other, so the only remaining problem is to find an orthonormal basis. In the previous section we saw that the basis $M^{1/2}E$ is orthonormal. Therefore, a thermalization of the electric fields consists of choosing the $M^{1/2}E$ Gaussian, inverting $M^{1/2}$ to find $E_{\alpha}$, and applying the orthogonal projection by repeated application of the continuum time dissipative algorithm of the last section.

The inversion can be carried out as follows: note that $M$ is $(7/6)I$ plus a small off diagonal term, $M_{\alpha\beta}=7/6(\delta_{\alpha\beta}+m_{\alpha\beta})$; write

[TABLE]

in other words, expand $M^{-1/2}$ in a Taylor series. The absolute magnitude of the largest possible eigenvalue of $m$ is 1/7, so the Taylor series converges rather well; and the inefficiency of the algorithm does not matter because it is only called occasionally, with several (order 100) lattice updates inbetween.

It is not completely trivial to carry this algorithm over directly to the discrete time setting. For one thing it is not clear that randomizing $E_{i}(x,t+\Delta t/2)$ is the appropriate thing to do. If we did, the thermalization algorithm would not be quite time reversal invariant, as ${\rm Tr}-i\tau^{a}E_{i}(x,t+\Delta t/2)$ would then be uncorrelated with $DF_{i}^{a}(x,t)$, but backwards evolving, we would find that ${\rm Tr}-i\tau^{a}E_{i}(x,t-\Delta t/2)$ was correlated with $DF_{i}^{a}(x,t)$. Also, the $E$ are matricies in SU(2), not Lie algebra elements, and while there is a mapping from SU(2) to the Lie algebra, given by $1/2{\rm Tr}-i\tau^{a}$, which is one to one from half of SU(2) and not difficult to invert in this patch, the determinant of the map is not 1, ie it does not preserve the measure of SU(2). Also, the energy depends not on $1/2$ the square of the Lie algebra element, but on $1-1/2{\rm Tr}$ of the group element, which, in terms of $E_{i}^{a}\equiv 1/2{\rm Tr}-i\tau^{a}E$, is $1-\sqrt{1-E^{a}E^{a}}\neq E^{a}E^{a}/2$. In terms of the Lie algebra element, nonquadratic terms will appear in the Hamiltonian, both from the measure and from the way the group element is traced to get the energy. Therefore the thermalization algorithm I have used is, strictly speaking, wrong.

However, I can define a Lie algebra element “electric field” at integer time steps, for which all of these ambiguities would appear as order $(\Delta t)^{2}$ or $(\Delta t)^{2}/\beta_{L}$ corrections, and can in practice be ignored. (It is not even clear that there is a meaning for temperature beyond leading order in $(\Delta t)^{2}$ in the discrete time system.) The definition is

[TABLE]

(all parallel transports to be conducted with respect to the links $U_{i}$ at time $t$), from which it is simple to continue the time evolution of the system. The Gauss constraints on either side of time $t$ will be satisfied provided that $E_{i}^{a}(x,t)$ obeys

[TABLE]

I will thermalize this $E_{i}^{a}$ by drawing it from the thermal distribution of the Hamiltonian $E_{\alpha}M_{\alpha\beta}E_{\beta}/2$. The algorithm is precisely the continuum time algorithm presented above.

To check that this thermalization algorithm is reasonable, in the sense that it gives thermodynamic properties which change only very weakly with stepsize when the stepsize is small, I have evaluated the magnetic energy and Wilson loops at different stepsizes. Some results for a $16^{3}$ lattice at lattice temperature $\beta_{L}=5$ are presented in Table 1. For each column I have used the algorithm to thermalize an initial configuration and then continued to implement the algorithm for 80 $E$ field randomizations, recording each observable once each $E$ field randomization. The error bars are from statistics, corrected for correlations because one $E$ field randomization is insufficient to fully rethermalize the system. There are clear stepsize effects which are consistent with the expected $(\Delta t)^{2}$ behavior. The changes are significant for the two larger stepsizes in the table, so $\Delta t$ cannot be taken to be this large; but the zero $\Delta t$ extrapolated magnetic energy differs from the $\Delta t=0.05$ value by only $0.5\%$, which will lead to a systematic error in Lyapunov exponents of order $0.5\%$ and in the motion of $N_{CS}$ under chemical potential of order $1.5\%$, which will be less than statistical measuring error. $\Delta t=0.05$ is therefore sufficiently small that the remaining stepsize systematics can be neglected. Note also that these systematics will apply to all measurements, so they will not affect any conclusions about the existence of a fine lattice spacing limit.

Based on this analysis, and because the evolution of the system will have stepsize errors of similar size, I have used $\Delta t=0.05$ everywhere in what follows, unless noted otherwise.

4 Numerical results

I have implemented the improved Hamiltonian system developed in the previous section numerically, and used it to study the motion of $N_{CS}$ under a chemical potential. It has become conventional in the field to quote the rate of $N_{CS}$ violation in the symmetric phase in terms of a dimensionless quantity $\kappa$. The linear response coefficient to a chemical potential is

[TABLE]

which in lattice units is

[TABLE]

By a detailed balance argument [15, 9] this is half the rate of $N_{CS}$ diffusion,

[TABLE]

measured in [8].

I have evaluated the rate of $N_{CS}$ motion under a chemical potential for several lattice spacings (lattice reciprocal temperatures), using in every case a chemical potential $\mu=2/\beta_{L}$, which previous work shows is easily small enough to fall in the linear response regime [9]. I used lattice sizes such that $N\geq 3\beta_{L}$ in most cases, because it has been shown that finite lattice size effects vanish already at $N=2\beta_{L}$ [8]. The evolutions were each of length $t=400$ in lattice units; for each lattice coarseness I used several such evolutions from independent initial configurations drawn from the thermal ensamble by the thermalization algorithm. The error bars are statistical; when there were several runs they were determined by the fluctuations between independent runs. I also computed them by assuming the statistical error comes from Brownian motion in $N_{CS}$, with a rate given by Eq. (111); the two estimates agree within error, so I used the latter estimate when there were few independent runs.

My results are presented in Table 2 and plotted, along with some results of [8, 9] for the unimproved Hamiltonian, in Figure 2. As the lattice becomes finer, $\kappa_{\mu}$ appears to converge to a lattice spacing independent result. The rate of convergence is slightly worse for the improved Hamiltonian than for the unimproved one, so the physics which establishes the rate of $N_{CS}$ diffusion depends on the Debye screening length being well shorter than the length scale of nonperturbative physics. In fact, from the data presented it is not entirely clear that the improved lattice results have converged to the fine lattice spacing limit; but the numerical demands rise as the 4’th power of $\beta_{L}$, and the improved Hamiltonian requires 5 times as many flops to evolve the same lattice 4-volume as the KS Hamiltonian , so I was unable to get useful data on finer lattices than those presented. I will assume here that the answer has converged at $\beta_{L}=10$, since even if the convergence rate is Debye mass limited, then the convergence should be as good here as in the unimproved system at $\beta_{L}\simeq 8$.

For the current level of statistics the results of the improved system are compatible with those found using the Kogut-Susskind Hamiltonian. However, the still substantial error bars leave open the possibility that the rates differ by order $5\%$ to $10\%$. To compare the rates to $1\%$ would require approximately 50 times the integration time used here, or about $2\times 10^{15}$ floating point operations, which is prohibitive.

To examine whether the two Hamiltonians really do produce different infrared dynamics, or whether I just have insufficiently refined lattices and insufficient statistics, I have measured another quantity, the maximal Lyapunov exponent, for the improved Hamiltonian system. The maximal Lyapunov exponent is defined as the rate at which two almost identical initial configurations diverge, as measured by some gauge invariant phase space metric. It is a better quantity to use to compare the improved and unimproved Hamiltonian, because it depends on less equipment than the measurement of $\kappa_{\mu}$, requiring only the evolution and thermalization algorithms, and because the statistics improve more quickly, rising roughly as $\sqrt{N^{3}t/\beta_{L}^{4}}$, whereas for $\kappa_{\mu}$ the statistics improved as $\mu\sqrt{\kappa_{\mu}N^{3}t/(\beta_{L}\pi)^{4}}/2$. It will therefore be possible to compare the maximal Lyapunov exponents of the two theories to order $1\%$ without excessive numerical demands.

I will base my approach on that of Krasnitz [7]. I draw an initial condition from the thermal ensamble and then perturb it slightly; then I allow the original configuration and the perturbed version to evolve under the classical equations of motion, periodically measuring the metric distance. At first the separation is dominated by transients; then there is an epoch of exponentially growing separation, and finally the separation saturates. It is the exponential epoch which I use. I repeat for several thermal initial configurations to determine and improve the statistics.

In my implementation I start by getting a thermal initial condition by applying my thermalization algorithm at single precision and a stepsize of $\Delta t=0.05$; I then convert to double precision, projecting the fields to the $SU(2)$ manifold and quenching the Gauss constraint violations to the level of double precision roundoff error. I then copy this initial configuration and evolve it one time step, both with and without a tiny ($10^{-13}$) chemical potential, to generate initial differences between configurations at the level of double precision roundoff error. I let the configurations evolve, tracking the metric distances

[TABLE]

as proposed in [7]. The choice of metric does not matter much as long as it is gauge invariant; the two choices above give almost exactly the same Lyapunov exponent for a run, and cannot be used as independent measurements. To eliminate transients I throw out all data until the metric distance has grown by a factor of 500; I then follow the divergence of the two paths, stopping when the metric distance is within a factor of 500 of saturation. This gives a very conservative fiducial period of evolution, over which the metric distance grows by a factor of $10^{9}$. (This is the advantage of using double precision.) I also use lattice sizes well larger than $\beta_{L}$ to prevent finite size effects. To get the Lyapunov exponent I make a linear regression fit of $\ln D_{E}$ or $\ln D_{M}$ against time for the data in this fiducial period. To test the technique, I repeated one of the measurements of Krasnitz, using the unimproved Hamiltonian and $\beta_{L}=4$, $N=12$; my value of $\lambda_{max}=.348\pm.004$ is within statistical error of the result presented in [7].

The technique thus tested, I repeated it for several lattice spacings, using the improved Hamiltonian; the results are presented in Table 3, which also contains the results of Krasnitz [7] for comparison. It is very clear from the table that $\beta_{L}\lambda_{max}$ has a good fine lattice limit, both for the improved and the unimproved Hamiltonian. It reaches its limit much faster in the improved case. The limit for the improved Hamiltonian differs from that for the KS Hamiltonian by about $5\%$, significant at about $5\sigma$, confirming that while the fine lattice limit of the infrared dynamics exists within a given lattice regulation, its value depends on the specifics of that regulation.

As a final test of stepsize dependence, I recomputed the $\beta_{L}=4$ datapoint using $\Delta t=0.1$ and $0.20$ (for both the thermalization algorithm and the evolution); the results for $\beta_{L}\lambda_{max}$ are $1.147\pm.008$ and $1.200\pm.010$, respectively. This agrees to within error with a quadratic dependence in stepsize. Back-extrapolating to zero stepsize changes the results in the table by less than statistical error.

Note also that the limit of the Lyapunov exponent in both regulations falls well below twice the damping rate of a plasmon at rest, which is [16]

[TABLE]

which would correspond to a value for $\beta_{L}\lambda_{max}$ of $1.408$.

5 Conclusion

An improved Hamiltonian can be constructed for the evolution of the classical Yang-Mills equations of motion on a lattice. Using it to compute the maximal Lyapunov exponent, I find a much faster approach to the continuum limit, but a value of that limit which differs by about $5\%$ from the value in the Kogut-Susskind Hamiltonian system. For the rate of $N_{CS}$ diffusion, the convergence to a fine lattice limit is apprarently no faster, and the answers agree to within error, but the statistical power of the test is limited, leaving open the possibility that, as in the case of the maximal Lyapunov exponent, there are lattice regulation artefacts which do not vanish in the $a=0$ limit.

Such an artefact presumably arises through interactions between the infrared modes, which determine the quanitites of interest, and the ultraviolet excitations, which move under lattice dispersion relations which differ between the two regulations, and are in neither case rotationally invariant or ultrarelativistic. This possibility arises because the “hard thermal loop” effects which the ultraviolet modes induce are linearly divergent in the classical thermal theory, so most of the contribution comes from the momentum region where lattice effects are significant; and because the functional form the hard thermal loop effects take is not strongly restricted by symmetries of the theory, so the lattice dispersion relations of the ultraviolet modes can cause the induced hard thermal loops to have the wrong dependence on $\vec{q}/q_{0}$. Bodecker et. al. have explicitly verified that this is the case for unimproved abelian Higgs theory [11]; in fact they show that the induced corrections to the infrared propagator are not even rotationally invariant. There is every reason to believe that the same thing happens in lattice Yang-Mills theory, in both the improved and unimproved cases.

What might at first seem surprising is that the overall magnitude of the hard thermal loop effects, which in thermal units diverges as $1/a$, apparently does not affect the infrared dynamical properties considered here, but the functional form of the hard thermal loops may. This is also true for one of the few infrared dynamical properties which is known for the quantum theory, the damping rate of a plasmon at rest; in its calculation the presence of hard thermal loops is essential, and their dependence on $\vec{q}/q_{0}$ determines the result, but their overall magnitude turns out to cancel [16]. Something similar seems to be happening for the maximal Lyapunov exponent. Since neither Hamiltonian considered to date produces the correct hard thermal loops, we cannot take either limit to be correct. The fact that they differ only by a modest amount, and that both are correctly accounting for the thermodynamic distribution of initial conditions, suggests that we can take the answers determined from either Hamiltonian as a rough $(\sim 10\%)$ estimate of the correct Lyapunov exponent and rate of $N_{CS}$ motion. However, to improve on this accuracy in the case of the rate of $N_{CS}$ motion it will be necessary not only to determine with greater resolution what the rate is in the Kogut-Susskind Hamiltonian system, but to check with the same resolution that it is the same in the improved system. If the answers do differ at the level of $5-10\%$, then to get the correct rate, some technique must be developed which incorporates or generates hard thermal loops of the correct form.

Acknowledgements

I would like to thank Alex Krasnitz, for sparking my interest in this matter and for helpful correspondence, and Chris Barnes, for useful conversations. I would also like to thank Neil Turok and Carl Edman, for donating computer resources.

All of the code used for the results presented here and in [9] is written in reasonably clear and well documented $c$, and is available upon request to my electronic address, guymoore@puhep1.princeton.edu.