Published for SISSA by Springer

Received: March 21, 2011 Accepted: April 12, 2011 Published: April 26, 2011

Non-renormalizability of the HMC algorithm

Martin Lüscher and Stefan Schaefer

CERN, Physics Department, 1211 Geneva 23, Switzerland

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

Abstract: In lattice field theory, renormalizable simulation algorithms are attractive, because their scaling behaviour as a function of the lattice spacing is predictable. Algorithms implementing the Langevin equation, for example, are known to be renormalizable if the simulated theory is. In this paper we show that the situation is different in the case of the molecular-dynamics evolution on which the HMC algorithm is based. More precisely, studying the theory, we find that the hyperbolic character of the molecular-dynamics equations leads to non-local (and thus non-removable) ultraviolet singularities already at one-loop order of perturbation theory.

Keywords: Lattice QCD, Lattice Quantum Field Theory, Renormalization Regulariza-tion and Renormalons

ArXiv ePrint: 1103.1810

Open Access

doi:10.1007/JHEP04(2011)104

Contents

1 Introduction

2 Stochastic molecular dynamics

2.1 Evolution equations

2.2 Solution of eq. (2.5) to leading order in g0

2.3 Iterative solution of the evolution equation

3 Autocorrelation functions

3.1 Feynman rules

3.2 Computation of the two-point function to one-loop order

3.3 Computation of the diagrams 2 and 3

4 Relation to the ordinary correlation functions

4.1 Two-point function

4.2 Four-point function at leading order

4.3 Four-point function at one-loop order

5 Non-renormalizability of the stochastic molecular dynamics

5.1 Parameter renormalization

5.2 Non-renormalizability of the four-point function

5.3 Implications for the HMC algorithm

5.4 The Langevin limit

6 Concluding remarks

A Properties of the Green function K(t,x)

A.1 Definition

A.2 Time-momentum representation

A.3 Explicit expression for K(t, x) in four dimensions

B Proof of eq. (3.18)

7 •n

12 ■—■

1 Introduction

Numerical simulations in lattice field theory are based on stochastic processes that produce random sequences of representative field configurations. It is often useful to interpret the simulation time in these calculations as a further space-time coordinate. The n-point autocorrelation functions of the local fields then formally look like the correlation functions in a field theory with an extra dimension and they are, in fact, sometimes representable in this

way. Depending on the simulation algorithm, and if the simulated theory is renormaliz-able, the autocorrelation functions may conceivably be renormalizable as well. The scaling properties of such algorithms (which, for brevity, will be referred to as renormalizable) are encoded in the continuum theory and thus become predictable to some extent.

In the pure SU(N) gauge theory, for example, simulation algorithms that integrate the Langevin equation are known to be renormalizable [1, 2]. The integrated autocorrelation times Tint of physical observables have dimension [length]2 in this case. Moreover, the standard renormalization group analysis and a one-loop calculation in perturbation theory [3-5] imply that they scale according to [6]

Tint = Cg09/11 {1 + o(#2)} r2 (1.1)

at small lattice spacings a, where C is an observable-dependent constant, g0 the bare gauge coupling and r0 the Sommer radius [11]. In lattice units, the autocorrelation times thus increase like 1/a2 as a — 0 up to a logarithmically decreasing factor.1

Most simulations of lattice QCD performed today are based on some variant of the HMC algorithm [12]. The form of the underlying molecular-dynamics equations and free-field studies [13] suggest that the simulation time has physical dimension [length] in this case and that the autocorrelation times consequently scale essentially like 1/a. As far we know, the renormalizability of the algorithm has however never been studied and its scaling properties in presence of interactions thus remain unknown.

In this paper, the issue is addressed in the framework of perturbation theory. For simplicity the 4>4 theory is considered, but our main result (the non-renormalizability of the molecular-dynamics equations) no doubt extends to most theories of interest. A slightly generalized version of the HMC algorithm is studied, which was introduced many years ago by Horowitz [14-16] (see sections 2 and 3). The non-renormalizability of the associated stochastic equation is then established by showing that the four-point autocorrelation function of the fundamental field has a non-removable ultraviolet singularity at second order in the coupling (sections 4 and 5).

•n o

2 Stochastic molecular dynamics

In order to simplify the discussion as much as possible, we consider the 04 theory with a single scalar field 0 and dimensional instead of a lattice regularization. The action of the field in D = 4 — 2e Euclidean dimensions is given by

S= f dDx ^d.md.m + ¿m&ix)2 + §</>(*)4} , (2.1)

