Published for SISSA by <£] Springer

Received: May 30, 2011 Accepted: June 23, 2011 Published: July 7, 2011

Lattice QCD without topology barriers

Martin LUscher and Stefan Schaefer

CERN, Physics Department,

1211 Geneva 23, Switzerland

E-mail: Martin.Luescher@cern.ch, Stefan.Schaefer@cern.ch

Abstract: As the continuum limit is approached, lattice QCD simulations tend to get trapped in the topological charge sectors of field space and may consequently give biased results in practice. We propose to bypass this problem by imposing open (Neumann) boundary conditions on the gauge field in the time direction. The topological charge can then flow in and out of the lattice, while many properties of the theory (the hadron spectrum, for example) are not affected. Extensive simulations of the SU(3) gauge theory, using the HMC and the closely related SMD algorithm, confirm the absence of topology barriers if these boundary conditions are chosen. Moreover, the calculated autocorrelation times are found to scale approximately like the square of the inverse lattice spacing, thus supporting the conjecture that the HMC algorithm is in the universality class of the Langevin equation.

Keywords: Lattice QCD, Lattice Gauge Field Theories, Stochastic Processes ArXiv ePrint: 1105.4749

Open Access

doi:10.1007/JHEP07(2011)036

Contents

1 Introduction

2 QCD with open boundary conditions

2.1 Boundary conditions in the continuum theory

2.2 Topology of the classical field space

2.3 Renormalization and stability of the boundary conditions

2.4 Lattice formulation

2.5 Quantum-mechanical representation

2.6 On-shell O(a) improvement

2.7 Other lattice formulations of the theory

3 Dynamical properties of QCD simulations

3. 1 Autocorrelations

3.2 Topology-changing transitions

3.3 Renormalizable algorithms

3.4 Scaling behaviour of the HMC algorithm

3.5 Making QCD simulations safer

8 8 9 9

Œ M •n o

Numerical studies

4.1 Simulation algorithms

4.2 Stochastic equation, parameter scaling and the SMDY algorithm

4.3 Observables

4.4 Lattice and simulation parameters

11 12 12 14

Simulation results

5.1 Scaling properties of the autocorrelation functions

5.2 Autocorrelation times

5.3 Dependence on the time-like extent of the lattice

6 Conclusions

A Notational conventions

B Calculation of integrated autocorrelation times

18 18 19

1 Introduction

One of the central goals in numerical lattice QCD is the computation of the properties of the light mesons and baryons with controlled errors. While the most important systematic errors in these calculations (finite volume and lattice spacing effects) are theoretically well understood, the relevant time scales in QCD simulations remain unpredictable. In practice, the correctness of the simulations within the quoted statistical errors can therefore only be established through empirical tests and thus only to a limited level of confidence.

In order to preserve the translation symmetry, the lattice theory is usually set up with periodic boundary conditions in all space-time directions. A side-effect of this choice of boundary conditions is the emergence of disconnected topological sectors in the continuum limit. On the lattice the sectors are not strictly separated from each other, but the relative weight in the functional integral of the gauge fields "between the sectors" decreases with a high power of the lattice spacing [1]. As a consequence, transitions from one sector to another tend to be suppressed in the simulations and may eventually become so rare that a proper sampling of the sectors would require far longer runs than are practically feasible [2-5].

In this paper we address both issues, the very long autocorrelation times caused by the emergence of the topological sectors and the lack of theoretical control over the simulations. The first of them we propose to avoid by imposing open boundary conditions on the gauge field in the time direction (see section 2). With this choice, the field space in the continuum theory becomes connected and the topological charge can smoothly flow in and out of spacetime through its boundaries. All statistically relevant parts of field space are therefore expected be accessible to the simulation algorithms without having to cross higher and higher topology barriers as the lattice spacing is reduced.

When properly renormalized, some algorithms may even converge to a well-defined stochastic process in the continuum limit (section 3). In asymptotically free theories, such algorithms have a predictable scaling behaviour as a function of the lattice spacing and are thus theoretically controlled to some extent. The HMC algorithm [6] recently turned out to be non-renormalizable in perturbation theory and is therefore not of this kind [7], but the algorithm may conceivably fall in the universality class of the Langevin equation (whose renormalizablity was established long ago [8, 9]). The empirical tests reported in sections 4 and 5 partly serve to verify that the topology barriers are indeed absent if open boundary conditions are imposed and partly to find out whether the HMC algorithm scales like an algorithm that integrates the Langevin equation.

2 QCD with open boundary conditions

Open boundary conditions are easily imposed in QCD and do not give rise to important theoretical complications. While the discussion in this section is more generally valid, the gauge group is taken to be SU(3) from the beginning and we assume there is a multiplet of quarks in the fundamental representation of SU(3). Our notational conventions are summarized in appendix A.

2.1 Boundary conditions in the continuum theory

In the continuum limit, the gauge and quark fields live on a four-dimensional space-time M with Euclidean metric, time extent T and spatial size L x L x L. Time thus runs from 0 to T, while space is taken to be a three-dimensional torus, i.e. all fields are required to satisfy periodic boundary conditions in the space directions. At time 0 and T, the boundary conditions imposed on the gauge potential AM(x) are

Fok(x)lxo=o = Fok(x)Ixo=T = 0 for all k = 1,2,3, (2.1)

where (x) denotes the gauge-field tensor. Note that these conditions preserve the gauge symmetry and therefore do not constrain the gauge degrees of freedom of the field. In perturbation theory, the boundary conditions on the latter instead derive from the gauge-fixing procedure. If the usual Lorentz-covariant gauge is chosen, for example, the time M

and space components of the gauge potential are found to satisfy Dirichlet and Neumann boundary conditions, respectively.1

In the case of the quark and antiquark fields tp(x) and tp(x), we require that

P+№) l,o=o = P-1>(x)\X0=T = 0, P± = 1(1 ±7o), (2.2)

^)P-|.o=o = ^)P+l.o=T = °- (2-3)

