Acceleration of gauge field dynamics

Volume 176, number 1,2

PHYSICS LETTERS B

21 August 1986

ACCELERATION OF GAUGE FIELD DYNAMICS Simon D U A N E Physics Department, Imperial College, London S W 7 2BZ, UK

Richard KENWAY, Brian J. P E N D L E T O N and Duncan R O W E T H

Departmentof Physics,

University of Edinburgh, The King's Buildings, Ma)field Road, Edinburgh, Scotland EH9 3JZ, UK

Received 19 April 1986; revised manuscript received 3 June 1986 It is shown that the success of acceleration for abelian gauge field dynamics need not depend on any choice of gauge. The discussion leads to proposing a particular scheme for acceleration in non-abelian theories which is also gauge independent.

1. Introduction As simulations of lattice field theories are performed nearer the continuum limit, the problem of critical slowing down becomes more urgent. This problem arises because the usual algorithms are constructed in a way which is independent of the physics of the system being considered. "Acceleration" is an attempt to tailor the dynamics of a simulation in some way to get around this. The general idea is easy to understand in the context of simulations based on molecular dynamics [1], for example of a real scalar field q,. The usual hamiltonian

essentially CPU time rather than physical time, which is represented by one of the euclidean lattice directions.) In the case of a free field having action =

-

(1.4)

the characteristic frequencies w(p) of the plane wave modes are given by ~02(p) = m 2 +p2.

(1.5)

defines a classical dynamics by the equations of motion

In simulations based on a discrete time version of this dynamics, the maximum possible step-size dt is inversely proportional to the highest frequency ~max, SO the minimum number of steps required for the slowest modes of the system to complete one cycle goes like

4, = 7r,

~0max/~min

H(4,, ~r) = ½~r2 + S(,~)

~" = - d S / d + ,

(1.1)

(1.2)

and time averages will coincide with ensemble averages F---(~- lim f r d t F(eo(t)), T ~ e~ ao

(1.3)

(F(q,)~ - fDq~ F(q,) e x p [ - S ( ~ ) ] , provided only that the conjugate momentum ~- is randomized occasionally from its distribution e x p ( - ~ r 2 / 2 ) [2]. (The time referred to here is

= [(?~2 + P .2. . . ) / ( m 2

{_Pmin 2 )J ,]1/2 -

(1.6) Critical slowing down occurs when this ratio diverges, i.e. when the number of degrees of freedom within one correlation length 1/m goes to infinity. This is a problem because all the interesting physics is associated with these slow modes. Suppose instead of (1.1) we take the hamiltonian to be H(q), .7r)=½'/r(M2-(} 2)

0370-2693/86/$03.50 © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)

l'rrq-S(~).

(1.7) 143

Volume 176, number 1,2

PHYSICS LETTERSB

21 August 1986

Then the free field frequencies w ( p ) are given by

namics is based on the hamiltonian

w2(p) = (m 2 + p 2 ) / ( M 2 + p 2 )

H(U. P) = ¼tr(P 2) -BS(U).

(1.S)

and critical slowing down is avoided completely by taking M 2= rn 2, so that all modes have a common frequency ,0--1. Mechanically, all we have done is make sure that the modes having a stiff potential (large spring constant) also have a large mass (inertia). The cost of this solution to the problem is that the equations of motion arising from (1.7), i.e. ~ ) = ( M 2 - 3 2 ) ~Tr, ~ = - d S / d ~

(1.9)

are more awkward to solve. (The most efficient method is probably to use Fast Fourier Transforms, hence the name, Fourier acceleration.) It is by no means obvious that this trick is helpful when interactions are introduced, indeed it is not even clear what value one should take for the parameter M 2, although intuition suggests that for a scalar field theory the renormalized mass-squared might be appropriate. Our prime interest here is in the extension of this idea to gauge theories. The essential feature of the simulation algorithm which allows the implementation of Fourier acceleration is that it should be based on a differential equation, so that the configuration evolves continuously. Whether it is based on molecular dynamics or the Langevin equation [3], or as we prefer on a hybrid between the two [2], is inessential. We choose the hybrid scheme because it includes the others as special cases and can be tuned to be more efficient than both [4,5]. For the sort of theoretical questions addressed here (gauge invariance and the effective equilibrium distribution) we find the hamiltonian formulation of molecular dynamics a little more transparent. Since the hybrid scheme is determined (up to tuning) by the hamiltonian, no generality is lost. Fourier acceleration for gauge theories has been discussed by Batrouni et ai. [6] in the context of Langevin simulations. The present paper grew out of our attempts to understand some of their comments. The starting point for our analysis in later sections is the (unaccelerated) molecular dynamics formulation of lattice gauge theory, which we recall briefly here. The simplest continuous time dy144

(1.10)

which leads to the equations of motion

~J=iPU,

]'=ifldS/dU,