where m0 denotes the bare mass parameter and g0 the bare coupling constant. The generalized HMC algorithm [14-16] integrates a stochastic version of the molecular-dynamics

1 Equation (1.1) is expected to hold on the infinite lattice and on finite lattices with open boundary

conditions [10]. On lattices with periodic boundary conditions, topology-changing tunneling transitions can

give rise to very large autocorrelation times [7-10]. Such transitions are lattice artifacts and are therefore

not included in the analysis that leads to eq. (1.1).

equations that derive from the action (2.1). In the following subsections, we briefly discuss these equations and solve them in powers of the coupling g0.

2.1 Evolution equations

As usual the molecular dynamics evolves the field $(t, x) together with its momentum n(t, x) as a function of a fictitious time t. The stochastic evolution equations [14-16]

dt $ = n, (2.2)

dtn = -—- 2/X07T + r?

= - m2)</> - g</>3 - 2Atovr + rj, (2.3)

involve another mass parameter, > 0, and a Gaussian noise n(t, x) with vanishing expectation value and variance

(n(t, x)n(s, y)) = 4^cS(t - s)S(x - y). (2.4)

Evidently, the ordinary molecular dynamics is recovered in the limit ^ 0. Moreover, in the second-order form,

9^ + 2/xo dt<l> = -^+v, (2-5)

and after substituting t ^ 2^0t, the evolution equations are seen to coincide with the Langevin equation up to a term that goes to zero at large

Since its introduction by Horowitz [14-16], the generalized HMC algorithm has been occasionally studied in the literature, where it is referred to as the Kramers equation or the L2MC algorithm (see refs. [13, 17, 18], for example). In practice, one starts from the firstorder equations (2.2), (2.3) and implements the algorithm using symplectic integrators and o acceptance-rejection steps. For the theoretical analysis in this paper, we however prefer to proceed with the second-order equation (2.5).

2.2 Solution of eq. (2.5) to leading order in g0

The leading-order equation

= n, (2.6)

D = + 2^odt - d^ + m0, (2.7)

coincides with the Klein-Gordon equation in D + 1 dimensions except for the term proportional to which tends to damp the time evolution of the field. At large ^0 and after a rescaling of t, the equation actually turns into the heat equation. The Green function

K(t, x) = / e-i"t+ipxK(u,p), (2.8)

K(w,p) = (-w2 - 2^0w + p2 + ml)-1, (2.9)

of the differential operator D is discussed in some detail in appendix A. Here and below, the notational convention

/.-/£■ H

is used. It is then straightforward to show that the solution of the wave equation (2.6) at time t > to with prescribed initial data at time to is given by

0o(t,x) = [ ds f dDy K(t - s,x - y)n(s,y) (2.11)

J to J

+ J dDy {K(t - to, x - y)(dt0o)(to, y) + (dt + 2^o)K(t - to, x - y)0o(to, y)} .

Note that the dependence on the initial data dies away exponentially with increasing time (see appendix A). The stochastic molecular-dynamics evolution thus thermalizes and eventually loses all memory of the initial values of the field.

In the following, we shall only be interested in the behaviour of the autocorrelation functions after thermalization. We therefore move the thermalization phase to time to = -to and are then left with the solution

0o(t,x)= f ds f dDy K(t - s,x - y)n(s,y) (2.12)

of eq. (2.6) that describes the stationary situation.

2.3 Iterative solution of the evolution equation

Equation (2.5) may be written in the form

V<f> = rl-^<t>3 (2.13)

or, equivalently, as an integral equation

<f>(t,x)=Mt,x)-^ f ds JdDyK(t-s,x-y)<f>(s,y)3. (2.14)

Iteration of the latter then yields the solution 0(t, x) in powers of go.

Each term in this expansion may be represented by a tree diagram with directed lines, four-point and one-point vertices (see figure 1). In frequency-momentum space,

= J dt dDx e^t—ipx0(t,x), (2.15)

the lines represent the Green function

-- = K(u,p), (2.16)

while the one-point vertices

w, p -> -h8> = r)(u,p) (2.17)

Figure 1. In perturbation theory, the solution of the integral equation (2.14) is given by a series of directed tree diagrams. The diagrams up to second order in g0 are shown in this figure. All diagrams have a single external line (labeled by a little square) with ingoing frequency-momentum (u,p). The arrows on the internal lines all point in the direction towards the external line.