I (rxti{t\lv{x)t\lv{x)\ +

where g0 and M0 are the bare coupling and quark mass matrix.

2.2 Topology of the classical field space

Since M is contractible to a three-dimensional torus, all SU(3) principal bundles over M are trivializable. Smooth classical gauge potentials may therefore be assumed to be globally defined differentiable fields. In view of the non-linearity of the boundary conditions (2.1), the classical field space is however not a linear space.

We now show that the field space is connected. Evidently, any given gauge potential AM(x) satisfying A0(x) = d0Ak(x) = 0 at x0 = 0 and x0 = T can be smoothly contracted to zero, without violating the boundary conditions, by multiplication with a scale factor. These fields are therefore continuously connected to the classical vacuum configuration.

1 As is the case with periodic boundary conditions, perturbation theory in finite volume is complicated by the presence of non-trivial gauge-field configurations with vanishing action (the constant Abelian fields). The remarks made here and below on perturbation theory refer to the situation at L = to and finite T, where the minimum of the action is unique up to gauge transformations.

These boundary conditions are familiar from the discussion of the QCD Schrodinger functional [10, 11]. Many of the theoretical results obtained in that context can actually be reused here. In particular, as explained in ref. [12], one is practically forced to choose the boundary conditions (2.2), (2.3) if parity and the time reflection symmetry are to be

preserved. The action of the theory (without gauge fixing terms) is then given as usual by o

S = TT~2 I d4.rtr{F^(.ï)F^(.ï)} + / d4.r Mx) (7»D» + M0) 0(.r), (2.4) 2go J m J M

On the other hand, if one starts from an arbitrary field AM(x) in the classical field space, a smooth curve of gauge transformations As(x), 0 < s < 1, may be defined through the differential equation

(do + sAo(x))As(x)-i = 0, As(x)|x0=o = 1- (2.5)

When applied to the potential AM(x), the transformation generates a curve in field space (parametrized by s) along which the field is continuously deformed to another field at s = 1 with vanishing time component. The transformed field can then be contracted to zero, as explained above, which proves that the field space is connected.

The absence of disconnected topological sectors goes along with the fact that the topological charge

is not quantized. When an instanton on M is contracted to the vacuum configuration, for example, the charge flows away through the boundaries and Q smoothly varies from 1 to

0. It may be worth noting here that the massless Dirac operator has no zero modes in the space of complex quark fields satisfying the boundary conditions (2.2). Moreover, the eigenvalues A of the Hermitian Dirac operator 75 (y^D^ + m) are all in the range |A| > m (see ref. [12], section 2.2, for a proof of these statements). As long as the quark masses are non-negative, the quark determinant has therefore a definite sign and never passes through zero even if some masses vanish.

Q = J d4xq{x), q{x) =-^^e^p<jti{F^(x)Fp<T(x)}, (2.6)

2.3 Renormalization and stability of the boundary conditions

The renormalization of quantum field theories on space-time manifolds with boundaries in general requires the usual (bulk) counterterms to be added to the action as well as further counterterms that are localized at the boundaries [13]. In the present case, the symmetries of the theory, power counting and the fact that t/np vanishes at the boundaries however exclude such boundary counterterms. The renormalization of the theory thus proceeds as in infinite volume by renormalizing the coupling, the quark masses and the fields in the correlation functions considered.

Boundary conditions are subject to renormalization too and sometimes require a fine-tuning of boundary counterterms. Neumann boundary conditions in scalar field theories, for example, are known to be unstable under quantum fluctuations [13]. The situation in QCD is safe from this point of view, because there are no relevant or marginal boundary counterterms with the required symmetries. In particular, the boundary conditions (2.1)-(2.3) are stable under quantum fluctuations (see ref. [12], section 3, for a broader discussion of the subject).

2.4 Lattice formulation

The lattice theory is set up on a hypercubic lattice with spacing a, time-like extent T + a and spatial size L x L x L, where T and L are integer multiples of a. Periodic boundary

conditions are imposed on all fields in the space directions, while time runs from 0 to T inclusive, the terminal time-slices being the boundaries of the lattice.

As usual the gauge and quark fields reside on the links and points of the lattice. In particular, the link variables U(x, € SU(3) live on all links (x, with both endpoints in the specified range of time. The Wilson gauge action is then given by the sum

SG = \^2w{p)tT{l-U{p)} (2.7)

over all oriented plaquettes p on the lattice, U(p) being the ordered product of the link variables around p. Only those plaquettes are included in the sum whose corners are in the time interval [0,T]. The weight w(p) is equal to 1 except for the spatial plaquettes at time 0 and T, which have weight 4.

In the case of the fermion fields, a possible choice of the action is

Sf = a4 + M°) (2-8)

xo=a x

Dw = (Vi + VM) - ±aV*VM (2.9)

denotes the Wilson-Dirac operator and the fields are assumed to satisfy the boundary

conditions (2.2), (2.3). Since the action (2.8) depends on the quark fields at times 0 < x0 < T only, one may then just as well set all components of the fields at time 0 and T to zero. The dynamical components of the quark fields are thus those residing in the interior of the lattice.

The functional integral and the basic correlation functions are now defined in the standard manner. Evidently, only the dynamical components of the fermion fields are integrated over in the functional integral. Note that the quark determinant is a product of real factors, one for each quark flavour, since the Wilson-Dirac operator is Y5-hermitian with the chosen boundary conditions. The established QCD simulation algorithms can therefore be applied straightforwardly.

It may not be completely obvious at this point that the fields satisfy the boundary conditions (2.1)-(2.3) in the continuum limit. As already mentioned in subsection 2.3, these boundary conditions are stable under quantum fluctuations, i.e. it suffices to check that they emerge at tree-level of perturbation theory when the lattice spacing is sent to zero. The explicit expression for the quark propagator obtained in ref. [14] and a similar computation of the gluon propagator in the standard covariant gauge actually show this to be so.

2.5 Quantum-mechanical representation

The formulation of the lattice theory described above admits a quantum mechanical description in terms of a Hilbert space H of physical states and a bounded, positive-definite transfer matrix T [15]. In particular, the partition function of the theory is given by the expectation value

Z = (Q, TT/a Q) (2.10)

of a power of the transfer matrix in a state Q € H that encodes the chosen boundary conditions at time 0 and T.

A relatively direct way of introducing the transfer matrix formalism starts from a representation of the physical states through wave functions depending on a gauge field V(x, k) and the components

X-{x) = P-X(x), x+{x)=x(x)P+ (2.11)

of a quark field on the spatial lattice (see ref. [11], for example). The boundary state Q is represented by the wave function

Q(V, X-,X+) = {det(l + aM0 - ±a2VA*Vfc)} (2.12)

in this language, the covariant derivatives being evaluated in presence of the gauge field V. Note that this expression is manifestly invariant under the gauge symmetry, the lattice symmetries and the (vector-like) flavour transformations. In other words, Q has the quantum numbers of the vacuum state.

Correlation functions of gauge-invariant fields have a quantum mechanical interpretation as well. The two-point function of a local scalar field 0(x), for example, is given by

(<t){x)<t){y)) = -(Q,T(T--To)/a0(.T)T(-T°-№)/a(/)(y)T№/aQ), (2.13)

