the lattice gauge action agrees with its well known continuum counterpart1
SG=− 1 2g2
Z
d4xTrFµνFµν. (4.4)
It is also convenient to define the inverse lattice coupling β = 2N
g2 . (4.5)
Once the action has been defined, the expectation value of observables is given by the path integral
hOi= 1 Z
Z
d[U]Oe−SG[U], Z = Z
d[U]e−SG[U], (4.6) where d[U] = Πx,µdUµ(x) is defined in terms of the SU(N) invariant Haar measure dUµ(x).Z is called the partition function in analogy to the one in statistical physics.
Notice that because we are working on Euclidean space-time; the exponential factor is real and so it resembles the Boltzmann factor in a statistical ensemble average.
The resemblance suggests that one can borrow the widely developed techniques from statistical physics to compute expectation values on the lattice. In particular, one can use importance sampling and Monte-Carlo techniques to efficiently computehOi in Eq. 4.6.
4.2 Simulation algorithms
The way in which Eq. 4.6 is made manifest in the day to day calculations of lattice QCD practitioners, is through Monte-Carlo simulations. In simple terms, a sequence of gauge configurations U(i) is generated, such that
hOi= lim
n→∞
1 n
Xn i=1
O(U(i))≡O¯. (4.7)
The sequence is taken from a Markov chain generated by a stochastic process, with a given transition probability W(U →U0). The transition probability must be constructed so that it satisfies the following properties
1. Stability: P
Ue−SG[U]W(U →U0) = e−SG[U0]
2. Ergodicity: ∀U, U0,∃n <∞,s.t.Wn(U →U0)>0.
1Note that from now on, we will use this convention in contrast to the one in used in Eq. (2.1).
30 4.2. SIMULATION ALGORITHMS Clearly, when performing a simulation, one is limited by a finite amount of re-sources, and the summation in Eq. 4.7 has to be truncated for a large but finite value of n. In this sense, a more illuminating way to express the idea behind importance sampling is to write
Oˆ = 1 n
Xn i=1
O(U(i)) = ¯O+ O 1
√n
, (4.8)
which states that one can define an estimatorOˆ for O¯, whose error decreases as the number of measurements n grows. We will make this statement more precise in the following.
4.2.1 Autocorrelation times
In principle, the chain of configurations generated through the Monte-Carlo process are not all independent. This is evident from the fact that U(i) is generated from U(i−1) through W(U(i−1) → U(i)). Moreover, Eq. (4.8) is valid when the set of configurations U(i) is distributed according to the stationary probability density P(U) = e−SG[U]/Z. Let us discuss this issue first.
Assume the Markov chain starts from a given configuration U(1), then, the prob-ability distribution of the k-th configuration is given by Pk(U) = Wk−1U(1). It can be shown that Pk(U)approaches the stationary P(U)as
P(U) = lim
k→∞
Pk(U) + O(e−k/τexp)
, (4.9)
where τexp is the exponential autocorrelation time. Eq. 4.9 tells us that a given configuration U(i) obeys the right stationary probability density if i τexp. Thus, in a simulation, it is customary to start measuring only after the previous condition has been fulfilled. The time before an actual measurement is performed is known as thermalization time and it has to be taken into consideration when planning any Monte-Carlo simulation.
Concerning the degree of independence of each configuration in the Markov chain, it is convenient to define the autocorrelation function Γ(k) as
ΓO(k) =D
O(U(i))O(U(i+k))E
−D
O(U(i))E D
O(U(i+k))E
, (4.10)
which for i τexp is independent of i. The average h·i is taken over an ensemble of chains with independent random numbers and initial states. The autocorrelation function basically estimates the correlation between any two measurements sepa-rated by a distance k in the Markov chain. Using the autocorrelation function one can give an explicit formula for the error σO in the estimatorOˆ
4.2. SIMULATION ALGORITHMS 31
σO2 =D
( ˆO −O¯)2E
= 2τint(O)
n ΓO(0) + O(n−2), (4.11) where τint is the integrated autocorrelation time defined as
τint(O) = 1 2ΓO(0)
k=∞X
k=−∞
ΓO(k). (4.12)
This quantity encodes the autocorrelation in the Markov chain, and its effect in Eq. (4.11) is to reduce the number of independent measurements from n to n/2τint. Notice that for independent configurations τint → 1/2. Unlike τexp, which depends only on the update algorithm, τint depends also on the particular observable which is being considered, but it is never larger than τexp. The formulas in Eqs. (4.11) and (4.12) can be generalized to more complicated functions of several observables f({Oα}), and for the details, the reader is invited to look at a more complete derivation of the error formula as the one in Ref. [66]. Basically, for the estimator F of the function f one has
σF2 = 2τint,F is given in terms of the correlation function
Γαβ(k) =D
Oα(i)−O¯α Oβ(i+k)−O¯β
E , (4.15)
and the derivatives of the function f evaluated at the true value of the observable Oα. In practice, the autocorrelation function cannot be summed up to arbitrarily large values of k, and the derivatives are computed numerically. Throughout this work, we compute the errors using the method described in Ref. [66].
In the following, we briefly describe some of the most frequently used algorithms to generate the transition probability W(U →U0).
4.2.2 The hybrid overrelaxation algorithm
In the case of the pure gauge theory, most algorithms make use of the locality of the action in Eq. (4.3). A full update of the gauge configuration is then a composition of elementary link updates visited in a particular order. One of the most efficient algorithms is a combination of heatbath and overrelaxation updates, which takes the name of hybrid overrelaxation.
32 4.2. SIMULATION ALGORITHMS Heatbath algorithm
In the heatbath algorithm [67], each gauge link Uµ(x) is replaced by Uµ0(x) chosen at random from the entire group with a probability density given by
p(Uµ0(x))dUµ0(x) =dUµ0(x)e−SG(Uµ0(x)) ∝dUµ0(x) exp β
NRe Tr
Uµ0(x)Σµ(x) , (4.16) where SG(Uµ0(x)) is the gauge action depending only on the link Uµ0(x), while the rest remain fixed. The stapleΣµ(x)is made from the sum of products of gauge links adjacent toUµ0(x). A full lattice update consists in a sweep over all gauge links up-dated one at a time. Although the idea is formally simple, the actual implementation is very challenging forN >2. Considering this, when performing a heatbath update inSU(N), one generally uses a method known as the Cabibbo-Marinari method [68], which updates SU(2) subgroups of SU(N) sequentially.
In the case of SU(2) there are very efficient ways to draw the gauge link Uµ0(x) with the right probability distribution [69, 70]. Basically, the probability distribution in Eq. (4.16) can be rewritten as
p(U0)dU0 ∝dU0exp β
2ζTrU0ω
, (4.17)
where ζ = √
detΣ and ω = Σ/ζ, and we have omitted the x and µ dependence of the variables for simplicity2. Considering that A = U0ω ∈ SU(2) can be expressed asA=a01+ia·σ, where σ are the Pauli matrices and a= (a0,a) is a set of 4real numbers, the problem reduces to generating a0 with probability density
p(a0)da0 = 1−a201/2
eαa0da0, (4.18)
and ai uniformly distributed on a two-sphere of radius (1−a20)1/2. The parameter α = β2ζ encodes the dependence on the inverse lattice coupling β and the staple Σ. Once A has been generated, the gauge link U0 can be reconstructed by right multiplication withw†.
As mentioned previously, the case ofSU(N)relies on theSU(2)updates through the Cabibbo-Marinari method. Briefly, one defines a subset F of SU(N) such that A(i,j)∈F is labelled by the two indices i6=j ,(i, j)∈[1, N]2, and
A(i,j)kl =
(δkl for k, l6=i, j
akl otherwise, (4.19)
2Notice thatω defined in this way is inSU(2), and also that detΣ>0.
4.2. SIMULATION ALGORITHMS 33 where the submatrix akl ∈SU(2). In addition, we require that no subset of SU(N), except SU(N) itself, is invariant under left multiplication by F. The minimal sub-set F which satisfies this property has N − 1 elements, however, it is more effi-cient to choose F with N(N −1)/2 elements given by F =n
A(i,j)o
i<j [71]. Then, the update proceeds by choosing A ∈ F with probability density proportional to exp Re TrAUµ(x)Σµ(x)
. Due to the way in whichAhas been defined in Eq. (4.19), one can show that the problem reduces to generating the submatrix akl ∈ SU(2), so one can use the efficient algorithms for the SU(2) gauge theory. The new link is then given by Uµ0(x) =AUµ(x), and the process is repeated for all A∈F.
Overrelaxation algorithm
Unlike the heatbath update, an overrelaxation update does not perform a random walk, and is known to explore the space of field configurations faster [72, 71]. The idea is the following: for each link, suppose there is a way to find the element Uµ0(x) which minimizes the local action SG(Uµ(x)), then, the overrelaxation update consists in choosing the new element as the link which is in the opposite direction of the original Uµ(x) with respect to the minimum [73, 74]
Uµ0(x) =Uµ0(x)Uµ(x)†Uµ0(x). (4.20) For SU(2), this is achieved by picking Uµ0(x) = q
det
Σµ(x)
Σ−µ1(x), which allows to write Eq. (4.20) in a more convenient way as
Uµ0(x) = 2 Σ†µ(x)Uµ(x)†Σ†µ(x) Tr
Σµ(x)Σ†µ(x) . (4.21)
This algorithm can easily be shown to be microcanonical, so it is common to per-form several overrelaxation updates followed by one heatbath update to guarantee ergodicity.
As with the heatbath update, it is convenient to use a Cabibbo-Marinari strategy to updateSU(N)matrices. We must also point out that a full SU(N)overrelaxation update without considering the SU(2) subgroups can be done [75], and promising results have been seen also for reduced models [71].
4.2.3 Critical slowing down
Lattice simulations are naturally limited by the computational resources available, and the situation only gets worse when approaching the continuum limit. At a fixed physical volume, the increase in the number of points in the lattice is proportional to a−4, but this is not the only difficulty when approaching the continuum. The algorithms that we have presented before suffer from a problem known as “critical
34 4.3. THE YANG-MILLS GRADIENT FLOW