stand for the insertion of the noise field. As in ordinary Feynman diagrams, there is a frequency-momentum conservation ¿-function

(2n)D+15(u 1 + W2 + W3 + U4)S(pi + P2 + P3 + P4) (2.18)

associated to each vertex

1 = -go (2.19)

with ingoing frequency-momenta (w1,p1),..., (u4,p4). The lines in these diagrams are directed, because K(u,p) is not invariant under a change of sign of u. In position space, the arrows point in the direction of increasing simulation time.

3 Autocorrelation functions

The n-point autocorrelation functions of the field <(t, x) are usually defined by taking the time average of the product <(t1, x1)... <(tn, xn) at fixed time lags U — tj. In the present setup, the translation symmetry in time allows the time average to be replaced by the average (...) over the noise field n(s, y). We are thus led to consider the correlation functions

An(U1,P1; ... ; Un,Pn) = (<(U1,P1) ...<(un,Pn)) (3.1)

in frequency-momentum space, which may be computed in perturbation theory by expanding the fields <(uk,Pk) in powers of the coupling g0, following the lines of the previous section, and by contracting the noise fields using Wick's rule. As a result one obtains a sum of Feynman diagrams for the autocorrelation functions similar to the ones for the ordinary (field-theoretical) correlation functions.

3.1 Feynman rules

The one-point vertices in the tree diagrams that represent the terms in the expansion of <(uk,Pk) are connected to the rest of the tree through a directed line. When the noise fields at any two such vertices are contracted, an undirected line

- = G(u,p) (3.2)

Figure 2. One-loop vertex diagrams contributing to the two- and four-point correlation functions of the basic field. Up to permutations of the external lines, these are all one-loop vertex diagrams with less than six legs.

is obtained, where

= 4poK(u,p)K(-u,p). (3.3)

In the case of the two-point autocorrelation function, for example, the contraction of the lowest-order diagram

w, p v,q

leads to

(t(u,p)t(v, q)) = (2n)D+15(u + V)6(p + q)G(u,p) + O(go). (3.5)

The Feynman diagrams for the autocorrelation functions thus involve directed lines (2.16), undirected lines (3.2) and the vertices (2.19). Given these lines and vertices, the Feynman rules are the usual ones except for the following special features:

(a) Both kinds of lines (directed and undirected) can connect to the vertices. The value of the latter is the same in all cases.

(b) Each vertex has exactly one outward-directed line attached to it. This line may be an internal or an external line.

(c) There are no diagrams with loops of directed lines.

(d) External lines may be outward-directed lines or undirected lines. There are no inward-directed external lines.

The structure of the Feynman diagrams is otherwise the same as in any field theory. In particular, the diagrams can be decomposed into one-particle irreducible parts and the lines connecting them. Note that the external legs of the irreducible parts may be directed or undirected (see figure 2). The two types of legs must be distinguished and there are thus many more irreducible diagrams than in the case of the ordinary correlation functions.

3.2 Computation of the two-point function to one-loop order

The decomposition of the two-point autocorrelation function into one-particle irreducible parts reads

A2(w,p; v, q) = (2n)D+15(u + v )S(p + q){G(u,p) (3.6)

+K (u,p)E(u,p)G(u,p) + G(u,p)E(-u,p)K (-w,p) + ■ ■ ■ },

where the self-energy £ is given to one-loop order by the diagram 1 in figure 2. A short calculation, using eq. (A.7), then shows that

S(w,p) = -f/i+0(9§), (3.7)

»1)^-^+0(1). (3.8)

To this order, the self-energy thus coincides with the familiar tadpole diagram that contributes to the ordinary two-point correlation function. In particular, the poles at e = 0 will later be seen to cancel when the mass parameter m0 is renormalized.

3.3 Computation of the diagrams 2 and 3

Apart from diagram 1 (which may be inserted in the external lines), the diagrams 2 and 3

in figure 2 are the only one-particle irreducible diagrams that contribute to the four-point autocorrelation function at one-loop order. Up to a factor and the statistical factor, they are given by

I2(u,p) = I K(v + w,q + p)G(v,q), (3.9)

Is(w,p) = / G(v + w,q + p)G(v,q), (3.10)

where (w,p) is the external frequency-momentum that flows through the diagrams from left to right.

The integrals I2 and I3 can be transformed to a useful alternative form by inserting the time-momentum representation of the propagators (see appendix A). Taking the support properties of the Green function into account, one first notes that