xo>yo Z

■. I-1

where 0(x) denotes the operator field associated to 0(x) [15]. Since the transfer matrix and

the space of physical states are independent of the boundary conditions at time 0 and T,

the hadron masses and many other physical quantities can, in principle, be extracted from

such correlation functions in much the same way as on lattices with periodic boundary o

conditions in time.

2.6 On-shell O(a) improvement

The O(a) improvement of the lattice theory follows the lines of refs. [16, 17]. There is actually very little difference with respect to the case of the Schrodinger functional discussed in the second of these papers. In particular, all bulk O(a) counterterms and their coefficients are exactly the same as those required for the on-shell improvement of the theory on the infinite lattice.

We wish to emphasize at this point that a further improvement is not needed if one is exclusively interested in the correlation functions of fields localized far away from the boundaries of the lattice, where the effects of the latter are exponentially suppressed. The improvement of correlation functions involving fields close to or at the boundaries however requires the addition of the O(a) boundary counterterms

SSG,b = A(cg - l)£tr{l - U(ps)}, (2.14)

2go Ps

¿SF,b = (cf - l)a3 J2{$(xMx)\xo=a + MxMx)\xo=T_a}. (2.15)

In these equations, ps runs over all space-like oriented plaquettes at the boundaries of the lattice and the coefficients

CX = 1 + cX1)g02 + cX2) g4 + ... (2.16)

must be adjusted so as to cancel the boundary effects of order a (since the boundary conditions are not the same, there is no reason to expect these coefficients to coincide with those needed to improve the Schrodinger functional).

2.7 Other lattice formulations of the theory

Lattice QCD with open boundary conditions can be set up in many different ways. Uni- C-i

versality actually suggests that the details of the lattice theory become irrelevant in the continuum limit if the gauge, space-time and flavour symmetries are respected. Lattice formulations that preserve chiral symmetry away from the boundaries exist as well, but some care is required in this case in order to guarantee the locality of the lattice Dirac operator near the boundaries [12, 18, 19].

3 Dynamical properties of QCD simulations

The interpretation of simulation data requires good control over the simulation dynamics. In this section, the relevant notions are briefly discussed and some specific issues are addressed, which arise when studying the scaling behaviour of QCD simulations.

3.1 Autocorrelations

QCD simulation algorithms produce random sequences of gauge-field configurations recursively, where the next configuration is obtained from the current one according to some transition probability. The simulation algorithms considered in this paper are the HMC algorithm [6] and the closely related SMD (stochastic molecular dynamics, or generalized HMC) algorithm [20-22]. In both cases the simulation time is proportional to the molecular-dynamics time in lattice units and will, for simplicity, be identified with the latter in the following.

Let Oi be a set of real-valued unbiased observables labeled by an index i. Their values Oj(i) measured at simulation time t are statistically correlated to some extent, i.e. the connected parts of the n-point autocorrelation functions