(1.11)

where d S / d U is the derivative on the group manifold [7]. The coordinates U~(x, t) are group valued fields living on the links of a hypercubic lattice, between sites x and x +/~. The discrete "time" takes values n d t, where n counts the number of updates made. The conjugate momenta P,(x, t d t / 2 ) are algebra valued fields defined on the lattice sites living at intermediate times ( n 1/2) dt. Evolution through one time step involves the replacements

P,( x, t + dt/2) = P,(x, t - d t / 2 ) + dt P , ( x , t),

(1.12)

U~,(x. t + d / ) = exp[idt P.(x, t + d t / 2 ) ] Us(x, t),

(1.13)

where for example the usual Wilson plaquette action

S(U) =

Re{tr[V.(x,

t)U,(x + 9, t)

× Uf(x + ~t, t ) U / ( x , t)]}

(1.14)

gives P , ( x , t ) = - E 2 proj{V,(x, t ) V , ( x + b ,

t)

v~p,

xuy(x

t)U/(x,

+ uT(x -

t)u,(x -

x U,(x -

+

,)u;(x,

t) ,)

t)}, (1.1s)

and proj( X} - ½ ( X - X +) - ( 1 / 2 N ) t r ( X - X +) (1.16) in the case of SU(N). In this formulation of hamiltonian dynamics on Lie groups, there is no need to introduce Lagrange multipliers to preserve the constraint on U, which is satisfied by virtue of the linear projection onto

PHYSICS LETTERSB

Volume 176, number 1, 2

the algebra which occurs in the definition of the derivative with respect to U. It is the molecular dynamics version of the discrete Langevin equation based on exponential updates proposed in ref. [7]. A theoretical analysis of unaccelerated hybrid algorithms is contained in ref. [5], where a different formalism based on constraints and Lagrange multipliers [1] is used. As far as we can see, the introduction of Fourier acceleration makes the solution to the equations determining the Lagrange multipliers too expensive, and so we prefer to use the multiplicative updates in the present analysis. Corresponding to the local gauge invariance of the action S there is a symmetry of the dynamics. The gauge transformation is restricted to be independent of the extra time parameter t and acts as follows:

+,),

u.(x, t)--.g(x)U.(x, P.(x, t + dt/2)

(1.17)

g(x)P.(x, t + dt/Z)g+(x).

(1.1S) and from eqs. (1.15) and (1.16), (1.19)

g(x)P~(x, t +dt/2)g+(x)] +

= g ( x ) exp[idt P . ( x , t + d t / 2 ) ] ×

U.(x, t)g+(x+~).

(1.20)

After the introduction of occasional randomizations of P, the result is an algorithm which generates configurations from the usual Boltzmann distribution with integration measure dU d P. (Strictly speaking, this statement applies to the distribution of U and P at a common time, and not the staggered variables U(t) and P ( t - d t / 2 ) . This means that one has to be somewhat careful in the randomization of P ( t - d t / 2 ) , but mostly one is interested in the distribution of the U ' s alone. Technicalities of this sort are dealt with at length in ref. [5].)

2. Abelian gauge theory We begin by extending the acceleration prescription given in the introduction for a free scalar field to one for compact abelian pure gauge theory on a lattice. In the limit of infinite 13 this reduces to free lattice photons, and complete acceleration is possible. Instead of (1.10), we take

H= ½P~[(I--K)--,~ O2]-'P.+BS(U)

because P involves untraced plaquettes which transform like P, i.e. according to the adjoint representation. This is compatible with the evolution (1.13), because the result of a particular gauge transformation followed by an update is identical to the that of an update followed by the same gauge transformation, as shown in fig. 1 because

U (x,t)

exp[idt

21 August 1986

(2.1)

(02 = O, 0, where 0~ denotes the lattice finite difference operator) so that r = 0 corresponds to the original choice, and we shall show that ~ ~ 1 gives complete acceleration for free photons. Note that because the adjoint representation is trivial for an abelian group, P is actually'gauge invariant and so this hamiltonian still has the full gauge symme-

P (x,t+dt/2)

~

~

eidtPu

update

trans g augeIn g(x) g(x)U g+(x+~)~ .... g ( x ) P u g

1

+ (x)

"~.

eidtgPg+(X)g(x)U

g+(x+~)

Fig. 1. Commutativity of updates and gauge transformations in unaccelerated pure gauge theory. 145

Volume 176, number 1,2

PHYSICS LETTERS B with

try. We substitute Us = exp(ieAs) ,

(2.2)

and take the limit of infinite/3 = (2e2) -~. Then

/3S(U)--~½As(x)(Ss O -3s~ OZ)A~(x)

(2.3)

and the equations of motion are

}}s(x)=-[a s3-6s,a2]A~(x).