I2(u,p) = --j dte^tdt{G(t,q+p)G(t,q)}. (3.11)

Partial integration then yields

I2(u,p) = ^Jo(p) + Y'h(u,p), (3.12)

Jo(p) = / {((q + p)2 + m2)(q2 + m2)} 1, (3.13)

Ji(w,p) = / dteiWtG(i,q + p)G(t, q). (3.14)

The integral I3 is similarly given by

Is(w,p) = Ji(w,p) + Ji( w,p) (3.15)

so that it suffices to work out the integrals J0 and J1.

The integral Jq coincides with the familiar one-loop diagram contributing to the ordinary four-point correlation function in the theory. It is logarithmically divergent and thus has a pole,

MP) =oI^ + o(i), (3.16)

at e = 0. In four dimensions, the integral J\ has negative dimension of mass and may therefore be expected to be absolutely convergent. The following explicit calculation however shows that this is not so.

The integral over t in eq. (3.14) can be performed analytically starting from the time-momentum representation (A.7) of the field propagator. After some algebra and the substitution q —> q — 7jp, the integration leads to the expression

Ji(u,p) = (2^o - iu) J j2(V + \'P2 + mo) + (2Ato - iu)( 4/x0 -x j 4(qp)2 + (2^q - iw)2

4( q2 + jP2 + rn'o j ~ - M

Q + + mo)(((?-^) + ™o)j • (3-17)

The integrand in this formula is a singularity-free function of w,p and q, but at large q the integral is logarithmically divergent in four dimensions. A somewhat lengthy calculation then shows that

-X s -1

=o ^ (2,0 - M ll + jl + \ + 0(1) (3.18)

k V V J)

(see appendix B; the branch of the square root to be taken is the principal one).

4 Relation to the ordinary correlation functions

Since the stochastic molecular dynamics simulates the field theory with action (2.1), the equal-time autocorrelation functions

Cn(pi,---,pn) = An(wi,pi;...; w„,p„) (4.1)

must coincide with the ordinary correlation functions of the fundamental field in momentum space [14-16]. In this section, we show that the two- and the four-point autocorrelation functions do have this property at one-loop order of perturbation theory. Partly the calculation serves as a consistency check, but some of the intermediate results will be helpful in section 5 as well, where we discuss the non-renormalizability of the stochastic molecular dynamics.

4.1 Two-point function

Recalling the results obtained in subsection 3.2, the two-point autocorrelation function is given by

A2(w,p; v, q) = (2n)D+1^(w + v)S(p + q)

x |g(W>P) + |I^G^p) + 0((/o) j • (4.2)

Using the time-momentum representation (A.7) of the field propagator, the integrals over the frequencies are easily worked out and one recovers the familiar expression

C2(p, q) = (2vt)d5(P + q) {(p2 + m2)"1 - |h(p2 + m2)"2 + O(g2)} (4.3)

for the two-point correlation function in the theory.

À(4]{ui,pi-,... ;w4,p4) i (4.8)

torjLm

4.2 Four-point function at leading order

The leading-order contribution

(wi ... ; W4,p4)=\/ + Y + + y (4.4)

to the four-point autocorrelation function is a sum of four diagrams. If one substitutes o

2n ¿(wi + ■ ■ ■ + W4) = dt e-ii(w1+-+W4) (4.5)

for the frequency-conservation ¿-function, the integrals over the frequencies that lead from the autocorrelation to the ordinary correlation functions factorize and after a few further steps one obtains

C40)(pi,... ,P4) = (2n)D¿(pi + ■ ■ ■ + P4)g0 r° didfc{G(i,pi).. .(?(t,P4)} (4.6)

for the leading-order four-point correlation function. Use has here been made of the identity (A.6) and of the fact that the Green function K(t,p) vanishes at negative times t. Performing the time integration in eq. (4.6), the correlation function

C40)(pi,... ,P4) = (2n)D¿(pi + ■ ■ ■ + p4)(-g0)(p? + m0)-i... (p4 + m0)-i (4.7)

is then seen to coincide with the expected expression. 4.3 Four-point function at one-loop order

The second-order contribution Ai4i)(wi,pi;... ; w4,p4) to the four-point autocorrelation function is a sum of terms proportional to the integrals Ii, J0 and Ji. There are 28 diagrams with an insertion of diagram 1 in one of the external lines. The sum of all these contributions to the four-point function is

As discussed in subsection 3.3, diagram 2 is a linear combination of the integrals J0 and J1, while diagram 3 is expressed through J\ alone. Collecting all terms proportional to J0, their sum is found to be given by