A(ti,... ,tn)i!...in = ((Oi! (tl) ... Oin (tra)» (3.1)

in general do not vanish. In this equation, the bracket ((...»» stands for the average over infinitely many statistically independent parallel simulations, which is the same as the average over time translations if the simulation is ergodic.

The connected parts of the autocorrelation functions tend to fall off exponentially at large separations in simulation time. In particular, the two-point autocorrelation functions

rij(t) = ((Oi(t)Oj(0)»» - ((Oi(t)»»((Oj(0)»» (3.2)

can be shown to have a spectral decomposition of the formo

Tj (t) = J]Re{ cinCjn Ají1}, |A„ | = e-1/Tn, (3.3)

where t0 > t1 > ... are the so-called exponential autocorrelation times of the algorithm. While these are independent of the observables considered, the coefficients cin measure how strongly the observables Oi couple to the eigenmode number n of the transition probability. Note that neither the spectral values An nor the coefficients cin are guaranteed to be real, except in the case of the HMC algorithm and the Langevin limit of the SMD algorithm.

In practice, the integrated autocorrelation times C-i

^ r (t)

rmt(Oi) = ±Aí +Aí^p^Aí), Pi(t) = (3.4)

fc=i iii(0)

of the observables of interest play an important rôle, where At is the separation in simulation time of the observable measurements. The sum in eq. (3.4) amounts to a numerical integration of the normalized autocorrelation function pi(t) using the trapezoidal rule. In particular, in the Langevin limit of the SMD algorithm or if the HMC algorithm is used, the formula

TWO*) = (3-5)

Sri=0(cira )

and thus the bound Tint(Oi) < t0 hold up to integration errors.

3.2 Topology-changing transitions

On lattices with periodic boundary conditions, the probability per unit simulation time for a HMC or an SMD trajectory to pass from one topological sector to another is a rapidly decreasing function of the lattice spacing [2-5]. Such topology-tunneling transitions are non-perturbative lattice artifacts that may informally be described as "an instanton falling through the lattice". The integrated autocorrelation time of the topological charge consequently tends to become very large, sometimes to the extent that the correctness of the simulation is compromised.

With open boundary conditions, the situation is different, because the topological charge can change smoothly along a molecular-dynamics trajectory by flowing in and out of the lattice via its boundaries. A catastrophic slowdown of the algorithms as in the case of periodic boundary conditions is therefore not expected.

3.3 Renormalizable algorithms

The n-point autocorrelation functions of gauge-invariant local fields formally look like the correlation functions in a field theory in five dimensions, where the simulation time is the

2The HMC and the SMD algorithm both evolve the gauge field U(x,^) together with its canonical momentum n(x,^) (cf. subsection 4.1). Equation (3.3) is partly a consequence of detailed balance in phase space and only holds for observables that do not depend on the momentum.

fifth space-time coordinate. When the lattice spacing is taken to zero, the autocorrelation functions may then conceivably have a continuum limit, provided the fields and the parameters (of both the theory and the simulation algorithm) are properly renormalized.

Algorithms that integrate the Langevin equation are known to be renormalizable in this sense to all orders of perturbation theory [8, 9]. An example of an algorithm of this kind is provided by the SMDY algorithm (cf. subsection 4.2). The simulation time has physical dimension [length]2 in this case and must be renormalized according to

t = ZttR/a2 (3.6)

where t is the simulation time in lattice units, Zt(g0) a renormalization constant and tR the renormalized simulation time in some physical units. Further renormalization is not required apart from the usual field and parameter renormalization.

Beyond perturbation theory, the renormalizability of an algorithm (and thus the associated scaling laws) can break down as a result of non-perturbative lattice artifacts. On lattices with periodic boundary conditions, topology-changing transitions have this effect in the case of the SMDY algorithm. However, if open boundary conditions are chosen, there is currently no reason to expect that the renormalizability of the algorithm does not extend to the non-perturbative level.

3.4 Scaling behaviour of the HMC algorithm

Free-field studies of the HMC algorithm suggest that the exponential autocorrelation times scale linearly (like a-1) if the length of the molecular-dynamics trajectories is scaled accordingly [23]. The algorithm however turns out to be non-renormalizable in perturbation theory [7] and its scaling behaviour in the presence of interactions may therefore be completely different.

The empirical studies reported later actually show that the HMC algorithm (on lattices with open boundary conditions) appears to fall into the universality class of the Langevin equation. In particular, the autocorrelation times scale approximately like a-2 rather than linearly. From this point of view, the non-renormalizability of the HMC algorithm in perturbation theory merely reflects the fact that the leading-order theory is in the wrong dynamical universality class and therefore not a suitable starting point for the perturbation expansion.

3.5 Making QCD simulations safer

In practice, numerical simulations should be much longer (by, say, a factor 100 at least) than the longest exponential autocorrelation time t0 , as otherwise a proper sampling of the functional integral is not guaranteed and the simulation may consequently be biased in an unpredictable way. Usually the integrated autocorrelation times of the quantities of interest are monitored, but it should be noted that the correctness of the simulation results (within the estimated statistical errors) cannot be taken for granted if only these autocorrelation times are much smaller than the total simulation time.

Tínt " 80-

Figure 1. Autocorrelation time of the density E at physical time L/2 in the SU(3) gauge theory, plotted as a function of the flow time t (cf. subsect. 4.3). The simulation data (points) were obtained on a lattice of size 324 with spacing a = 0.05 fm and open boundary conditions, using the SMD0.3 algorithm. The line is a fit to the data of the form Tint = c0 — cie-C2t with c0 ~ 94, while the leading exponential autocorrelation time in the even-parity sector is found to be about 100 in these simulations.

Integrated autocorrelation times of both physical and other observables can in fact be very much smaller than t0. In particular, the autocorrelation times of noisy quantities (large Wilson loops, for example) tend to be practically unrelated to the exponential autocorrelation times. To illustrate this point, consider an observable

Oo = Oi + cn,

where c is a constant and n a statistically independent Gaussian noise with mean zero and unit variance. O0 has the same expectation value as O1 and its autocorrelation function is given by

roo(t) = c%o + rn(t). (3.8)

At large c, i.e. when the added noise term is large, the integrated autocorrelation time Tint(Oo) decreases like 1/c2 and can therefore be made arbitrarily small. Nothing is gained in this way, but the example shows that integrated autocorrelation times may not be representative of the true autocorrelations in the simulation.

Exponential autocorrelation times are difficult to determine reliably if very long simulations are impractical. In this case, a pragmatic way to proceed is to look for observables with large integrated autocorrelation times and to take the maximum of the latter as an estimate of t0. The observables that provide the best probes for autocorrelations should be sensitive to the smooth modes of the gauge field, since these tend to be updated least efficiently. Moreover, for the reasons given above, good probes are likely to have small statistical fluctuations. Quantities obtained by integrating the Wilson flow [1], such as the average action density E at positive flow time, satisfy both criteria and are therefore recommended probes (see figure 1).

EC M •n o

4 Numerical studies

In order to verify and complement the theoretical discussion in the previous sections, we performed extensive simulations of the SU(3) gauge theory with open boundary conditions. The algorithms, observables and simulation parameters used in these studies are described in this section.

4.1 Simulation algorithms

Both the HMC and the SMD algorithm operate in phase space, i.e. on the gauge field U(x,^) and its canonical su(3)-valued momentum field The O(a) boundary coun-

terterm (2.14) is not included in the Hamilton function

H(it,U) = ±(IT,IT) + Sg(U) (4.1)

of the system, partly for simplicity and partly because the term is unimportant in the present context.

The HMC algorithm proceeds in cycles, where in each cycle one first chooses the momentum field randomly, with normal distribution, and then evolves the fields according to the molecular-dynamics equations that derive from the Hamilton function (4.1). In our simulations, the equations were integrated from molecular-dynamics time 0 to t using n0 iterations of the 4th-order Omelyan-Mryglod-Folk (OMF) integrator defined through eqs. (63) and (71) in ref. [24]. At the end of the evolution, the fields are submitted to an acceptance-rejection step that corrects for the integration errors. This algorithm has two parameters, t and n0, and requires the derivative of the gauge action to be calculated 5no times per cycle.

In the case of the SMD algorithm, one proceeds in essentially the same way, but the momentum field is only partially refreshed according to

n(x,^) ^ c1n(x,p) + c2v(x,p), (4.2)

ci =e-Y<5r, C2 = (1 - c?)1/2, (4.3)

where u(x, is a randomly chosen momentum field with normal distribution, while y and 5r are parameters of the algorithm. The molecular-dynamics equations are then integrated from 0 to 5r by applying a single iteration of the 4th-order OMF integrator and the fields are finally submitted to an acceptance-rejection step. When rejected, the fields are reset to their values before the integration, except for a change of sign

n(x,^) ^ —n(x,^) (4.4)

of the momentum field (see ref. [25] for a straightforward proof of the correctness of the algorithm). Note that the simulation time t elapsed after n SMD cycles is, by definition, equal to n5r, irrespectively of the rejection rate.

Since the OMF integrator is applied only once, the molecular-dynamics evolution time 5r is usually set to a value much smaller than 1 in order to guarantee a high acceptance

rate Pacc. Otherwise the SMD algorithm is frequently backtracking, on average after every period of time equal to

taCc = 6r facc , (4.5)

1 P acc

and thus tends to become inefficient. With respect to the leapfrog and the 2nd-order OMF integrator, the 4th-order OMF integrator has the advantage that very high acceptance rates can be achieved with a moderate computational effort.

4.2 Stochastic equation, parameter scaling and the SMDY algorithm

In the limit ¿t — 0, the SMD algorithm amounts to solving the stochastic molecular-dynamics equations

dsUs(x,^)= (4.6)

ds = —Ta(dX;MSo)(Us) — + (4.7)

where ns a Gaussian random noise with mean zero and variance o

(n,a(x, v)) = 4^%* ¿(s — r)a-X. (4.8)

In these equations, the evolution time s and the mass are related to the simulation time t and the parameter y through

s = ta, = Y/2a, (4.9)

respectively. Evidently, eqs. (4.6), (4.7) reduce to the standard molecular-dynamics equations if is set to zero (see appendix A for the definition of the derivative of the gauge action).

When the continuum limit is approached, the scaling behaviour of the simulation algorithms depends on how their parameters are scaled. The fact that the evolution time in eqs. (4.6), (4.7) has dimension [length] suggests to scale the HMC trajectory length t proportionally to 1/a [23]. For the same reason, one can argue that should be scaled like a physical mass up to a logarithmically varying renormalization factor perhaps. This choice of the parameter scaling (which, however, leads to non-removable ultra-violet singularities in perturbation theory [7]) will be referred to as free-field scaling.

Alternatively, if 7 is held fixed, and if 5t is such that the continuous-evolution time tacc is on the order of the exponential autocorrelation times (or larger), the SMD algorithm effectively performs a numerical integration of the Langevin equation [7]. For clarity, we use the acronym SMDY for the SMD algorithm with this parameter scaling.

4.3 Observables

As already noted in subsection 3.5, observables based on the Wilson flow probe the slow modes of the gauge field and are therefore well suited for studying autocorrelations in QCD simulations. A review of the Wilson flow is beyond the scope of this paper and we merely write down the differential equation

díVí(x,^) = -ag2Ta(da;MSo)(Ví)Ví(x,^), Vt(x,^)|t=o = U (x,^), (4.10)

L/a /3 a [fm] t0/a2

16 5.96 0.1000(6) 2.698(3)

20 6.09 0.0802(5) 4.203(5)

24 6.21 0.0667(5) 6.086(7)

32 6.42 0.0500(4) 11.045(15)

40 6.59 0.0402(3) 17.49(4)

Table 1. Lattice parameters and reference flow time to/a2 calculated at physical time L/2 on the (L/a)4 lattices.

that generates the flow V(x,^), t > 0, in the space of gauge fields (see refs. [1, 5, 26] for a comprehensive discussion of the flow and some of its surprising properties).

In the course of the simulations, the observables are evaluated at fixed separations M

in simulation time. Starting from the current gauge-field configuration U(x,^), we first integrate the flow equation (4.10) numerically up to some flow time t. The field tensor G^v (x) of the gauge field V(x,^) generated in this way is defined through the clover formula, i.e. through the four plaquette Wilson loops in the v)-plane that start and end at x (at the boundaries xo = 0 and xo = T we set G0k (x) = 0). The primary observables considered are then the time-slice averages

E{xo) = -^3 Çtr{G^(x)G^(x)} (4.11)

of the action density and the time-slice sums

Q{xo) = -32^2 zl el-wpa^T{G^ (x)Gpa (x) } (4.12)

of the topological charge density. Evidently, the autocorrelations of the total charge

Q = aJ2 Q(x o) (4-13)

are studied as well. In all these equations, the dependence on the flow time has been suppressed for simplicity.

At positive flow time t, the expectation values of arbitrary (finite) products of the observables E(xo), Q(xo) and Q do not require renormalization and are expected to have a well-defined limit when the lattice spacing is taken to zero [1, 26]. While these expectation values do not have any obvious interpretation in terms of glueballs or colour flux tubes, for example, they are properties of the continuum theory which reflect the dynamics of the smooth modes of the gauge field (the smoothing radius being roughly equal to V81). In particular, as explained in ref. [1], on lattices with periodic boundary conditions, the topological charge Q (as defined here) converges to an integer-valued observable in the continuum limit, which labels the topological sectors of field space.

Lattice T n0 p 1 acc At Acnfg

164 2.0 6 0.953 6 30697

204 2.5 9 0.975 10 25713

244 3.0 12 0.979 15 25625

324 4.0 20 0.985 24 24041

Table 2. Parameters of the HMC algorithm.

Lattice ÖT tacc At Acnfg

164 0.1410 516(2) 5.92 35093

204 0.1128 748(2) 9.14 20209

244 0.0940 1009(3) 14.5 20521

48 > < 243 0.0818 1205(2) 13.7 70000

80 > < 243 0.0809 964(2) 13.6 40991

324 0.0705 1633(3) 23.7 20473

404 0.0564 2368(5) 38.6 15976

Table 3. Parameters of the SMD0.3 algorithm.

4.4 Lattice and simulation parameters

In table 1 we list the spatial sizes and the inverse gauge couplings @ = 6/g^ of the lattices that we have simulated. The number of lattice points in the time direction (which is equal to T/a + 1) coincides with L/a in most cases, but lattices with larger time extent have been considered too. For the conversion to physical units, we use the Sommer radius r0 = 0.5fm [27] and the results obtained for r0/a by Necco and Sommer [28]. The values of the lattice spacing determined in this way (3rd column of table 1) are such that all lattices have the same physical size L, as is desirable for a scaling study.

As a reference for the Wilson flow time t, we prefer to use the scale t0 determined through the implicit equation [1]

{t2(E(L/2)}}t=to=0.3. (4.14)

At flow time t0, the Wilson flow has a smoothing range approximately equal to r0, i.e. this point in flow time is about where the non-perturbative regime sets in. Since L is quite small in physical units, the values of t0/a2 quoted in table 1 are probably affected by finite-volume effects and they are, in fact, a few percent lower than those previously obtained in ref. [1] at L ~ 2.4fm. In the present context, the effect can however be safely ignored since L is the same on the lattices simulated.

The trajectory length t in the HMC simulations was scaled according to the free-field parameter scaling, and the number n0 of integration steps (each consisting of one iteration of the 4th-order OMF integrator) was then tuned to achieve fairly high acceptance rates Pacc (see table 2). In a second set of simulations, we used the SMDY algorithm. Some experimenting suggests that the autocorrelation times of the observables considered have a flat minimum near 7 = 0.3 and we therefore decided to stick to this value of 7. The

1 p(t) l 1 1 1 1 | 1 1 1 : k <i , | , , E 11 1 1 1 1 1 1 1 1 1 1 1 1 Q2 " 1 1 1 I I 1 1 I I L ' , Q2'-

- ► « - ► A T 0

T + * t * * . ♦ n t

0.1 - ■ 404 : ► 324 - „ 244 - 204 > 1 ^ t 1 1 1 1 1 1 1 1 1 1 1 1 1 + 'V + + + j H 1 \ 1 " 1 1 1 1 1 1 1 1 ll 1 1 1

- 164 1 1 1 1 1 1 1 1 t , l , ,

0.1 0.2 0 0.1 0.2 0

t (a /L)2

Figure 2. Normalized autocorrelation functions of the observables E(L/2), Q(L/2)2 and Q2 at flow time t0, plotted as a function of the simulation time lag t given in units of (L/a)2. The SMD0.3 algorithm was used all cases shown here. For better legibility, the data points obtained on the coarsest lattices (164 and 204) are coloured in grey, while the black points are those from the other lattices (244, 324 and 404).

other parameter of the algorithm, 6r, was adjusted to ensure acceptance over an average simulation time tacc significantly larger than the exponential autocorrelation times (see table 3).

The observables were measured at the separations At in simulation time quoted in tables 2 and 3. On each lattice, a fairly large number Ncnfg of configurations were analyzed, the total length of the simulations thus being equal to NcnfgAt.

5 Simulation results

5.1 Scaling properties of the autocorrelation functions

While the chosen observables do not require renormalization, the flow time at which they are evaluated should be scaled like a physical quantity of dimension [length]2 in the continuum limit. In the following, the flow time is set to the reference time to, the results at other values of the flow time being similar as long as the short-time regime (where lattice effects are large) is avoided.

The SMD0.3 algorithm is renormalizable to all orders of perturbation theory since it effectively integrates the Langevin equation (cf. subsects. 3.3 and 4.2). Moreover, with open boundary conditions, the topology barriers that otherwise slow down the algorithm are absent. It is therefore not unreasonable to expect that the normalized autocorrelation functions of the selected observables converge to universal functions in the continuum limit, provided the simulation time is scaled according to eq. (3.6).

The autocorrelation functions plotted in figure 2 in fact behave as expected if one assumes that the renormalization constant Zt varies only slightly on the lattices considered. Note that all points obtained on a given lattice are statistically correlated. In particular, the seemingly systematic deviation of the measured autocorrelation functions on the 404 lattice

T,nt /z:

1500 0

1500 0

1000 1500

(L /a)2

Figure 3. Integrated autocorrelation times of the observables E(L/2), Q(L/2)2 and Q2 at flow time t0, as obtained on the (L/a)4 lattices using the HMC algorithm (open circles, scale factor Z = 1.32) and the SMD0.3 algorithm (full circles, Z = 1). Many HMC points lie on top of the SMD0.3 points and thus mask the latter. The curves are straight-line fits of the SMD0.3 data.

from those on the 324 and 244 lattices may very well be a statistical fluctuation. Large deviations are however seen in the case of the time-slice and the total topological charge on the coarser lattices, where topology-tunneling transitions are not totally suppressed and thus reduce the autocorrelations. Langevin scaling then sets in as expected once these lattice artifacts become unimportant.

On physically large lattices, the four-point autocorrelation function of the topological charge Q is dominated by its disconnected parts. The normalized two-point autocorrelation function of Q2 is then related to the one of Q through

PQ2 (t) ^ÍPq(Í)}2

Although the simulated lattices are not very large in physical units, we found that eq. (5.1) is accurately satisfied. In particular, the autocorrelation functions of Q and Q2 scale in practically the same way.

5.2 Autocorrelation times

Similarly to the energy spectrum in finite volume, the exponential autocorrelation times depend on the symmetry sector considered. In particular, eq. (5.1) suggests that the longest autocorrelation time in the odd-parity sector is larger, by a factor 2 perhaps, than the one in the even-parity sector. On the basis of the data shown in figure 2, we estimate that the latter is about 1.2 x (r0/a)2 (thus ranging from 30 to 187) in the case of the SMD0.3 algorithm and the (L/a)4 lattices we have simulated.3 As usual, such estimates should be taken with a grain of salt, because the slowest modes in the system may not couple sufficiently strongly to the measured observables for their effects to be seen in the available data.

In accordance with the conventions adopted in section 3, all autocorrelation times are quoted in units of simulation time (i.e. molecular-dynamics time in lattice units).

EC M •n o

Figure 4. Integrated autocorrelation times of E(xq) and Q(xn)2 at flow time to, plotted as a function of the physical time x0 in lattice units. The data were obtained on the 244, 48 x 243 and 80 x 243 lattices using the SMD0 3 algorithm. For better legibility, the data points obtained on the two smaller lattices are coloured in grey.

The integrated autocorrelation times plotted in figure 3 and their errors were calculated following the lines of appendix B. As is evident from the figure, the autocorrelation times all scale linearly in 1/a2 and thus as expected for algorithms that integrate the Langevin equation. From the point of view of the continuum limit, the intercepts at L/a = 0 of the straight lines in figure 3 are O(a2) lattice corrections to the Langevin scaling, while the ratios of their slopes are universal properties of the simulation dynamics.

Figure 3 also shows that the HMC algorithm (with free-field parameter scaling) scales like the SMD0.3 algorithm. The matching of the autocorrelation times requires a renormal-ization of the simulation time by the factor Z ~ 1.32, but in terms of computer time, HMC simulations tend to be faster than SMD0.3 simulations, because a very accurate integration of the molecular-dynamics equations is not needed.

EC M •n o

5.3 Dependence on the time-like extent of the lattice

In practice, the time extent of the simulated lattices will often have to be larger than the one of the (L/a)4 lattices in our scaling studies. Autocorrelations in general depend on the physical situation and thus also on the lattice geometry. For illustration, the autocorrelation times of E and Q2 calculated on three lattices with the same spacing and spatial size, but different time extent T, are plotted in figure 4. Close to boundaries of the lattice, the autocorrelation times shown in these plots are thus practically independent of T, while well inside the lattices they rapidly converge to a constant value when T is increased. The behaviour of the autocorrelation times of these observables discussed in subsection 5.2 is therefore expected to be representative of the situation on larger lattices as well.

The total topological charge Q is a special case, because it can only change (at small lattice spacings) by flowing in and out of the lattice. In the course of a simulation, the measured values of Q fluctuate around the origin with a standard deviation that increases proportionally to VTL3 on large lattices. The charge however flows through the boundaries

with a rate proportional to \fl? only. The simulation time required for a significant change in Q must therefore be expected to grow with T (proportionally to T if Q performs a random walk). On the 244, 48 x 243 and 80 x 243 lattices, we actually find that the autocorrelation times of Q2 (42.1(2.5), 113(6) and 148(10), respectively) grow roughly linearly with T.

We wish to conclude this discussion by emphasizing that the autocorrelation times on lattices of a given physical size are expected to scale linearly in 1/a2. Independently of the chosen geometry, the computational effort for HMC simulations with the standard leapfrog integrator, for example, thus scales approximately like 1/a7.

•n o

6 Conclusions C-I

The theoretical and empirical results presented in this paper show that the topology barriers in the SU(3) gauge theory can be avoided by choosing open boundary conditions in the time direction. Moreover, on lattices with these boundary conditions, the HMC and the SMDy simulation algorithm both appear to fall in the dynamical universality class of the Langevin equation, i.e. simulations based on these algorithms slow down proportionally -J

the square of the lattice spacing when the continuum limit is approached.

In our numerical studies, the autocorrelation times of the topological charge (as well as those of observables unrelated to the latter) went up to values greater than 100 in units of molecular-dynamics time. While such autocorrelations may be affordable in a given case, the experience suggests that there is ample room for algorithmic improvements. A separate treatment of the high-frequency and the smooth modes of the gauge field, for example, might be worth considering at this point.

Open boundary conditions can easily be imposed in QCD with a non-zero number of sea quarks. We do not foresee any technical issues when these boundary conditions are chosen, but an interesting theoretical question is whether the Langevin equation remains renormalizable in the presence of the pseudo-fermion fields that need to be introduced to be able to simulate the theory [29, 30].

All simulations reported in this paper were performed on a dedicated PC cluster at CERN. We are grateful to the CERN management for funding this machine and to the CERN IT Department for technical support.

A Notational conventions

The Lie algebra su(3) of SU(3) may be identified with the linear space of all traceless anti-hermitian 3 x 3 matrices. We choose the generators Ta, a = 1,..., 8, of the Lie algebra to be such that

ti{TaTb} = -\5ab. (A.l)

The general element X of su(3) is then given by X = XaTa with real components Xa (repeated indices are automatically summed over). The Euclidean Dirac matrices , ^ = 0,..., 3, are assumed to be hermitian.

Gauge potentials AM(x) take values in su(3) and are normalized such that the field tensor and the covariant derivatives that appear in the Dirac operator are given by

= - dvAM + Av], (A.2)

D = + (A.3)

On the lattice, the gauge-covariant forward and backward difference operators in presence of a lattice gauge field U(x,/) act on the quark field ^(x) according to

v„0(x) = - {U(x, n)ip(x + afi) - tp(x)} , (A.4)

Vft0(x) = - {ip(x) — U(x — a,/i,/x) lip(x — a,/t)} , (A.5)

where a denotes the lattice spacing and / the unit vector in direction /.

The scalar product of any two vector fields w(x,/) and u(x,/) with values in su(3) is normalized such that

(w, u) = -2a4 Y^ tr{w(x, /)u(x, /)}. (A.6)

If F(U) is a differentiable function of the gauge field, its derivative with respect to the link variable U(x, /) in the direction of the generator Ta is defined by

di4lHU) = or^HUt)

TTf s )etTa U (X,/) if (y,v) = ft-S

, Ui(y,V ) = )rrs , ■ (A.7)

t=o IU (y,v) otherwise.

In particular, in the case of a scalar function F(U), the combination TadX;jUF(U) is a vector field with values in su(3) that transforms under the adjoint representation of the gauge group.

B Calculation of integrated autocorrelation times

The integrated autocorrelation times of the selected observables O were obtained as usual from the empirical estimates

A i ttot — t

rnit) = --— V (Oi(s) - Oi) (Ot(s + t)~ Oi) (B.l)