(2.4)

The hamiltonian describes coupled harmonic Oscillators, whose eigenmodes are plane waves of momentum q. Those with transverse polarization have w2(q) = q2/(1 - K + Kq 2),

(2.5)

while the longitudinal modes have zero frequency. The latter are precisely the unphysical gauge modes, which make no contribution to any observables. Complete acceleration occurs when ~ --. 1, since all the modes which contribute to physical quantities then have a common frequency. (The limiting value r = 1 may in practice be excluded: for example in the system with periodic boundary conditions the acceleration operator in (2.4) would have a zero mode.) For finite/3 the presence of gauge field self-interactions makes the effectiveness of this type of acceleration less clear. However, in the abelian pure gauge theory the effectiveness of acceleration is not dependent on gauge fixing. This is because, even with accelerated dynamics, updates commute with gauge transformations. A discrete time step with the accelerated dynamics is given by (1.13)

u

(x,t)

gauge ~ transn g(x)U g+(x+~)

P

P (x,t+dt/2)

replaced by ,~ 32] ~-~l]'s(y, t).

Since P~ and P~ are both gauge invariant and exp[idt P,(x, t+dt/2)] commutes with gauge transformations, the result follows as shown in fig. 2. It is therefore not true that one has to maintain a gauge condition when accelerating pure gauge theory. Now suppose we introduce charged matter fields, for example Higgs scalars living on the sites of the lattice:

S(U, ~) =/3S(U) + (Dsq~) + (Ds~),

(2.7)

where the covariant difference D s is defined by

Ds~(x) =

U~(x)~(x+~)-~(x).

If we try and improve on the simplest kinetic energy for ~ by replacing ~r+Tr --* ~r+ [(1 - K ' ) - , ~ '

32] -'~..

(2.8)

then this term and hence the new hamiltonian is no longer gauge invariant. It is essential to fix a gauge in this case. Another way to cure this would be to replace the ordinary difference in (2.7) by a covariant difference. This restores gauge invariance but would spoil the U distribution. In fact, integrating out the 7r's would modify the effective probability distribution to det[(1 - ~') - ~:'D 2 ] e x p [ - / 3 S ( U ) ] .

(2.9)

A further difficulty is that while Fast Fourier Transforms make it relatively cheap to apply the

eitPU I , eidtPg(x)U g+ (x+~)

Fig. 2. Commutativityof accelerated updates and gauge transformations for abelian pure gauge theory. 146

(2.6)

Y

>

update

P,(x, t)

~2[(1 - , ~ ) -

(x,t+dt/2) ~

21 August 1986

Volume 176, number 1,2

PHYSICS LETTERS B

inverse of an operator based on the ordinary difference, they are no help when covariant differences are involved. Both of these problems arise in the nonabelian pure gauge theory and in the next section we suggest a method for avoiding them. A minor awkwardness when it comes to discretizing time for these equations is that ,# and rr naturally live at staggered times t and t - dt/2, and U does not naturally live at both of these.

3. Nonabe#an gauge theory Yang-Mills particles differ from photons in that they carry charge themselves and so the difficulties found in trying to accelerate the dynamics of charged matter arise already in accelerating the dynamics of the pure gauge theory. However, the need for acceleration is so great that it is worthwhile considering how to overcome these difficulties. The simplest way is to introduce extra scalar fields q,~, in the adjoint representation, one for each P~,, so that the total hamiltonian becomes H=Jtr{Pu[(1-x)-xD

2] ~ P , ) - ~ S ( U )

21 August 1986

dynamics hamiltonian ~, (m 2 ,,+

)-l~r~ + ~+~b

(3.2)