Â(l\ui,pi-,... ;w4,î94) = -^-{M'Pi +P2) + <MPi +P-s) + <Mpi + p4)}

xA40)(wi,pi;... ; u4,p4). (4.9)

Both expressions (4.8) and (4.9) are easily integrated over the frequencies, because the only factor that depends on w1,... ,u4 is the tree-level four-point function. One then discovers that the expected result for C^ (p1,..., p4) is obtained already from these two contributions to the four-point function.

In order to show that the remaining terms (i.e. those proportional to J1) vanish when integrated over the frequencies, first consider the channel where the frequency-momentum M

combination

(w,p) = (wi + W2,pi + P2) (4.10)

flows through the diagrams 2 and 3. Dropping the terms proportional to J0 and using the identity (A.5), one obtains

>0<=- !•«-> x{X->0<-><} (4.11) g

for the contribution of diagram 2 in this channel. The contribution of the other diagram is similarly given by

XX = -f [JiKiO + Ji(-w,p)l X j X }, (4.12)

where the prefactor includes the symmetry factor of the diagram. In total there are six further diagrams that differ from the diagrams 2 and 3 by a different distribution of the arrows to the external lines, each of them being given by the corresponding expression (4.11) or (4.12) with the proper assignment of arrows. The sum of all these contributions to the four-point function is then equal to

- g0 Ji(u,p) X { X + X + X + X } + {Wfc ^ -Wfc} (4.13)

i.e. all terms where the arrows are both ingoing or both outgoing cancel in the sum.

It is not difficult to show that the terms in eq. (4.13) vanish when integrated over the frequencies. The integral over the first term, for example, is proportional to

/ Ô(Ui +-----h U4)Ji(u,p)K(Ui,Pi)G(U2,P2)G(U3,P3)K(-U4,P4). (4.14)

Eliminating w4, the integral becomes

/ Ji(ui + U2,p)K(ui,pi)G(u2,P2)G(u3,p3)K(ui + U2 + ^3,^4), (4.15)

and if one first integrates over w1, the term is seen to vanish, because the integrand has no singularities in the half-plane Im w1 > 0 and falls off at least like w-4 at large w1. Exactly the same happens for all other terms in eq. (4.13) and also those in the other two frequency-momentum channels.

go = M^g{l + ^L- + 0(g2)}, (5.1)

5 Non-renormalizability of the stochastic molecular dynamics

We now address the question whether the ultraviolet singularities of the autocorrelation functions can be canceled by the addition of local counterterms to the evolution equation (2.5).

5.1 Parameter renormalization

Evidently, the list of counterterms must include those corresponding to the usual parameter M

and field renormalization that is required for the renormalization of the ordinary correlation functions. In the minimal subtraction scheme, the bare coupling and mass are related to the renormalized parameters g and m through

where M denotes the normalization mass. To one-loop order, the fundamental field does not need to be renormalized in this theory.

Recalling eqs. (4.2), (4.8) and (4.9), it is then immediately clear that the parameter renormalization cancels the singularities of the two- and four-point autocorrelation functions which derive from the poles (3.8) and (3.16) of the integrals I1 and J0.

5.2 Non-renormalizability of the four-point function

The four-point autocorrelation function has further singularities proportional to the divergent part (3.18) of the integral J1. As explained in section 4, the ordinary four-point correlation function does not receive any contributions from this integral (and is therefore finite after the parameter renormalization), but the terms proportional to J1 do contribute to the autocorrelation function at non-zero time separations.

The residue of the pole in eq. (3.18) is the Fourier transform of a distribution

e-2^o t

32^0mt2 ~ X'2) {5-3)

supported on the light cone t = |x|. Both diagrams 2 and 3 thus have a non-local singularity that cannot be canceled by including local counterterms in the stochastic molecular dynamics. The latter is therefore not renormalizable.

The presence of the singularity (5.3) can be understood by noting that the integrand of the integral

J1 (w,p) = J dty dDxeiwt-ipxG(t,x)2 (5.4)

random initial momentum

Figure 3. The HMC algorithm moves the fundamental field < through field space along a piecewise smooth curve. In the smooth segments of the curve, the field is evolved from time t = 0 to some time t = t according to the molecular-dynamics equations, starting from the current field < and Gaussian random values for its momentum n.

5.3 Implications for the HMC algorithm