ttot - * s=At

of the autocorrelation functions ^(t), where ttot = NcnfgAt denotes the total simulation time of the run and Oi the average of the measured values of Oi. In all cases, the autocorrelation functions are found to decay exponentially at large time separations with remarkably consistent values of the exponential autocorrelation times. The estimate

Tm\{Oi) \At + At Pi(kAt), pt(t) = (B.2)

therefore rapidly approaches a constant value when the "summation window" W = kmaxAt is sufficiently large.

On the (L/a)4 lattices considered, the summation window for even-parity observables was set to

, , x2 f 6.0 (HMC runs), , x

W = (ro/a)2 x 60 ( )' (B.3)

[4.5 (SMD0.3 runs).

Given the measured exponential autocorrelation times (subsection 5.2), the systematic error that derives from the truncation of the sum (B.2) is estimated to be at most 3% with this choice. The statistical errors of the autocorrelation functions and the integrated autocorrelation times were determined using the Madras-Sokal approximation [31] (see ref. [32], appendix E, for a detailed description of the procedure).

References O

Open Access. This article is distributed under the terms of the Creative Commons Attribution Noncommercial License which permits any noncommercial use, distribution, and reproduction in any medium, provided the original author(s) and source are credited.

[1] [2]

[8] [9;

[11 [12 [13

M. Luscher, Properties and uses of the Wilson flow in lattice QCD, JHEP 08 (2010) 071 [arXiv: 1006.4518] [SPIRES].

L. Del Debbio, H. Panagopoulos and E. Vicari, d-dependence of SU(N) gauge theories, JHEP 08 (2002) 044 [hep-th/0204125] [SPIRES].

S. Schaefer, R. Sommer and F. Virotta, Investigating the critical slowing down of QCD .simulations, PoS LAT2009 (2009) 032 [arXiv:0910.1465] [SPIRES].

S. Schaefer, R. Sommer and F. Virotta, Critical slowing down and error analysis in lattice QCD simulations, Nucl. Phys. B 845 (2011) 93 [arXiv:1009.5228] [SPIRES].

M. Luscher, Topology, the Wilson flow and the HMC algorithm, PoS Lattice 2010 (2010) 015 [arXiv: 1009.5877] [SPIRES].

S. Duane, A.D. Kennedy, B.J. Pendleton and D. Roweth, Hybrid Monte Carlo, (T)

Phys. Lett. B 195 (1987) 216 [SPIRES].

M. Luscher and S. Schaefer, Non-renormalizability of the HMC algorithm, JHEP 04 (2011) 104 [arXiv: 1103.1810] [SPIRES].

J. Zinn-Justin, Renormalization and stochastic quantization, Nucl. Phys. B 275 (1986) 135 [SPIRES].

J. Zinn-Justin and D. Zwanziger, Ward identities for the stochastic quantization of gauge fields, Nucl. Phys. B 295 (1988) 297 [SPIRES].

M. Luscher, R. Narayanan, P. Weisz and U. Wolff, The Schrodinger functional: A renormalizable probe for non-Abelian gauge theories, Nucl. Phys. B 384 (1992) 168 [hep-lat/9207009] [SPIRES].

S. Sint, On the Schrddinger functional in QCD, Nucl. Phys. B 421 (1994) 135 [hep-lat/9312079] [SPIRES].

M. Luscher, The Schrodinger functional in lattice QCD with exact chiral symmetry, JHEP 05 (2006) 042 [hep-lat/0603029] [SPIRES].

K. Symanzik, Schrodinger representation and Casimir effect in renormalizable (quantum field theory, Nucl. Phys. B 190 (1981) 1 [SPIRES].

[14] M. Lüscher and P. Weisz, O(a) improvement of the axial current in lattice QCD to one-loop order of perturbation theory, Nucl. Phys. B 479 (1996) 429 [hep-lat/9606016] [SPIRES].

[15] M. Lüscher, Construction of a selfadjoint, strictly positive transfer matrix for Euclidean lattice gauge theories, Commun. Math. Phys. 54 (1977) 283 [SPIRES].

[16] B. Sheikholeslami and R. Wohlert, Improved continuum limit lattice action for QCD with Wilson fermions, Nucl. Phys. B 259 (1985) 572 [SPIRES].

[17] M. Lüscher, S. Sint, R. Sommer and P. Weisz, Chiral symmetry and O(a) improvement in lattice QCD, Nucl. Phys. B 478 (1996) 365 [hep-lat/9605038] [SPIRES].

[18] Y. Taniguchi, Schrodinger functional formalism with Ginsparg-Wilson fermion, JHEP 12 (2005) 037 [hep-lat/0412024] [SPIRES].

[19] Y. Taniguchi, Schrodinger functional formalism with domain-wall fermion, JHEP 10 (2006) 027 [hep-lat/0604002] [SPIRES].

[20] A.M. Horowitz, Stochastic quantization in phase space, Phys. Lett. B 156 (1985) 89 [SPIRES].

[21] A.M. Horowitz, The second order Langevin equation and numerical simulations, Nucl. Phys. B 280 (1987) 510 [SPIRES].

[22] A.M. Horowitz, A generalized guided Monte Carlo algorithm, Phys. Lett. B 268 (1991) 247 [SPIRES].

[23] A.D. Kennedy and B. Pendleton, Cost of the generalised Hybrid Monte Carlo algorithm for free field theory, Nucl. Phys. B 607 (2001) 456 [hep-lat/0008020] [SPIRES].

[24] I.P. Omelyan, I.M. Mryglod, R. Folk, Symplectic analytically integrable decomposition algorithms: classification, derivation, and application to molecular dynamics, (quantum and celestial mechanics simulations, Comp. Phys. Commun. 151 (2003) 272.

[25] K. Jansen and C. Liu, Kramers equation algorithm for simulations of QCD with two flavors O of Wilson fermions and gauge group SU(2), Nucl. Phys. B 453 (1995) 375

[hep-lat/9506020] [SPIRES].

[26] M. Lüscher and P. Weisz, Perturbative analysis of the gradient flow in non-abelian gauge theories, JHEP 02 (2011) 051 [arXiv:1101.0963] [SPIRES].

[27] R. Sommer, A New way to set the energy scale in lattice gauge theories and its applications to the static force and as in SU(2) Yang-Mills theory, Nucl. Phys. B 411 (1994) 839 [hep-lat/9310022] [SPIRES].

[28] S. Necco and R. Sommer, The Nf = 0 heavy quark potential from short to intermediate distances, Nucl. Phys. B 622 (2002) 328 [hep-lat/0108008] [SPIRES].

[29] A. Ukawa and M. Fukugita, Langevin simulation including dynamical quark loops, Phys. Rev. Lett. 55 (1985) 1854 [SPIRES].

[30] G.G. Batrouni et al., Langevin simulations of lattice field theories, Phys. Rev. D 32 (1985) 2736 [SPIRES].

[31] N. Madras and A.D. Sokal, The pivot algorithm: a highly efficient Monte Carlo method for selfavoiding walk, J. Statist. Phys. 50 (1988) 109 [SPIRES].

[32] M. Lüscher, Schwarz-preconditioned HMC algorithm for two-flavour lattice QCD, Comput. Phys. Commun. 165 (2005) 199 [hep-lat/0409106] [SPIRES].