leads to the (square of the) fermion determinant without a prohibitive computational cost: each update of the system requires one sequence of conjugate gradient iterations (per sweep). The extra doubling of the number of flavours, which was originally introduced to make the hamiltonian real, was subsequently removed [9] by restricting the scalar to live on every other site, which is consistent with staggered fermions because the operator in (3.2) does not mix even and odd sites. This system can be simulated using the hybrid algorithm, randomizing either % or ~ (in addition to P). The randomization of % is achieved by generating uncorrelated gaussian noise 1~ on all sites and taking % to be the projection onto even sites of (m + D)~. Thus one avoids a pseudofermion style Monte Carlo by taking advantage of the factorization

mZ--O2=(m+~)+(m+O).

+ ¼tr{ %[(1 - x ) - xD2] ~r, } + ¼tr(q}uga,). (3.1) Now the integral over the extra fields ¢r introduces an inverse determinant which cancels that from the gauge momentum integral, and the equilibrium distribution is correct. As always, ergodicity is guaranteed in the hybrid algorithm by the occasional randomization of one from each pair of canonically conjugate variables. For the gauge fields this is complicated by the fact that the distribution of P ' s depends on U in a non-trivial way, and it would appear necessary to randomize the P ' s using a Monte Carlo, somewhat analogous to that in the pseudofermion method. This is straightforward using the change of variable P~ = [(1 - ~) - KD2]V~. There is a superficial similarity between the determinant which results from the integral over P,, to that which represents the effects of dynamical fermions (however (2.8) lacks any spin structure). Polonyi and Wyld [8] showed how introducing an extra bosonized fermion with molecular

The operator which appears in the gauge kinetic energy term in (3.1) is not factorizable in the same way, so we are apparently left with no alternative but to generate V by Monte Carlo. In contrast to that in the pseudoferrnion algorithm the randomization here is only required every 100 sweeps or so. This permits the Monte Carlo to be continued for long enough to reach equilibrium for the P ' s in the given U background. For the scalar fields it is simpler to randomize the ,/,'s. The equation of motion for the scalar field is straightforward but unaccelerated: & = [(1 - K) - xD~D,] rr , i = - q'~,"

(3.3)

The equation for the gauge field involves an inversion somewhat similar to that required in simulations which include dynamical fermions [8], for which conjugate gradient iteration is tolerably efficient: U~ = {i[(1 - K) - KD~D.] -IPo } U~,

P, = i d n / d U , .

(3.4) 147

Volume 176, number 1,2

PHYSICS LETTERS B

Thus the o v e r h e a d associated with trying to accelerate the gauge d y n a m i c s in this w a y w o u l d be relatively less severe for the full system including fermions, i.e. two conjugate g r a d i e n t iterations i n s t e a d of one. It is unclear how effective this acceleration scheme would be as neither the scalar field dyn a m i c s nor that of any d y n a m i c a l fermions is accelerated. Since it a p p e a r s that fermion effects are of the o r d e r of 10% in the case of Q C D , this partial acceleration m a y still be useful. As an alternative to m a k i n g further m o d i f i c a t i o n s to the d y n a m i c s in o r d e r to perfect the acceleration, it m a y be simpler a n d as useful to find some p r e c o n d i t i o n i n g of the o p e r a t o r s used in the conjugate g r a d i e n t iterations which speeds convergence [7]. A n y m a n a g e a b l e p r e c o n d i t i o n i n g m u s t be b a s e d on the o r d i n a r y difference, a n d this would only be helpful if the c o n f i g u r a t i o n were p u t into a " s m o o t h " gauge. T h e fact that the d y n a m i c s b a s e d

148

21 August 1986

on (3.1) is c o m p a t i b l e with the gauge s y m m e t r y m e a n s that any such gauge c o n d i t i o n can be imp o s e d without spoiling the d i s t r i b u t i o n of gauge invariant quantities. W e t h a n k K e n Bowler a n d Christine Davies for useful discussions.

References [1] D. Callaway and A. Rahman, Phys. Rev. Lett. 49 (1982) 613. [2] S. Duane, Nucl. Phys. B257 [FS14] (1985) 652. [3] G. Parisi and Wu Yongshi, Sci. Sin. 24 (1981) 483. [4] J.B. Kogut, UIUC preprint. [5] S. Duane and J.B. Kogut, in preparation. [6] G. Batrouni et al., Phys. Rev. D32 (1986) 2736; A.S. Kronfeld, report Wuppertal meeting (1985). [7] I.T. Drummond et al., Nucl. Phys. B220 [FS8] (1983) 119; A. Guha and S.C. Lee, Phys. Rev. D27 (1983) 2412. [8] J. Polonyi and H.W. Wyld, Phys. Rev. Lett. 49 (1983) 2257. [9] J. Polonyi et al., Phys. Rev. Lett. 53 (1984) 363.