In practice, the HMC algorithm involves a numerical integration of the (ordinary) molecular-dynamics equations and acceptance-rejection steps to correct for the integration errors. For simplicity the integration is assumed to be exact in this section. No acceptance-rejection steps are then required and whether one uses the first- or the second-order form of the molecular-dynamics equations makes no difference.

The molecular-dynamics trajectories generated by the algorithm are smooth segments of a continuous curve in field space (see figure 3). Along the trajectories, the n-point autocorrelation functions in the time-momentum representation,

An(ti,pi;...; tn ,Pn) = (<(ti,pi) ...<j)(tn,Pn)), 0 < tk < t, (5.5)

may be defined, where the bracket (...) stands for the average over all trajectories in an infinitely long simulation. The autocorrelation functions (5.5) only describe the dynamical properties of the algorithm in the specified range of times, but the discussion in the following paragraphs shows that already these correlation functions are not renormalizable.

The average over trajectories in eq. (5.5) amounts to taking the average over the initial values of the field < and its momentum n = dt<. Since these are distributed according to the equilibrium distribution (a Gaussian in the case of the momentum), the average coincides with the ordinary expectation value. In perturbation theory, the correlation

has a non-integrable singularity in D = 4 dimensions proportional to (t — |x|)-1 (see subsection A.3). Such light-cone singularities are a characteristic feature of Green functions of hyperbolic wave equations and the non-renormalizability of the stochastic molecular dynamics is thus seen to be related to its hyperbolic nature. O

In ordinary field theory, one-loop integrals do not have non-local ultraviolet singularities, because they can be Wick rotated to Euclidean space where the propagators are singular at coinciding points only. The spectral condition and the locality of the theory guarantee that no singularities stand in the way of the Wick rotation [19]. In the case of the diagrams 2 and 3, however, the integrands have poles in all quadrants of the complex frequency plane and the integrals (3.9) and (3.10) consequently cannot be Wick rotated without generating additional terms.

functions can therefore be calculated by solving the (non-stochastic) molecular-dynamics equations in the range 0 < t < t with prescribed initial data at t = 0 and by computing the expectation value of the product </»(t1,p1)... 0(tn,pn) using the standard Feynman rules for the correlation functions of the initial data.

In the case of the stochastic molecular dynamics, the computation of the autocorrelation functions in the time-momentum representation can be organized in the same way. A notable difference is that the contractions of the noise field give rise to additional diagrams, but since all these diagrams disappear in the limit ^ 0, it is clear that the autocorrelation functions (5.5) are given by

(t1,p1;...; t„,p„) = lim / e-i(uitl+-+untn) AUw1,p1; •••; w„,p„), (5.6) Ml0 J

where the autocorrelation functions on the right are those discussed in the previous sections. M

Note that the frequency integrals must be performed before is taken to zero, as otherwise one may run into infrared-singular intermediate expressions.

In view of its relation to the stochastic molecular dynamics, as expressed through eq. (5.6), and since the distribution (5.3) remains non-local at = 0, we are thus led to conclude that also the HMC algorithm is not renormalizable.

5.4 The Langevin limit

As already mentioned in subsection 2.1, the stochastic molecular-dynamics equation (2.5) is equivalent to the Langevin equation in the limit ^ to up to a rescaling of the simulation time.2 The associated n-point autocorrelation functions,

A°(w1,p1;...; w„,p„ )= lim (2^0)-nAlra(w1/2^0,p1;...; wra/2^0,p«), (5.7)

are known to be renormalizable to all orders of perturbation theory [1].

It may be instructive to see how exactly the renormalizability of the autocorrelation functions gets restored at one-loop order when is sent to infinity. To this end, first note that the Green function

lim K(u/2m>,p) = . , \ ,-2 (5-8)

uo^^ —iw + p2 + m0

assumes the expected form, which is smooth in position space except for a singularity at the origin. At one-loop order, the limit of the two-point function and the renormalizable parts of the four-point function is then easily determined starting from the identities (4.2), (4.8) and (4.9). These remain valid in the limit and only the tree-level autocorrelation

functions are replaced by the corresponding expressions involving the propagator

lim (2^io)~1G(uj/2^o,p) = 9 . 9 ,-oT? (5-9)

uo^^ w2 + (p2 + m2)2

and the Green function (5.8).

2 In lattice field theory, the Langevin limit can also be reached together with the continuum limit by

setting ^0 to some fixed value in units of the lattice spacing.

The contribution to the four-point function proportional to the integral J1 (which contains the non-removable ultraviolet singularity in the case of the stochastic molecular dynamics) however changes its character, because the integral

lim (2№)_1Ji(w/2№,p) = / I -iu + '2q2 + \p2 + 2mÛ

Jg I 2 J

1 \2 \// 1 , 1 * ,2 \ I I „ I

Q + ^Pj (5.10)

turns out to be absolutely convergent. In the Langevin limit, all ultraviolet singularities

at one-loop order are thus canceled by the parameter renormalization, as is expected to be the case in this theory [1].

6 Concluding remarks

The HMC algorithm is currently the preferred simulation algorithm in lattice QCD. In the past two decades, various improvements were included in this algorithm, many of them with the aim of reducing the computational effort required at small sea-quark masses (see ref. [20] for a recent review). Its scaling behaviour with respect to the lattice spacing has not received as much attention so far, but rapidly becomes an important issue when the continuum limit is approached.

While the dynamical properties of the HMC algorithm are well understood in free field theory [13], the situation in the presence of interactions tends to be rather more complicated. In particular, certain lattice artifacts (topology-changing tunneling transitions, for example, or unphysical critical points in the space of bare couplings) can cause large autocorrelations. The results obtained in this paper show that even in the absence of such effects there is no reason to expect that the HMC algorithm scales essentially as in a theory of free fields. Evidently, the non-renormalizability of the algorithm does not imply that it is invalid or unusable close to the continuum limit, but without further insight its scaling behaviour is unpredictable in interacting theories.

The HMC algorithm and the stochastic molecular dynamics may conceivably fall into the universality class of the Langevin equation. Independently of whether this is the case or not, it may be worth looking for renormalizable algorithms where the simulation time has scaling dimension less than 2. Eventually such algorithms might turn out to be faster than the HMC algorithm and they would have the advantage that their efficiency at small lattice spacings is predictable.

A Properties of the Green function K(t, x) A.l Definition

The Fourier transform (2.9) of the Green function is a smooth function of u and p that satisfies

\K(u,p)\2 < C (u2 + p2 + m0)-1 (A.1)

for some (mass-dependent) constant C. K(w,p) is therefore a tempered distribution and so is the Green function in position space. However, the Fourier integral (2.8) is not absolutely convergent and is to be understood in the sense of distributions. All these comments also apply to the propagator (3.3) of the basic field.

As a distribution, K(t, x) satisfies the wave equation

DK (t,x) = 5(t)5(x), (A.2)

where D is given by eq. (2.7). Since the polynomial representing D in frequency-momentum space is nowhere equal to zero, the Green function is the unique tempered distribution that solves eq. (A.2). In particular, one does not have the freedom of specifying retarded or advanced boundary conditions.

K(t,p) = I e-iUtK(w,p) (A.3)

A.2 Time-momentum representation

Using the residue theorem, it is possible to work out the Green function

= <,(i)e-»<5!2M £p=(p2 + m§-^)1/2. (A.4)

in the time-momentum representation. Note that ep is imaginary at small momenta if > m0, but the Green function always decays exponentially in the time direction. In the case of the field propagator, the identity

iwG(w,p) = K (w,p) — K (—w,p) (A.5)

and thus o

dtG(t, p) = K(-1, p) — K(t, p) (A.6)

may be used to show that

A/ x e-u0 |t| f . . sin(e„|t|)l

G i,p = -2--2 1 cosM) + № v pl 17 ^ • (A.7)

p2 + m0 L ep J

An immediate consequence of these results is that K(t, x) and G(t, x) become smooth functions of t at all t = 0 when smeared with a test function in x. Moreover,

K(t, x) ^ 0, dtK(t, x) ^ 5(x), (A.8)

as t | 0.

A.3 Explicit expression for K(t, x) in four dimensions

In dimension D = 4, the Green function in position space, K(t, x), can be calculated analytically. While the expression is of some interest in the context of our discussion of the non-renormalizability of the stochastic molecular dynamics in section 5, its exact form is not needed and the proof of the results quoted below is therefore omitted.

First note that the product ei№ K(t, x) is invariant under Lorentz transformations in 5 dimensions. As a consequence, and since the Green function vanishes at negative times, its support must be contained in the forward light cone t > |x|. Inside the light cone, i.e. at all t > |x|, the formula

e-Mot r cos(eos) £o sin(eos))

4n I s s J s=(t2_x2)1/2

holds, where e0 = (m2 — ,2)1/2. While the Green function has no singularities there, the formula shows that it diverges proportionally to

e_Moi(t° — x2)_3/2 (A.10)

along the light cone t = | x| .

M •n

B Proof of eq. (3.18)

First note that the integral

Ji(-u,p) = Ji(u,p)*, (B.2)

it suffices to calculate t/1(w,p) at all positive u. The imaginary part of the factor

(qp)2 + (2,o - iu)2(q2 + m0) (B.3)

is then strictly negative and the Feynman parameter representation

Ji(u,p) = 7(2,0 - iu>) ¡^ ck/.ch; f B)] e"^2+mo) (B.4)

2 J0 Jq

is therefore well defined.

The Gaussian integral over the momentum q can now be performed and leads to the formula

= o(A \D,2(^0 - M / dvdve-^mo 2(4n)D/2 Jo

x {u(z + ip2) + v}_1/2 {uz + v}_(D_1)/o , (B.5)

z = i(2,0 - iu)2 = 4,0u + i(4,2 - u2). (B.6)

2 I [(qp)2 + (2№ ~ M2(<?2 + m2)] (1,2 + m2} (B-1)

has the same divergent part as J1, because the integrand of the difference J1 — J1 falls off like (q2)-3 at large q and is therefore absolutely integrable in four dimensions. In view of the reality property

In eq. (B.5), the phases of the expressions in the curly brackets range from —\ir to and it is understood that the corresponding branch of their powers is taken.

The integral (B.5) is absolutely convergent for D < 4 and a holomorphic function of z in the half-plane Re z > 0. Its values along the real axis z > 0 can be computed by first substituting u ^ u/z and subsequently

u = ts, v = t(1 - s), (B.7)

0 < t < to, 0 < s < 1. (B.8)

The integral then factorizes into analytically calculable integrals and one finds that it is given by

This expression analytically extends to the half-plane Re z > 0 and therefore coincides with the integral (B.5) at all these values of z. Inserting eq. (B.6), the result

= {(2"° - (l+Vl+} (B10)

is thus obtained, where the branch of the square root with positive real part is to be taken. This formula holds for both positive and negative w.

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.

References

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

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

[3] A. Muñoz Sudupe and R.F. Alvarez-Estrada, Renormalization constants for the propagator of the stochastically quantized Yang-Mills field theory, Phys. Lett. B 164 (1985) 102 [SPIRES].

[4] A. Munoz Sudupe and R.F. Alvarez-Estrada, ß-function for Yang-Mills field theory in stochastic quantization, Phys. Lett. B 166 (1986) 186 [SPIRES].

[5] K. Okano, Background field method in stochastic quantization, Nucl. Phys. B 289 (1987) 109 [SPIRES].

[6] L. Baulieu and D. Zwanziger, QCD4 from a five-dimensional point of view, Nucl. Phys. B 581 (2000) 604 [hep-th/9909006] [SPIRES].

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

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

[9] 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].

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

[11] 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].

[12] S. Duane, A.D. Kennedy, B.J. Pendleton and D. Roweth, Hybrid Monte Carlo, Phys. Lett. B 195 (1987) 216 [SPIRES].

13] A.D. Kennedy and B. Pendleton, Cost of the generalised Hybrid Monte Carlo algorithm for

"""......................EC

free field theory, Nucl. Phys. B 607 (2001) 456 [hep-lat/0008020] [SPIRES].

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

r , . hd

[15] A.M. Horowitz, The second order Langevin equation and numerical simulations,

Nucl. Phys. B 280 (1987) 510 [SPIRES].

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

[17] M. Beccaria and G. Curci, The Kramers equation simulation algorithm. 1. Operator analysis, Phys. Rev. D 49 (1994) 2578 [hep-lat/9307007] [SPIRES].

[18] K. Jansen and C. Liu, Kramers equation algorithm for simulations of QCD with two flavors of Wilson fermions and gauge group SU(2), Nucl. Phys. B 453 (1995) 375 [Erratum ibid. B 459 (1996) 437] [hep-lat/9506020] [SPIRES].

[19] N.N. Bogoliubov and D.V. Shirkov, Introduction to the theory of quantized fields, Wiley-Interscience, New York U.S.A. (1959) [Intersci. Monogr. Phys. Astron. 3 (1959) 1] [SPIRES].

[20] M. Luscher, Computational strategies in lattice QCD, lectures given at the Summer school on "Modern perspectives in lattice QCD", Les Houches France August 3-28 2009

[arXiv: 1002.4232] [SPIRES].