- Email: [email protected]
Loop dynamics for gauge theories: A numerical algorithm
Nuclear Physics 13239(1984) 135-160 © North-Holland Pubfishing Company
LOOP DYNAMICS FOR GAUGE THEORIES: A NUMERICAL ALGORITHM
G. MARCHESINI Dipartimento di Fisica, Universiti~ di Parma, INFN, sezione di Milano, Italy
Received 31 October 1983
The loop equations defining the lattice U(N) gauge thco~ are recalled and a formal solution is presented (for N = 1 and N --* zc). Wilson's expectation function W(C) for a loop C is expressed as a matrix element of a resolvent in the space of all the loops between C and 0 (no loop). It is shown that such a solution provides a new numerical algorithm to compute physical quantities. This is based on a Ulam-von Ncumann type of method for computing inverse matrLx elements by introducing importance sampling for the paths in the space of the matrix indices, which in this case are the loops. As a result W(C) is obtained by summing over important paths in the loop space connecting C to 0. A Monte Carlo program is presented, for the N ---,~c case, where a ve~' simple form for the importance sampling is introduced so that the computer time for each step in the construction of the path is minimized.The rates for successful paths (i.e. path C --*0 within a given finite number of steps) are computed for D = 2 and D = 4. Both rates and computer time involved encourage us to attempt a large scale calculation. Here the numerical studies of the convergence and of the fluctuations are presented only for D = 2. Convergence is rather fast, but specially in the weak-couplingregion rare and large fluctuations appear thus suggestingthat a bctter tuning for the importance samphng is needed.
1. Introduction T h e present numerical algorithms [1] for c o m p u t i n g physical quantities in (lattice) gauge theories [2] are based on the Metropolis m e t h o d in which the gauge fields are treated as d y n a m i c a l variables a p p r o a c h i n g the e q u i l i b r i u m in a Markov process. W i t h i n this method one meets with difficulties of two kinds: (i) since the gauge fields are r e d u n d a n t variables one does not k n o w a priori whether two different field configurations, generated in the M o n t e Carlo process, are physically different or rather connected b y gauge rotation, a n d this may reduce the efficiency of the method. ( i i ) A more i m p o r t a n t difficulty from a practical point of view is that the l i m i t a t i o n in the n u m b e r of the d y n a m i c a l variables implies that the calculations can be d o n e in a space-time world extremely small. Because of these difficulties it seems desirable to find a n d test alternative numerical algorithms which could provide a different perspective a n d eventually overcome some of the difficulties. T h e physical quantities in (pure) gauge theories are the Wilson loops [2]. A n a t u r a l f o r m u l a t i o n of gauge theory is then in terms of loop variables rather then in terms of 135
136
G. Marchesini / Loop dynamics
gauge fields. Recently some progress [3-8] in this direction has been achieved. In particular for the U ( N ) gauge theory dynamical equations for the Wilson loops have been obtained. They correspond [9] to the Schwinger-Dyson equations of the theory. For the two extreme cases N = 1 and N ~ oc these equations involve single Wilson loops only. In this paper we obtain a formal solution for these loop equations and show that this solution provides a numerical algorithm to compute physical quantities. For the moment we discuss in detail only the N = 1 and N ~ ~ cases. The formulation is in the space of loops (which lie on the space-time lattice): for instance in QED a matrix element K c , c is defined which relates a given loop C to a next-neighbour loop C'. The Wilson expectation value W(C) for a loop C is expressed in terms of the matrix element of the resolvent of K between C on the zero loop configuration. In spite of such an exotic matrix inversion problem, numerical algorithms can be proposed. In this paper we explore a method based on the Neumann expansion, where one has to compute the matrix element of K" between the loop C and the zero loop. The intermediate indices of this power of K are of course loops. The matrix element of K" can then be viewed as given by the sum over all the paths in the loop space which connect the initial loop C to zero with steps defined by K. The number of paths in loop space is very large so that this sum can be practically done only by an importance sampling [10] method. Note that although the number of loops is large, with such types of algorithm one deals for each step with one or few loops so that even for a very large space-time world there is no problem for computer memory; in fact a loop is a variable in one dimension, the proper time. To study this method we have devised a Monte Carlo program for the U ( N ) gauge theory in D dimensions for N ~ oo. The most crucial point for the efficiency of the method is of course the choice of the a priori probabilities for the importance sampling. In fact one has to maximize the rate of successful paths i.e. of paths in the loop space which start from the original loop C and shrink to zero in a finite number of steps. For the moment we have made the most simple choice by defining a priori probabilities depending on the individual steps, so that the shrinking of a loop can be favoured only in a limited way. In this way one minimizes the computer time for each step in the development of a path, but one expects a low rate of successful paths especially for large loops (actually we will show on the basis of the formal solution and of the loop equations that a large-loop calculation can be subdivided into smaller parts). A first aim of our numerical analysis is to study these rates: the result is that even in D = 4 the rate of successful paths is relatively high at least for small loops. A second part of the numerical analysis is devoted to the study of the convergence of the algorithm both in the strong- and weak-coupling region. We shall limit ourselves here to report the results for D = 2. In the strong-coupling region the convergence is rather fast and the fluctuations are small. In the weak-coupling region instead there appear very rare (with a relative frequency of order 10 -~) but large fluctuations of indefinite sign which obscure the convergence (in fact if these
G. Marchesini / Loopdynamics
137
fluctuations are artificially removed a rather fast convergence is obtained). This situation is typical of an importance sampling strategy not perfectly tuned. In sect. 2 we set the notation apt to describe the loops in a D-dimensional space-time lattice and for any given loop we define the set of its next-neighbour loops. In terms of this notation we recall the form of the loop equations for the U ( N ) gauge theories and discuss the specific cases N -- 1 and N --o ~ , the planar theory. In sect. 3 we introduce the space of loops and present the formal solution of the loop equations for N = 1 and N ~ oo in terms of a resolvent in the loop space. We describe here how from the Neumann expansion of this resolvent a numerical algorithm can be constructed. In sect. 4 we describe the details of the specific algorithm we have constructed for the planar D-dimensional theory, and report the results on the analysis of efficiency and convergence. In sect. 5 some conclusions are drawn. In appendix A we obtain the solution of the loop equations by using the formalism of stochastic quantization [11, 9] which may be useful especially to study the extension of our method to the N-finite case. Finally in appendix B we discuss the single-plaquette planar model.
2. Loops and loop equations The loop equations set relations among different Wilson loop expectation functions. In this section we recall the form of these equations for the U ( N ) gauge theory and the simplifications which occur for N = 1 and N -o oo. Before this however, (i) we shall set the notation to describe loops in a D-dimensional space-time lattice and (ii) we shall define the set of next-neighbour loops of any given loop C.
2.1. NOTATION AND NEXT-NEIGHBOUR LOOPS
Consider an oriented closed curve (loop) C in the space time lattice
c:
(2.1)
are the lattice points in D dimension (the lattice unit is set to be one). Two successive points x k , X k + 1 belong to a link in the lattice
x = (nln2...no)
Xk - 1 = Xk + e,,
(2.2)
where e~, are unit vectors in the positive and negative directions:/~ = _+1, _+2 , . . . , _+D. Of course C has no spikes i.e. for any k xk+ 1 ~ Xk_ 1 .
(2.3)
G. Marchesini / Loop ~lynamics
138
XK Fig. 1. Transformation C --, C ' ( k , ~,) in eq. (2.4) for the link x k, Xk_ v
A first set of loops next-neighbour to C is obtained by considering the following two classes of elementary deformations of C: (a) C ' ( k , u) is a loop obtained from C when the two successive points of C, x k and Xk. ~ ( = Xk+ e~,) are substituted by the set of points (cf. fig. 1) X k ' X k + l =e, X k ,
X k Jr- elL,
x k -[- elz ~- e~,,
x k-'F e u,
Xk~
Xk! l,
(2.4) where e~ is in any direction orthogonal to e,, and 1 ~< k ~
xk+e~,
xk+e~+e,,,
xk~ t.
(2.5)
Note that such a deformation of C may generate spikes. In this case C-(k, ~,) is actually obtained by removing them. In the case in which C has intersections one introduces another set of nextneighbour loops obtained by splitting C. The following possibilities are considered. (a) Two points of C x j a n d x k ( j < k) have a link parallel intersection when Xj = X k ,
Xj+ 1 = Xk !
1"
(2.6)
An example is given in fig. 3 where the link is intersecting only twice. In this case the loop C can be split in two loops
% - (Xj.,...Xk), Ckj = (x x ...XjXk. 1 . ..x,,c).
(2.7)
(b) Two points of C xj and x~ ( j < k) have a link antiparallel intersection when (fig. 4) xj = Xk,
t
xj÷ 1 = x~ _ ~.
1
Fig. 2. Transformation C --, C - ( k , v) in eq. (2.5) for the link xk, xk, 1.
(2.8)
G. Marchesini / Loopdynamics
v
139
J
~kj
Fig. 3. Example of a splitting for a parallel intersection case: C ~ CjkCtj.
As shown in fig. 4 the loop C can be split in Cjk~-~ ( X j + l " " X k
-- 2)'
Ckj =-- (X1...Xjxk+l...XnC ).
(2.9)
These two new loops may have spikes which have to be removed. Intersection of sites, but not on links (cf. fig. 5), are not considered for the loop equations. The inverse operation of combining two loops to make a single loop is relevant only for N ~ 1 and finite and is not considered here (see appendix A).
_5
tI
),
i --° v
,k
i
Fig. 4. Example of a splitting for an antiparallel intcrscction case: C --* CjkCg ~.
Ii Fig. 5. Example of an intersection only on the site. This case does not contribute to the loop equations since the intersection is at 90 degrees.
140
G. Marchesini / Loop dynanlics
2.2. LOOP EQUATIONS
One is now in a position to recall the loop equations for the Wilson function 1
/
(2.10)
where U(C) is the usual U(N) gauge field associated to the loop C, i.e. the ordered group product of link variables around the loop. The expectation value in (2.10) is obtained for instance by assuming the Wilson form for the action S(U):
S ( U ) = ~ - E Tr(Up+ U~),
(2.11)
P
where the sum extends over all elementary squares or plaquettes of the lattice and Up is the usual plaquette field. For this form of the action the loop equations are nc
1
E
0 = H{ W(C)} = ncW(C) + X
k=l
w(c (k,.))} v
+ 2 E sjkW(Cjk,Ckj),
(2.12)
j
where W ( C I ' C2) -
N1
TrU(C1 ) ~1T r U ( C 2 ) ) •
(2.13)
In (2.12) the notation of the previous subsection is used. Here n c is the number of points of the loop C; the sum over u is restricted to the directions orthogonal to the link xk, Xk~ 1 (cf. eqs. (2.4) and (2.5)). Finally in the last term Sjk = + 1 for a link parallel intersection at the points xj and Xk; while Sgk = --1 for a link antiparallel intersection (of course sj~ = 0 otherwise). The boundary condition for these loop equations is that when C is shrunk to zero U(0) = 1 so that
W(0) = 1.
(2.14)
Eq. (2.12) relates the single-loop expectation function to the two-loop expectations in (2.13). However for N = 1 and N ~ ~ these loop equations involve only single-loop expectations. In the case of QED the gauge field is a phase and onc finds in (2.12)
w(Cjk, Ck ~) = W(C),
(2.15)
G Marchesini / Loop dynamics
141
so that the loop equations become rt C
(nc+2Esjk ,
j
W(C)+~ E ~
k=l
E'{W(C=(k,P))-W( C
(k,u))}=0.
u
(2.16) It is interesting to note that, if the loop is at fixed time, the coefficient of W(C) is just the electric energy of the loop. In fact ec = n c + 2 E Sjk = E ( q [ j
-- q[)2,
(2.17)
I
where the sum over l extends to all the links which are crossed by the loop C and qt ~ is the number of times the loop crosses the link I in one or the other direction. The interpretation of this fact is clear in the formalism of appendix A. In the case of the U(N) theory for N --* ~ the simplification of the loop equations arises because of the factorization property [12] W(C1,C2) = W(Ct)W(C2)(1 + O ( 1 / N 2 ) ) .
(2.1.8)
3. Wilson's loop and paths in loop space In this section we present the formal solution of the loop equations (both for N = 1 and N ~ oo) in terms of matrix elements of a resolvent in the space of loops. We then show how the Neumann expansion of this matrix element can be viewed as obtained by summing over paths in the space of loops, and illustrate the possibility of performing this sum by importance sampling methods [10]. Here we simply discuss the general features of the method and leave to sect. 4 the presentation of a specific algorithm. 3.1. F O R M A L
SOLUTION
To a given loop C in (2.1) we associate a ket IC) with the normalization (C'[C) = 6(C',C),
(3.1)
where 8(C',C)= 1 if C and C' coincide both in points and orientation, and 8 ( C ' , C ) = 0 otherwise. We denote by 10) the state corresponding to no loop; by [C1,C2) the state corresponding to the two loops C 1 and C2; and so on. In appendix A we show, within the formalism of stochastic quantization [11], that the Wilson loop expectation W(C) can be expressed as the matrix element W(C) = (01WIC),
(3.2)
G. Marchesini / Loop dynamics
142
and the expression for W is obtained. Here we simply give the form of W for N = 1 and N ~ oc and show that the loop equations are satisfied. Consider first the simpler case of QED where the loop equations are given in (2.16) with the boundary (2.14), and where, due to (2.15), only single-loop states need to be introduced. In this case W QE° is given by the resolvent of K Q ~ :
W QED= (1 - KQED) -1 ,
(3.3)
where K QED is defined as follows (for notation see sect. 2)
KQEDIO) = 0,
KQ~IC)= -
(3.4)
E~, IC)+x-~c F, E'[IC-(j,~))-IC*(j,~))]Cj
]
j'~l
(3.5)
v
In fact, due to (3.4), the boundary condition W(0) = (01(1 - K Q E D ) - l l 0 ) = 1,
(3.6)
is satisfied. Moreover, from the form (3.3) one has (01 wQLD(1 -- KQED)IC ) = (01C) ,
(3.7)
and because of (3.5) this is just the loop equation for W(C). Consider now the N--, ~ theory. In this case one needs to introduce also the many-loop states since the expectations W(Ct, C2) in (2.13) are involved. In appendix A it is shown that for a general U ( N ) theory
w(c) = (01 wl c), w ( c ~ , c : ) = (Ol WlCl,C,).
(3.8)
For N ~ ~ moreover the factorization property in (2.18) implies (01 w l c t , c 2 ) = (01WlC1)(Ol WLC2)(1 + O(1/N2)),
(3.9)
which sets the property of the action of W on a multi-loop state in terms of the action on single-loop states. Analogously to QED one requires W of the form W = (1 - X ) - l ,
(3.10)
143
G. Marchesini / Loop dynamics
where XlO ) = O, i, )>1 .
KIC) = - - - 2 ~ sjilCij,Cj~) + 1 - -
nc
hnc j~l =
j
,
(3.11) As in the previous case one show that with this form both the boundary conditions and the loop equations are satisfied. 3.2. Q E D
To obtain some familiarity with a resolvent in a loop space consider the Neumann expansion* of (3.3) (the index QED is here dropped)
w ( c ) = <01(1 - x ) - ' l c > oC
=
'Koc,,...Kc2c,Kclc,
(3.12)
n = 0 C1C2 ... C,,
where the sum is restricted to non-zero loops (C i ~=0, i = 1. . . . . n) and the matrix elements of K are given by (3.5): the non vanishing ones are
Kc, c = (C'IK[C) =
1 - ec/nc,
C'=C
Yr 1 / h n c ,
c,=
(3.13)
This matrix then connects C to itself and to the next-neighbour loops C - ( j , p) introduced in sect. 2. ec is defined in (2.17). According to (3.12) W(C) is obtained by summing over all possible ways in which the initial loop C is shrunk to zero via the elementary steps in (3.13); these will be called paths C -~ 0 in the loop space. This picture suggest a possible way for computing W(C) by using numerical algorithms of the type proposed by Ulam and von Neumann [10]. On the basis of (3.13) one defines an importance sample for the paths C --, 0 and W(C) is obtained as an expectation value corresponding to (3.12). The crucial point for the efficiency of such an algorithm is of course the definition of the importance sampling. In fact the matrix (3.13) connects a loop C to * This corresponds to the strong-coupling expansions. To study the weak-coupling region one needs a representation for the resolvent which allows the analytic continuation to small ~. This wSll be discussed in sect. 4, and not here, so as to leave the argument simple.
144
G. Marchesini / Loop dynamics
4(D - 1)n c modified loops, in D space-time dimensions, so that quickly the number of loops becomes considerably large, and one may expect that in a finite number of steps only a negligible fraction of paths leads to a zero loop. In order to obtain a contribution to W(C) one has then to favour those steps in (3.13) which lead to smaller loops, that means in general the transition C ~ C - ( j , u). Note moreover that the contribution for these transitions is positive and in general larger. This explains why "although the matrix elements in (3.13) have not a definite sign, the value of W(C) is positive for not-nested loops. To show that the transition C ~ C - ( j , u) leads to larger contributions it is convenient to isolate from K in (3.5) its diagonal part K D. Introducing then K ' --- (1 - KD)- I(K - KI)) one finds
W(C) = <01(1 - K ) l I C > = n~-ff:c ec @1(1 - K ' ) - l l C > '
(3.14)
where K~,c = T n c ,
1
C,=C±(j,v)"
(3.15)
£C' ~?/C '
Note then that (nc,/ec,) is in general larger for C' = C - ( j , u) than C ' - C ~(j, p). Anyhow if the loop C is large, this may not be enough to obtain a considerable fraction of paths C ~ 0. However the loop equations may be used in this case to overcome this difficulty by making shorter the minimum required path. In fact observe that any path C---, 0, before reaching the zero loop, must reach a loop around a plaquette (C = P). The path C ---, 0 may then be divided in the first part C ~ P and then P ---, 0. This last step leads to the expectation function for the single plaquette W1 = W (C = P), which can then be factorized out, and one needs only to study the paths C ---, P where the plaquette is reached only in the last step. This may be achieved formally by using the decomposition K = Kp + K, Kp ~ E *K[ P)(P[,
(3.16) (3.17)
P
where IP) is the state corresponding to the plaquette loop and the sum is over all oriented plaquettes. The matrix K acts as in (3.5) except that the subspace of plaquette loops is projected out: R IP) = 0.
(3.18)
G. Marchesini / Loop dynamics
145
By using (3.16) one finds in fact WIG) = ( 0 1 ( 1 - K ) - l f C ) = W l ~ * ( P l ( 1 - g ) - ~ l C ) .
(3.19)
P
Due to (3.18) the plaquette states will never appear as intermediate states in the Neumann expansion of (1 - K)-a. This means that one needs to consider now only the paths C --* P (for any P) which come to the single plaquettes only in the last step. This procedure can of course be iterated by projecting out of K also the subspace of states next-neighbour of the single plaquette and so on. In this way the calculation of a large Wilson loop expectation may be split into calculations involving smaller steps. 3.3. U(N) THEORY FOR N--+
The paths related to the Neumann expansion of W(C) for the U(N) theory are far more complicated than in the previous QED case. In fact here one has to take into account also disconnected loops. Due to the form of K in (3.11) starting from IC) after a few steps one may generate a state corresponding to two loops, say IC 1, C2). The subsequent evolution of this state is complicated by the possibility of an interaction of the two loops to make a single loop again (see appendix A eq. (A.6)). This interaction term however is of order 1/N 2, so that for N ~ oe the factorization (3.9) holds and the two loops move on independent paths. As a result the Neumann expansion of W(C) corresponds to the sum of all paths which branch as tree and each branch independently shrinks to zero. This can be seen also formally by singling out from K in (3.11) the splitting part:
K=K'+K~,
(3.20)
where KsIC ) -
(C']K'IC) =
2
~
nC j
Sii]Cji,C,j),
-T-1/)~nc,
(3.21)
C ' = C~-(j, v).
(3.22)
1Ks)-i(1 - K ' ) IIC ) .
(3.23)
In fact one can write
w(c) =
(01(1 - x ) -11 c )
= (0](1 - ( 1 - K ' ) Introducing the "propagator"
D(C', C) = (C'1(1 - K ' ) - 1]C),
(3.24)
G. Marchesini / Loop dynamics
146
one finds from (3.23)
Esj~W(Cj,)w(c,i)D(C,C),
w(c) = D(0,C)- 2E !
C' nC' j < i
(3.25)
where the sum extends over all possible splittings c ' --, c ; i % . The tree equations (3.25) is characteristic of the N ~ ~ limit: after the branch C ' ~ C)'iCi~ the paths C,'j ~ 0 and Cj; ~ 0 are independent. As in the QED case, when the single-plaquette function W1 = W (C = P) is known, the value of W(C) can be obtained considering only the paths C ~ P which go to a single plaquette state IP) only in the last step. In fact as in QED by projecting out the single-plaquette subspace in K one obtains
W(C) = W t ~ *(P[(1 - K)-IIC),
(3.26)
P
where the sum is over all oriented plaquettes and in K the subspace of plaquette-loops is projected out: K = KE~(1 - IP)(P I). Moreover from the factorization property in (3.8) one has Wl(1 - l ~ ) - l l C l , C 2 ) = WI(1 - K ) - l i c 1 ) w t ( 1 - K ) - I I C 2 ) .
(3.27)
These two equations give an alternative structure of the paths needed to compute W(C) if W~ is known: one has to consider all the paths which may branch and each branch ends when a plaquette state is reached. According to (3.26) and (3.27) each branch gives a factor W1 so that, in contrast to the QED case, W(C) is not proportional to W1. Note finally that the splitting,possibility, absent in QED, considerably increases the fraction of successful paths since smaller loops are generated from larger ones.
4. Numerical algorithm and first results In this section we present the Monte Carlo program we have constructed along the lines of sect. 3 for the calculation of Wilson's loop in the planar theory (N ---, ~ ) . To be able to explore the weak-coupling region we introduce a relaxation parameter. As explained earlier the crucial point for the method is the choice of the importance sampling, and we describe in detail the expected feature of our choice. We then report first the numerical results aimed (i) to study the rate of successful paths in D = 2 and D = 4; (ii) to study for D = 2 the convergence of the algorithm both in the strong- and weak-coupling region, together with the role of fluctuation.
G. Marchesini / Loopdynamics
147
4.1. THE ALGORITHM
As previously explained the expansion of the resolvent in (3.26) corresponds to the strong-coupling expansion and therefore, to study also the weak-coupling region, one needs on analytical continuation. For this reason we introduce a relaxation parameter r. This can be done by using the following expression for W(C) which is equivalent to (3.26):
W(C) = G }--'.*(P[ (1 - A )
11C) ,
(4.1)
P
where A, which replace the matrix K, is defined as
AI0) =0, A IC) = (1
-
r -
rK)l
C),
(4.2)
and for the factorization (3.21) WI(1 - A ) - I I C I , C 2 ) = WI(1 - A ) - I l C t ) W ~ ( 1 - A)--llCz).
(4.3)
The non-vanishing matrix elements of A are then
Ac,c =
{C'IA
IC> =
1 - r, r T ~n---7 ,
C' = C C' = C ! (
j, v),
(4.4)
and for the splitting
(Cji, CijlA IC) =
2r
- --sj,. r/c
(4.5)
The representation (4.1) holds for any value of r and by taking r = O(~) this may be used to extend the convergence of the Neumann expansion of W(C) to the weak-coupling region. In appendix B we test the convergence for the single-plaquette planar model. In spite of the fact that this model undergoes a phase transition [13] even at very weak coupling the expansion is quickly convergent when r = O(~). We come now to discuss the most important point for the Monte Carlo program: the a priori probabilities for the importance sampling. For these probabilities we made a very simple choice: they are defined for each step of the path and are largely independent of the structure of the loop to be transformed. Specifically we introduce p _+as the total probability for the loop transitions in the class C --, C -(j, v); p~ as the total probability for the loop splitting (of course if the loop has no intersection p~ = 0). Within each class the probabilities for the individual transitions are equally
G. Marchesini / loop dynamics
148
distributed. Finally the probability for no transformation, P0, is defined by the normalization
(4.6)
p+ + p _ + p.~ + p o = 1.
W(C) is now obtained as an expectation value for the product of the matrix element of A accordingly normalized. In table 1 in correspondence with the class of loop transitions we report the multiplicities of the class, the total probabilities, the related matrix elements and finally the corresponding factors to be used in the expectation values. Of course the last row is absent when C has no intersections, in this case, according to (4.6), P0 = 1 - p + - p _ . The procedure for computing W(C) is then the following: starting from the loop C and using the probabilities in table 1 one generates a new loop C'. The process is then repeated for C'. Eventually after a few steps the loop splits. In this case one continues the evolution for the two loops separately ( N - , m). The evolution of a branch ends when a plaquette loop is reached, and the full process (path) ends when all branches end. We will call successful a path which ends within a given finite number of steps. Suppose that one of such paths has the following rates for the related transformations: N . : C --> C ~-:(j, p), N/: C ~ Ci, C 0
(Sii = + 1 ) ,
N,":c-,cj,%
(,j,=-l),
No': C ~ C
(C with intersections),
No": C ~ C
(C without intersections).
(4.7)
According to table 1 the contribution to W(C) from such a path is TAm.~:1 Loop transitionsand dynamicalfactors Transitions
Multiplicities
Probability
Matrix element
Factors
C--,C C--*C=(j,: ,)
1 2(D-1)n c
Po p~
1-r T-t - ~,nc
(1 r)/p o _T_2(D-1)r ?tp r
Mc = ~ 1.5,1
P~
C --, Cj,Co
j
- 2rSji
nc
- 2rSJ~' Mc Ps nc
G. Marchesini / Loopdynamics
w(path) = W~%(
1-r
)N (1-r
1 -P+-P--Ps
149 NIT
1 -p+-p_
×
(4.8) Here Nr~ = 1 + N / + N~" is the total number of branches and each branch brings a factor W1 according to (4.1); and n c , Me, are the length and number of intersections of the loop C i that underwent the splitting. Finally the value of W(C) is given by
W(C) =
lim -~path -~ OC
1
• w(path).
Npath p a t h
(4.9)
According to (4.8) and (4.9) W(C) is computed when W1 is known. However even Wt can be computed with this method by using the loop equations. In fact they set a relation between W1 and the two-plaquette loop expectations. By computing these last functions with (4.8) and (4.9) one obtains a relation to determine W~. 4.2. DISCUSSION ON THE ALGORITHM (a) Owing to our choice of the a priori probabilities (i.e. satisfying the Markov property) the computer time for each step in the path is minimized, but the shrinking of a loop can be favoured only in a limited way: we require p _ > p + in order to inhibit the formation of very nested loops. Because of the expected low efficiency of this importance sampling we have used the formulation which requires path C ~ P rather than C ~ 0. However, in order to compute large Wilson loops, one probably must split the calculation into smaller steps by using the loop equations along the lines described in subsect. 3.2. (b) As clear from (4.8) the sign of w(path) is not definite, therefore in the sum (4.9) one expects cancellations. These cancellations are not important in the strongcoupling region where the paths are dominated by the ones in which at each step the transition is C ~ C (j, p). In the weak-couping region instead one expects cancellations since these paths are not so dominant. (c) Moreover, within this formulation, one also expects large fluctuations in the weak-coupling region although rare. In fact in this region one requires a small relaxation parameter so that the transition C ---, C contributes with a factor (1 - r)/po sizably larger than one. However when the numbers in (4.7) for the various transitions are around their average values the large factors coming for C--, C transitions are compensated by the factors coming from the other transitions (for
G. Marchesini / Loop dynamics
150
TABLE 2 Example of path in loop space for D = 2
6
Transition
1. C ~ C ' ( 3 , 1 )
3
Resulting loop 1
it 3['i:1 so -
2
2. C - , C
I I
3. C ---,C-(6,2)
4. C --, C-(5,1)
13~12 5. C ---, C - ( 4 , - 2)
5-
6. C ~ C2.1OCLO,2
7. C ---, C
Resulting loop 2
-4
:1 !: !
8. C --* C - ( 6 , - 1)
9. C--.C
]i
I
-5
8
7
I
I
I I
II
I
I
TABLE 2
(continued)
i0. C ~ C
/
/
11. C---,C
I I 5-
1/ II -4
II
13. C ~ C
]
/
I1
14. C---*C
/
/
II
15. C --* C (2, - 1)
plaquette 2-
16. C --, C - ( 4 , 1 )
1 2tl0
18. C ---, C - ( 2 , - 11)
6
!]
3
2
~
4
20. C ---"C
I I
21. C - . C
I I
22. C ~ C
II
23. C ---, C
I I
2,:1. C - - + C
I (2,-2)
9-5
1-
19. C ~ C - ( 7 , 2 )
25. C - - , C
7
11
12"
17. C --, C - ( 6 , - 1)
48
I
plaqucttc
G. Marchesini / Loop (l~,narnics
152
TABLE 3 Rate of successful path
D= 2 D= 4
C2
C4
90% 20%
66% 1%
,X~3
Q2b
F=0.11 o
0.1
;:oo
~,:~:o
:.,1__~_~2°-.°°--~=~-o~ o = : ° _ . ~ : .
....... ¢
¢¢
o
,-;=~2==:
oo
oe
o
0
.-.
I
.
.
.
.
.
...................... ,NPath
I
.
10
20
3()
40
50x10 ~
(a) 2
J
,.,~, = 2
0.~--
!
7" ~0.08
O# 03i ° o
u°
o
• o
0.2~
°
.. o
o
•
•
u°
°
o oo ' : :___ .;=:
Q*
•
°
oO
.......
•
=
i
i
10
_.L
I
20
30
o
o
oo
o
. . . . .
o
o°
I
40
50x103
(b) Fig. 6. Monte Carlo results for the calculation of W(C) for U ( N ) gauge theory in D dimensions for N ~ oc. These results refer to D = 2 and C given by a two-plaquettc loop (n c =6). Solid lines axe the exact results. The dots and the circle~ refer to results obtained with two different sets of random number generators. (a) Strong-coupling region. (b) The phase transition point. (c) Weak-coupling region. The crossed dots are obtained from the dots by removing the fluctuation at :'Vpath----35 000. The crossed circles are obtained from the circles by removing the two large fluctuations at ,¥path -- 30000 and at :Vpath= 40000.
G. Marchesini / Loop dynamicv
153
O O O O
.~=1
O
r =0.04
O O t
O
O
O =O
O
O
O O O O O
O
° o
O
*ix I °°.[o**
i
°
0.5~
O m
Q D
I
Q
elx
X X ~ X ~ ( X ) < X -~ X . x
~ X X
°°
0 0 0 0 0 0
9
I0 0 0 0
B
i
O
I
•
L
...........
10
I
I
20
30
......
_*
i
0I
t
40
..... j 50x103
fig.6c Fig. 6 (continued)
suitable relaxation parameters). The rare paths with a large fluctuation for these numbers, instead, may give a very large contribution. The effect of these fluctuations on the convergence wilt be analyzed later. 4.3. F I R S T N U M E R I C A L
RESULTS
As mentioned, the aim of the present calculations is: (i) to obtain a first indication about the efficiency of our importance sampling by evaluating the rate of successful paths; (ii) to study the convergence of the formalism in D = 2 and to explore the relative importance of the fluctuations. The reported results are obtained for the following probabilities p_ = 0.3,
p+ = 0.1,
Ps = 0 . 3 ,
P0 = 0 . 3 ,
( w i t h intersections),
Ps = 0,
P0 = 0.6,
(without intersections),
(4.10)
In table 2 we report an example of a path C ~ P in D = 2, where C is a loop which wraps around two plaquettes ( n c = 6 ) . According to (4.8) this path give the
(71.Marchesini / Loop dynamics
154
contribution 2r )10(
2r ]z(
2r
1
0 ,4)
In order to study the rate of successful paths C ~ P we have placed a cutoff in the numbers of actual transitions. In the notation of (4.7) we set N + + N + Ns'+ N,~" ~< 100.
(4.12)
We considered two initial loops: C 2 the loop around two plaquettes as in table 2 and C 4 the square loop around four plaquettes. We report in table 3 the rate of successful paths. This result shows that even for this simple form of importance sampling, the rates of successful paths are relatively high also for D = 4. Our present program executes over 5000 transitions in one second of CPU time on a CDC 7600 computer with an average n c < 40. We come now to discuss the convergence of the algorithm and its fluctuations. We report only the results for the D = 2 theory, where the solution is known [13]. Using the exact value for W1, W(C) is computed and then compared with its exact value. In fig. 6 we report the results for the calculation of W(C), where C is the two-plaquette loop (n c = 6). This is obtained from 105 paths generated with the a priori probabilities in (4.10). The cutoff in (4.12) is assumed and 90% of these paths are successful, i.e. contribute to (4.9). We show the results for three typical values of the couping constant: ~ = 3 in fig. 6a and h = 1 in fig. 6c for the strongand weak-coupling region respectively; and ~ = 2, the phase transition point, in fig. 6b. The relaxation parameter r is fixed at a value for which the contribution of a path in (4.8) is w(path)_< 1 when the rates of the various transitions in (4.7) are around their average values. To get a feeling about the occurrence of large fluctuations we report the data for two runs with Np~th up to 5 . 1 0 4 (tWO different seeds were chosen for random number generator): dots refer to run 1 and circles to run 2. In fig. 6a (strong-coupling region) the average approaches rapidly the exact value which is reached by Npath -- 10 000. For the run 1 the average is stable while for the run 2 fluctuations appear at Np~th= 30000 and Neath = 40000. In fig. 6b (~. at the critical value) the situation is similar except that the fluctuations in the run 2 at Npath ----30 000 and Np~th = 40000 are more pronounced. The two fluctuations are with opposite sign so that they compensate for Npath > 40 000. In the weak-coupling region these fluctuations become very large (fig. 6c). For the run 2 they are not yet washed out at Npath = 50000. As mentioned before such fluctuations are due to paths for which the rates of the different loop transitions are very far from their average. For instance the path responsible for the first large fluctuations of the run 2 (Npath=30000) has 28 steps with 25 loops with no intersections. The rates (4.7) for this path are N0' = 2,
No'' = 23,
U~ = 1,
N_ = 2,
Ns = 0,
(4.13)
G. Marchesini / Loop dynamics"
155
to be constructed with the average values for instance of No" = 15 and N_ = 8.5. As a result the large factors in (4.8) due to the C ~ C steps are not compensated for. These fluctuations of indefinite sign, although very rare, may undermine the convergence of the algorithm. To emphasize this in fig. 6c we have removed the paths with [w(path)[ > 1000.
(4.14)
For the run 1 there is a single path at Npath----35 000 which gives a contribution w(path) = 2600. The cross-dotted points are obtained by removing this path. For the run 2 the effect is more relevant. There are two paths which exceed the bound (4.14): one around Npath ~--30000 with w(path)-- - 2 1 600 and the other around Npath-40000 with w(path)-- +8800. When these two paths are removed the cross-circles are obtained which are remarkably close to the exact result. The picture we have described is rather typical also for other runs we have made: the convergence is rather stable in the strong-coupling region. In the weak-coupling region rare but large fluctuations obscure the convergence. When they are removed, even with a very large cut off, the exact result is approached. We have also performed this calculation for the square loop around four plaquettes and confirmed these features. 5. Conclusions In this paper we have proposed and tested a new numerical algorithm for pure gauge theories. The algorithm aims at solving the loop equations of the theory. The advantages with respect to more conventional Monte Carlo methods are clear: (a) No large computer memory is needed; the method is based on the motion of one or few loops in the space-time lattice. (b) A single path in the loop space gives a typical contribution in (4.8) with a simple analytical dependence on the parameters such as the coupling constant. This means that a single Monte Carlo run gives the full coupling constant dependence of a given W(C). (c) Finally the formalism is based only on physical quantities. The crucial point of the method is the importance sampling strate~" to increase the efficiency i.e. the rate of successful paths. We have shown that the simple one we have used for a first attempt gives a sufficient rate even at D = 4 at least for small loops. This choice however can be certainly improved by taking into account more information on the structure of the loop to be varied and on the path. We have fully discussed only the case N = 1 and N ~ oo. Along the lines of appendix A it is possible to generalize the formalism to other finite N. Some indications for the generalization have been pointed out in sect. 3 and appendix A and will be described in a future paper. It is a pleasure to thank M. Ciafaloni, F. Duimio and G. Veneziano for various discussions. I am also grateful to B. Webber for some advice on computer matters. I
156
G. Marchesini / Loopdynami~v
owe special thanks to E. Onofri for various discussions, advice and help especially in the last part of this work.
Appendix A In this appendix the formal solution of the loop equations of sect. 3 is deduced from the stochastic quantization [11,9] formalism. This may provide a useful technical framework to generalize the solution to U ( N ) theories with N finite and to other gauge theories. Moreover operators such as K in the text have here an explicit representation as differential operators. For a general loop C defined in (2.1) consider the loop function (D = 4) (p(C) = 1 Tr U(C),
(A.1)
u ( c ) = u , , ( < ) . . , g,,(x,,).
(A.2)
where
U~(x) are the usual U ( N ) gauge fields associated to the link x, x + e~. Consider now a five-dimensional electric field E~(x) associated to the link x, x + e~. For U ( N ) theory one has a = 1,2 . . . . . N 2,
(T~T ~ = ½N),
[e;(x),U~(x)]= r%(x),
(A.3)
and zero for non-corresponding links. A stochastic Fokker-Planck operator for this problem is
H=-~2 ~_~ e-S(V)E~,(x)eS(°lE~,(x).
(A.4)
x,lx.a
When for the action S ( U ) one takes the expression in (2.11) one finds
n r ( c ) = ,,c,~ (c) + 2 g
s,,r (%),~ (%)
j
nC
+X E E'['~(c+(J,~'))-'~(c-(J,"))], J=l
v
where the notation was introduced in sect. 2.
(A.5)
G. Marchesini / Loopdynamics
157
Of course for the U(1) theory in (A.5) one has ~(Cij)~(Cii ) = 9~(C). For N v~ 1 instead one has to consider the action of H on two-loop functions, and so on. Suppose that two loops Ct and C 2 overlap on a single rink. One has then . 2s1~
U(,p(C~)~(C~)) = (H~(C~)),p(C~) + ~(CI)(U~(C2)) + - ~ ( C ~
u C~), (A.6)
where C~ to C2 is a loop obtained by joining C~ and C 2 a t the overlapping link and st2 = _+1 according to whether the overlapping is parallel or antiparallel. The generalization to more complex or simple situations is obvious. For N --+ o¢ the last term in (A.6) drops and the action of H on multi-loop functions factorizes. The expression (A.4) shows that H has a non-negative spectrum and that the ground state wave function is q)(0)-- I corresponding to no loop. Introducing the abstract notation (see sect. 3)
~(c)=
~(o)=
(A.7)
one finds for the right 14,0 R) and left (4'}1 ground state of H I~R) = 10),
(A.8)
(ff~/I = (0J es(c'),
whose scalar product gives the partition function
z = ( ~,o~i~,o ~) = f[DU ]e s`~'.
(A.9)
Here DU represents the integration over all U(N) gauge fields U~(x). Under the hypothesis that the ground state (A.8) is unique one has the following expression for the ground state projection l i m e -u~
IxP~)(~P~}
(A.IO)
which allows the following representation for the Wilson loop function:
W(c)-l f [DU]es(v'l Tr U(C) =l(OfeS(U)[C) Z
= lira ( O l e - " ~ l C ) ,r __+ ~ C
(A.11) "
G. Marchesini / Loop dynamics
158
Similarly, for the two-loop Wilson function, one has
1 =-2fiDUleS(" '
w(cl,c2)
Tr V(C ) I Tr U(Cz)
= lim (01e-n~lC1,C2), '/" "--~ O(2
(A.12)
where (UICa, C2> = ( 1 / N ) T r U(C1)(1/N)Tr U(C2). It is clear from (A.7) that if one drops the last term of order ( 1 / N z) the action of e -n" factorizes and one proves then the factorization [12] W(C,,C2) = W(C1) W(C:)(1 + O ( 1 / N 2 ) ) .
(A.13)
It is interesting to observe, as was shown in ref. [9], that the loop equation recalled in sect. 2 can be readily deduced from the representation (A.11). In fact if the limit z ~ ~ exists one has 0 = lim (01e-mill C> = H ( W ( C ) ) .
(m.14)
The limit r ~ ~ in (A.11) can actually be done (formally) and W(C) can be represented in terms of a resolvent. To do th~s consider the identity (01e-talC> -- <01(T)zqC>,
(A.15)
where T = e -~n and e = "r/L. We shall keep e fixed while taking the limit r ~ ~ . Since TI0 > = [0> it is convenient to project out of T the state 10>:
T = 10><01 + T'.
(A.16)
Using (A.11) and (A.15) this gives W ( C ) = lim (01(T)zqC> = ]~ <0I(T')Jl C> L---, oc
=
(0l(1
j=0
-
T')-IlC>.
(A.17)
Note that this representation holds for any U ( N ) theory. From (A.16) one has T ' = e-~X(1 - 10>(01).
(A.18)
The value of e is actually immaterial as a consequence of the loop equations. As shown in sect. 3 these equations allow us to obtain a different representation for the resolvent in (A.17).
G. Marchesini / Loop dynamics
159
Appendix B The partition function of the single plaquette model is
Z=fdVeS' ,
(B.1)
where U is a matrix of the U(N) group and dU is the group invariant measure. The Wilson loops are ~ =(1TrUn
) = 21J drU e s w ~ l~ T r
U" .
(B.2)
These functions of ~. can be computed for any N and for N -o vc they undergo [13] a phase transition at h = 2. For instance for n = 1 { 1/~, ~b1= 1 - ¼ ~ ,
h> 2 2~<2.
(B.3)
When the formalism of sect. 3 is developed in this case, the Wilson loops are expressed as matrix elements of a resolvent. In this appendix we want to show that, when a relaxation parameter is introduced, the expansion of the resolvent converges also in the weak-coupling region. From this
,~.=~ r =0.05
'6 - - - -
_/
j
loo
~0 3
f
200
Fig. 7. Wilson loop for the singlc-plaquette model. Values of ~b.([) in (B.5) as function of I. The exact results are also indicated.
G. Marchesini / Loop dynamics
160
resolvent one can deduce the following expansion of the Gauss-Seidel type n-- 1
T,~(i+a)= (1 - r)~P U)' ~-~n (
-
1
--
r
E (.,.(i)a,(,) \~k 't'n
k
_~f)
,
1
(B.4)
~p(0i+a) = 1,
where r is the relaxation parameter and for instance ~(i=1)= 6n,0"In fig. 7 we report the results for X = 0.1. Even for this very small value of the coupling constant, the average /
7 converges quite fast to the exact values.
References [1] J.B. Kogut, Rev. MOd. Phys. 55 (1983) 775; M. Creutz, L. Jacobs and C. Rebbi, Phys. Reports 95 (1983) 201 [2] K.G. Wilson, Phys. Rev. D10 (1974) 2445; Phys. Reports 23 (1975) 331 [3] G.F. Dc Angelis, D. De Falco and F. Guerra, Nuovo Cim. Lett. 19 (1977) 55; F. Guerra, R. Marra and G. Immirzi, Nuovo Cim. Lett. 23 (19771 237 [4] Yu.M. Makeenko and A.A. Migdal, Phys. Lett. 88B (1979) 135 [5] A.M. Polyakov, Nucl. Phys. B164 (1979) 171 [6] D. Foerster, Phys. Lett. 87B (1979) 83 [7] T. Eguchi, Phys. Lett. 87B (1979) 91 [8] B. Sakita, Phys. Rev. D21 (1980) 1067: A. Jevicki and B. Sakita, Nucl. Phys. B165 (1980) 511 [9] G. Marehesini, Nucl. Phys. B191 (1981) 214 [10] (i.E. Forsythe and R.A. Leibler, Math. Tables Comp. 4 (1950) 127; F.J. Oswald, Matrix inversion by Monte Carlo methods, in Mathematical methods for digital computers vol. I (1960); J. Kuti, Phys. Rev. Lett. 49 (1982) 183; B. Berg and D. Foerster, Phys. Lett. 106B (1981) 323 [11] G-. Parisi and Y. Wu, Sci. Sinica 24 (1981) 483 [12] G. 't Hooft, Nucl. Phys. 72 (1974)461; G. Veneziano, Phys. Lett. B52 (1974) 220; E. Witten, Nucl. Phys. B160 (1979) 57 [13] D. Gross and E. Witten, Phys. Rev. D21 (1980) 446; G. Paffuti and P. Rossi, Phys. Lett. 92B (1980) 321
LOOP DYNAMICS FOR GAUGE THEORIES: A NUMERICAL ALGORITHM
G. MARCHESINI Dipartimento di Fisica, Universiti~ di Parma, INFN, sezione di Milano, Italy
Received 31 October 1983
The loop equations defining the lattice U(N) gauge thco~ are recalled and a formal solution is presented (for N = 1 and N --* zc). Wilson's expectation function W(C) for a loop C is expressed as a matrix element of a resolvent in the space of all the loops between C and 0 (no loop). It is shown that such a solution provides a new numerical algorithm to compute physical quantities. This is based on a Ulam-von Ncumann type of method for computing inverse matrLx elements by introducing importance sampling for the paths in the space of the matrix indices, which in this case are the loops. As a result W(C) is obtained by summing over important paths in the loop space connecting C to 0. A Monte Carlo program is presented, for the N ---,~c case, where a ve~' simple form for the importance sampling is introduced so that the computer time for each step in the construction of the path is minimized.The rates for successful paths (i.e. path C --*0 within a given finite number of steps) are computed for D = 2 and D = 4. Both rates and computer time involved encourage us to attempt a large scale calculation. Here the numerical studies of the convergence and of the fluctuations are presented only for D = 2. Convergence is rather fast, but specially in the weak-couplingregion rare and large fluctuations appear thus suggestingthat a bctter tuning for the importance samphng is needed.
1. Introduction T h e present numerical algorithms [1] for c o m p u t i n g physical quantities in (lattice) gauge theories [2] are based on the Metropolis m e t h o d in which the gauge fields are treated as d y n a m i c a l variables a p p r o a c h i n g the e q u i l i b r i u m in a Markov process. W i t h i n this method one meets with difficulties of two kinds: (i) since the gauge fields are r e d u n d a n t variables one does not k n o w a priori whether two different field configurations, generated in the M o n t e Carlo process, are physically different or rather connected b y gauge rotation, a n d this may reduce the efficiency of the method. ( i i ) A more i m p o r t a n t difficulty from a practical point of view is that the l i m i t a t i o n in the n u m b e r of the d y n a m i c a l variables implies that the calculations can be d o n e in a space-time world extremely small. Because of these difficulties it seems desirable to find a n d test alternative numerical algorithms which could provide a different perspective a n d eventually overcome some of the difficulties. T h e physical quantities in (pure) gauge theories are the Wilson loops [2]. A n a t u r a l f o r m u l a t i o n of gauge theory is then in terms of loop variables rather then in terms of 135
136
G. Marchesini / Loop dynamics
gauge fields. Recently some progress [3-8] in this direction has been achieved. In particular for the U ( N ) gauge theory dynamical equations for the Wilson loops have been obtained. They correspond [9] to the Schwinger-Dyson equations of the theory. For the two extreme cases N = 1 and N ~ oc these equations involve single Wilson loops only. In this paper we obtain a formal solution for these loop equations and show that this solution provides a numerical algorithm to compute physical quantities. For the moment we discuss in detail only the N = 1 and N ~ ~ cases. The formulation is in the space of loops (which lie on the space-time lattice): for instance in QED a matrix element K c , c is defined which relates a given loop C to a next-neighbour loop C'. The Wilson expectation value W(C) for a loop C is expressed in terms of the matrix element of the resolvent of K between C on the zero loop configuration. In spite of such an exotic matrix inversion problem, numerical algorithms can be proposed. In this paper we explore a method based on the Neumann expansion, where one has to compute the matrix element of K" between the loop C and the zero loop. The intermediate indices of this power of K are of course loops. The matrix element of K" can then be viewed as given by the sum over all the paths in the loop space which connect the initial loop C to zero with steps defined by K. The number of paths in loop space is very large so that this sum can be practically done only by an importance sampling [10] method. Note that although the number of loops is large, with such types of algorithm one deals for each step with one or few loops so that even for a very large space-time world there is no problem for computer memory; in fact a loop is a variable in one dimension, the proper time. To study this method we have devised a Monte Carlo program for the U ( N ) gauge theory in D dimensions for N ~ oo. The most crucial point for the efficiency of the method is of course the choice of the a priori probabilities for the importance sampling. In fact one has to maximize the rate of successful paths i.e. of paths in the loop space which start from the original loop C and shrink to zero in a finite number of steps. For the moment we have made the most simple choice by defining a priori probabilities depending on the individual steps, so that the shrinking of a loop can be favoured only in a limited way. In this way one minimizes the computer time for each step in the development of a path, but one expects a low rate of successful paths especially for large loops (actually we will show on the basis of the formal solution and of the loop equations that a large-loop calculation can be subdivided into smaller parts). A first aim of our numerical analysis is to study these rates: the result is that even in D = 4 the rate of successful paths is relatively high at least for small loops. A second part of the numerical analysis is devoted to the study of the convergence of the algorithm both in the strong- and weak-coupling region. We shall limit ourselves here to report the results for D = 2. In the strong-coupling region the convergence is rather fast and the fluctuations are small. In the weak-coupling region instead there appear very rare (with a relative frequency of order 10 -~) but large fluctuations of indefinite sign which obscure the convergence (in fact if these
G. Marchesini / Loopdynamics
137
fluctuations are artificially removed a rather fast convergence is obtained). This situation is typical of an importance sampling strategy not perfectly tuned. In sect. 2 we set the notation apt to describe the loops in a D-dimensional space-time lattice and for any given loop we define the set of its next-neighbour loops. In terms of this notation we recall the form of the loop equations for the U ( N ) gauge theories and discuss the specific cases N -- 1 and N --o ~ , the planar theory. In sect. 3 we introduce the space of loops and present the formal solution of the loop equations for N = 1 and N ~ oo in terms of a resolvent in the loop space. We describe here how from the Neumann expansion of this resolvent a numerical algorithm can be constructed. In sect. 4 we describe the details of the specific algorithm we have constructed for the planar D-dimensional theory, and report the results on the analysis of efficiency and convergence. In sect. 5 some conclusions are drawn. In appendix A we obtain the solution of the loop equations by using the formalism of stochastic quantization [11, 9] which may be useful especially to study the extension of our method to the N-finite case. Finally in appendix B we discuss the single-plaquette planar model.
2. Loops and loop equations The loop equations set relations among different Wilson loop expectation functions. In this section we recall the form of these equations for the U ( N ) gauge theory and the simplifications which occur for N = 1 and N -o oo. Before this however, (i) we shall set the notation to describe loops in a D-dimensional space-time lattice and (ii) we shall define the set of next-neighbour loops of any given loop C.
2.1. NOTATION AND NEXT-NEIGHBOUR LOOPS
Consider an oriented closed curve (loop) C in the space time lattice
c:
(2.1)
are the lattice points in D dimension (the lattice unit is set to be one). Two successive points x k , X k + 1 belong to a link in the lattice
x = (nln2...no)
Xk - 1 = Xk + e,,
(2.2)
where e~, are unit vectors in the positive and negative directions:/~ = _+1, _+2 , . . . , _+D. Of course C has no spikes i.e. for any k xk+ 1 ~ Xk_ 1 .
(2.3)
G. Marchesini / Loop ~lynamics
138
XK Fig. 1. Transformation C --, C ' ( k , ~,) in eq. (2.4) for the link x k, Xk_ v
A first set of loops next-neighbour to C is obtained by considering the following two classes of elementary deformations of C: (a) C ' ( k , u) is a loop obtained from C when the two successive points of C, x k and Xk. ~ ( = Xk+ e~,) are substituted by the set of points (cf. fig. 1) X k ' X k + l =e, X k ,
X k Jr- elL,
x k -[- elz ~- e~,,
x k-'F e u,
Xk~
Xk! l,
(2.4) where e~ is in any direction orthogonal to e,, and 1 ~< k ~
xk+e~,
xk+e~+e,,,
xk~ t.
(2.5)
Note that such a deformation of C may generate spikes. In this case C-(k, ~,) is actually obtained by removing them. In the case in which C has intersections one introduces another set of nextneighbour loops obtained by splitting C. The following possibilities are considered. (a) Two points of C x j a n d x k ( j < k) have a link parallel intersection when Xj = X k ,
Xj+ 1 = Xk !
1"
(2.6)
An example is given in fig. 3 where the link is intersecting only twice. In this case the loop C can be split in two loops
% - (Xj.,...Xk), Ckj = (x x ...XjXk. 1 . ..x,,c).
(2.7)
(b) Two points of C xj and x~ ( j < k) have a link antiparallel intersection when (fig. 4) xj = Xk,
t
xj÷ 1 = x~ _ ~.
1
Fig. 2. Transformation C --, C - ( k , v) in eq. (2.5) for the link xk, xk, 1.
(2.8)
G. Marchesini / Loopdynamics
v
139
J
~kj
Fig. 3. Example of a splitting for a parallel intersection case: C ~ CjkCtj.
As shown in fig. 4 the loop C can be split in Cjk~-~ ( X j + l " " X k
-- 2)'
Ckj =-- (X1...Xjxk+l...XnC ).
(2.9)
These two new loops may have spikes which have to be removed. Intersection of sites, but not on links (cf. fig. 5), are not considered for the loop equations. The inverse operation of combining two loops to make a single loop is relevant only for N ~ 1 and finite and is not considered here (see appendix A).
_5
tI
),
i --° v
,k
i
Fig. 4. Example of a splitting for an antiparallel intcrscction case: C --* CjkCg ~.
Ii Fig. 5. Example of an intersection only on the site. This case does not contribute to the loop equations since the intersection is at 90 degrees.
140
G. Marchesini / Loop dynanlics
2.2. LOOP EQUATIONS
One is now in a position to recall the loop equations for the Wilson function 1
/
(2.10)
where U(C) is the usual U(N) gauge field associated to the loop C, i.e. the ordered group product of link variables around the loop. The expectation value in (2.10) is obtained for instance by assuming the Wilson form for the action S(U):
S ( U ) = ~ - E Tr(Up+ U~),
(2.11)
P
where the sum extends over all elementary squares or plaquettes of the lattice and Up is the usual plaquette field. For this form of the action the loop equations are nc
1
E
0 = H{ W(C)} = ncW(C) + X
k=l
w(c (k,.))} v
+ 2 E sjkW(Cjk,Ckj),
(2.12)
j
where W ( C I ' C2) -
N1
TrU(C1 ) ~1T r U ( C 2 ) ) •
(2.13)
In (2.12) the notation of the previous subsection is used. Here n c is the number of points of the loop C; the sum over u is restricted to the directions orthogonal to the link xk, Xk~ 1 (cf. eqs. (2.4) and (2.5)). Finally in the last term Sjk = + 1 for a link parallel intersection at the points xj and Xk; while Sgk = --1 for a link antiparallel intersection (of course sj~ = 0 otherwise). The boundary condition for these loop equations is that when C is shrunk to zero U(0) = 1 so that
W(0) = 1.
(2.14)
Eq. (2.12) relates the single-loop expectation function to the two-loop expectations in (2.13). However for N = 1 and N ~ ~ these loop equations involve only single-loop expectations. In the case of QED the gauge field is a phase and onc finds in (2.12)
w(Cjk, Ck ~) = W(C),
(2.15)
G Marchesini / Loop dynamics
141
so that the loop equations become rt C
(nc+2Esjk ,
j
W(C)+~ E ~
k=l
E'{W(C=(k,P))-W( C
(k,u))}=0.
u
(2.16) It is interesting to note that, if the loop is at fixed time, the coefficient of W(C) is just the electric energy of the loop. In fact ec = n c + 2 E Sjk = E ( q [ j
-- q[)2,
(2.17)
I
where the sum over l extends to all the links which are crossed by the loop C and qt ~ is the number of times the loop crosses the link I in one or the other direction. The interpretation of this fact is clear in the formalism of appendix A. In the case of the U(N) theory for N --* ~ the simplification of the loop equations arises because of the factorization property [12] W(C1,C2) = W(Ct)W(C2)(1 + O ( 1 / N 2 ) ) .
(2.1.8)
3. Wilson's loop and paths in loop space In this section we present the formal solution of the loop equations (both for N = 1 and N ~ oo) in terms of matrix elements of a resolvent in the space of loops. We then show how the Neumann expansion of this matrix element can be viewed as obtained by summing over paths in the space of loops, and illustrate the possibility of performing this sum by importance sampling methods [10]. Here we simply discuss the general features of the method and leave to sect. 4 the presentation of a specific algorithm. 3.1. F O R M A L
SOLUTION
To a given loop C in (2.1) we associate a ket IC) with the normalization (C'[C) = 6(C',C),
(3.1)
where 8(C',C)= 1 if C and C' coincide both in points and orientation, and 8 ( C ' , C ) = 0 otherwise. We denote by 10) the state corresponding to no loop; by [C1,C2) the state corresponding to the two loops C 1 and C2; and so on. In appendix A we show, within the formalism of stochastic quantization [11], that the Wilson loop expectation W(C) can be expressed as the matrix element W(C) = (01WIC),
(3.2)
G. Marchesini / Loop dynamics
142
and the expression for W is obtained. Here we simply give the form of W for N = 1 and N ~ oc and show that the loop equations are satisfied. Consider first the simpler case of QED where the loop equations are given in (2.16) with the boundary (2.14), and where, due to (2.15), only single-loop states need to be introduced. In this case W QE° is given by the resolvent of K Q ~ :
W QED= (1 - KQED) -1 ,
(3.3)
where K QED is defined as follows (for notation see sect. 2)
KQEDIO) = 0,
KQ~IC)= -
(3.4)
E~, IC)+x-~c F, E'[IC-(j,~))-IC*(j,~))]Cj
]
j'~l
(3.5)
v
In fact, due to (3.4), the boundary condition W(0) = (01(1 - K Q E D ) - l l 0 ) = 1,
(3.6)
is satisfied. Moreover, from the form (3.3) one has (01 wQLD(1 -- KQED)IC ) = (01C) ,
(3.7)
and because of (3.5) this is just the loop equation for W(C). Consider now the N--, ~ theory. In this case one needs to introduce also the many-loop states since the expectations W(Ct, C2) in (2.13) are involved. In appendix A it is shown that for a general U ( N ) theory
w(c) = (01 wl c), w ( c ~ , c : ) = (Ol WlCl,C,).
(3.8)
For N ~ ~ moreover the factorization property in (2.18) implies (01 w l c t , c 2 ) = (01WlC1)(Ol WLC2)(1 + O(1/N2)),
(3.9)
which sets the property of the action of W on a multi-loop state in terms of the action on single-loop states. Analogously to QED one requires W of the form W = (1 - X ) - l ,
(3.10)
143
G. Marchesini / Loop dynamics
where XlO ) = O, i, )>1 .
KIC) = - - - 2 ~ sjilCij,Cj~) + 1 - -
nc
hnc j~l =
j
,
(3.11) As in the previous case one show that with this form both the boundary conditions and the loop equations are satisfied. 3.2. Q E D
To obtain some familiarity with a resolvent in a loop space consider the Neumann expansion* of (3.3) (the index QED is here dropped)
w ( c ) = <01(1 - x ) - ' l c > oC
=
'Koc,,...Kc2c,Kclc,
(3.12)
n = 0 C1C2 ... C,,
where the sum is restricted to non-zero loops (C i ~=0, i = 1. . . . . n) and the matrix elements of K are given by (3.5): the non vanishing ones are
Kc, c = (C'IK[C) =
1 - ec/nc,
C'=C
Yr 1 / h n c ,
c,=
(3.13)
This matrix then connects C to itself and to the next-neighbour loops C - ( j , p) introduced in sect. 2. ec is defined in (2.17). According to (3.12) W(C) is obtained by summing over all possible ways in which the initial loop C is shrunk to zero via the elementary steps in (3.13); these will be called paths C -~ 0 in the loop space. This picture suggest a possible way for computing W(C) by using numerical algorithms of the type proposed by Ulam and von Neumann [10]. On the basis of (3.13) one defines an importance sample for the paths C --, 0 and W(C) is obtained as an expectation value corresponding to (3.12). The crucial point for the efficiency of such an algorithm is of course the definition of the importance sampling. In fact the matrix (3.13) connects a loop C to * This corresponds to the strong-coupling expansions. To study the weak-coupling region one needs a representation for the resolvent which allows the analytic continuation to small ~. This wSll be discussed in sect. 4, and not here, so as to leave the argument simple.
144
G. Marchesini / Loop dynamics
4(D - 1)n c modified loops, in D space-time dimensions, so that quickly the number of loops becomes considerably large, and one may expect that in a finite number of steps only a negligible fraction of paths leads to a zero loop. In order to obtain a contribution to W(C) one has then to favour those steps in (3.13) which lead to smaller loops, that means in general the transition C ~ C - ( j , u). Note moreover that the contribution for these transitions is positive and in general larger. This explains why "although the matrix elements in (3.13) have not a definite sign, the value of W(C) is positive for not-nested loops. To show that the transition C ~ C - ( j , u) leads to larger contributions it is convenient to isolate from K in (3.5) its diagonal part K D. Introducing then K ' --- (1 - KD)- I(K - KI)) one finds
W(C) = <01(1 - K ) l I C > = n~-ff:c ec @1(1 - K ' ) - l l C > '
(3.14)
where K~,c = T n c ,
1
C,=C±(j,v)"
(3.15)
£C' ~?/C '
Note then that (nc,/ec,) is in general larger for C' = C - ( j , u) than C ' - C ~(j, p). Anyhow if the loop C is large, this may not be enough to obtain a considerable fraction of paths C ~ 0. However the loop equations may be used in this case to overcome this difficulty by making shorter the minimum required path. In fact observe that any path C---, 0, before reaching the zero loop, must reach a loop around a plaquette (C = P). The path C ---, 0 may then be divided in the first part C ~ P and then P ---, 0. This last step leads to the expectation function for the single plaquette W1 = W (C = P), which can then be factorized out, and one needs only to study the paths C ---, P where the plaquette is reached only in the last step. This may be achieved formally by using the decomposition K = Kp + K, Kp ~ E *K[ P)(P[,
(3.16) (3.17)
P
where IP) is the state corresponding to the plaquette loop and the sum is over all oriented plaquettes. The matrix K acts as in (3.5) except that the subspace of plaquette loops is projected out: R IP) = 0.
(3.18)
G. Marchesini / Loop dynamics
145
By using (3.16) one finds in fact WIG) = ( 0 1 ( 1 - K ) - l f C ) = W l ~ * ( P l ( 1 - g ) - ~ l C ) .
(3.19)
P
Due to (3.18) the plaquette states will never appear as intermediate states in the Neumann expansion of (1 - K)-a. This means that one needs to consider now only the paths C --* P (for any P) which come to the single plaquettes only in the last step. This procedure can of course be iterated by projecting out of K also the subspace of states next-neighbour of the single plaquette and so on. In this way the calculation of a large Wilson loop expectation may be split into calculations involving smaller steps. 3.3. U(N) THEORY FOR N--+
The paths related to the Neumann expansion of W(C) for the U(N) theory are far more complicated than in the previous QED case. In fact here one has to take into account also disconnected loops. Due to the form of K in (3.11) starting from IC) after a few steps one may generate a state corresponding to two loops, say IC 1, C2). The subsequent evolution of this state is complicated by the possibility of an interaction of the two loops to make a single loop again (see appendix A eq. (A.6)). This interaction term however is of order 1/N 2, so that for N ~ oe the factorization (3.9) holds and the two loops move on independent paths. As a result the Neumann expansion of W(C) corresponds to the sum of all paths which branch as tree and each branch independently shrinks to zero. This can be seen also formally by singling out from K in (3.11) the splitting part:
K=K'+K~,
(3.20)
where KsIC ) -
(C']K'IC) =
2
~
nC j
Sii]Cji,C,j),
-T-1/)~nc,
(3.21)
C ' = C~-(j, v).
(3.22)
1Ks)-i(1 - K ' ) IIC ) .
(3.23)
In fact one can write
w(c) =
(01(1 - x ) -11 c )
= (0](1 - ( 1 - K ' ) Introducing the "propagator"
D(C', C) = (C'1(1 - K ' ) - 1]C),
(3.24)
G. Marchesini / Loop dynamics
146
one finds from (3.23)
Esj~W(Cj,)w(c,i)D(C,C),
w(c) = D(0,C)- 2E !
C' nC' j < i
(3.25)
where the sum extends over all possible splittings c ' --, c ; i % . The tree equations (3.25) is characteristic of the N ~ ~ limit: after the branch C ' ~ C)'iCi~ the paths C,'j ~ 0 and Cj; ~ 0 are independent. As in the QED case, when the single-plaquette function W1 = W (C = P) is known, the value of W(C) can be obtained considering only the paths C ~ P which go to a single plaquette state IP) only in the last step. In fact as in QED by projecting out the single-plaquette subspace in K one obtains
W(C) = W t ~ *(P[(1 - K)-IIC),
(3.26)
P
where the sum is over all oriented plaquettes and in K the subspace of plaquette-loops is projected out: K = KE~(1 - IP)(P I). Moreover from the factorization property in (3.8) one has Wl(1 - l ~ ) - l l C l , C 2 ) = WI(1 - K ) - l i c 1 ) w t ( 1 - K ) - I I C 2 ) .
(3.27)
These two equations give an alternative structure of the paths needed to compute W(C) if W~ is known: one has to consider all the paths which may branch and each branch ends when a plaquette state is reached. According to (3.26) and (3.27) each branch gives a factor W1 so that, in contrast to the QED case, W(C) is not proportional to W1. Note finally that the splitting,possibility, absent in QED, considerably increases the fraction of successful paths since smaller loops are generated from larger ones.
4. Numerical algorithm and first results In this section we present the Monte Carlo program we have constructed along the lines of sect. 3 for the calculation of Wilson's loop in the planar theory (N ---, ~ ) . To be able to explore the weak-coupling region we introduce a relaxation parameter. As explained earlier the crucial point for the method is the choice of the importance sampling, and we describe in detail the expected feature of our choice. We then report first the numerical results aimed (i) to study the rate of successful paths in D = 2 and D = 4; (ii) to study for D = 2 the convergence of the algorithm both in the strong- and weak-coupling region, together with the role of fluctuation.
G. Marchesini / Loopdynamics
147
4.1. THE ALGORITHM
As previously explained the expansion of the resolvent in (3.26) corresponds to the strong-coupling expansion and therefore, to study also the weak-coupling region, one needs on analytical continuation. For this reason we introduce a relaxation parameter r. This can be done by using the following expression for W(C) which is equivalent to (3.26):
W(C) = G }--'.*(P[ (1 - A )
11C) ,
(4.1)
P
where A, which replace the matrix K, is defined as
AI0) =0, A IC) = (1
-
r -
rK)l
C),
(4.2)
and for the factorization (3.21) WI(1 - A ) - I I C I , C 2 ) = WI(1 - A ) - I l C t ) W ~ ( 1 - A)--llCz).
(4.3)
The non-vanishing matrix elements of A are then
Ac,c =
{C'IA
IC> =
1 - r, r T ~n---7 ,
C' = C C' = C ! (
j, v),
(4.4)
and for the splitting
(Cji, CijlA IC) =
2r
- --sj,. r/c
(4.5)
The representation (4.1) holds for any value of r and by taking r = O(~) this may be used to extend the convergence of the Neumann expansion of W(C) to the weak-coupling region. In appendix B we test the convergence for the single-plaquette planar model. In spite of the fact that this model undergoes a phase transition [13] even at very weak coupling the expansion is quickly convergent when r = O(~). We come now to discuss the most important point for the Monte Carlo program: the a priori probabilities for the importance sampling. For these probabilities we made a very simple choice: they are defined for each step of the path and are largely independent of the structure of the loop to be transformed. Specifically we introduce p _+as the total probability for the loop transitions in the class C --, C -(j, v); p~ as the total probability for the loop splitting (of course if the loop has no intersection p~ = 0). Within each class the probabilities for the individual transitions are equally
G. Marchesini / loop dynamics
148
distributed. Finally the probability for no transformation, P0, is defined by the normalization
(4.6)
p+ + p _ + p.~ + p o = 1.
W(C) is now obtained as an expectation value for the product of the matrix element of A accordingly normalized. In table 1 in correspondence with the class of loop transitions we report the multiplicities of the class, the total probabilities, the related matrix elements and finally the corresponding factors to be used in the expectation values. Of course the last row is absent when C has no intersections, in this case, according to (4.6), P0 = 1 - p + - p _ . The procedure for computing W(C) is then the following: starting from the loop C and using the probabilities in table 1 one generates a new loop C'. The process is then repeated for C'. Eventually after a few steps the loop splits. In this case one continues the evolution for the two loops separately ( N - , m). The evolution of a branch ends when a plaquette loop is reached, and the full process (path) ends when all branches end. We will call successful a path which ends within a given finite number of steps. Suppose that one of such paths has the following rates for the related transformations: N . : C --> C ~-:(j, p), N/: C ~ Ci, C 0
(Sii = + 1 ) ,
N,":c-,cj,%
(,j,=-l),
No': C ~ C
(C with intersections),
No": C ~ C
(C without intersections).
(4.7)
According to table 1 the contribution to W(C) from such a path is TAm.~:1 Loop transitionsand dynamicalfactors Transitions
Multiplicities
Probability
Matrix element
Factors
C--,C C--*C=(j,: ,)
1 2(D-1)n c
Po p~
1-r T-t - ~,nc
(1 r)/p o _T_2(D-1)r ?tp r
Mc = ~ 1.5,1
P~
C --, Cj,Co
j
- 2rSji
nc
- 2rSJ~' Mc Ps nc
G. Marchesini / Loopdynamics
w(path) = W~%(
1-r
)N (1-r
1 -P+-P--Ps
149 NIT
1 -p+-p_
×
(4.8) Here Nr~ = 1 + N / + N~" is the total number of branches and each branch brings a factor W1 according to (4.1); and n c , Me, are the length and number of intersections of the loop C i that underwent the splitting. Finally the value of W(C) is given by
W(C) =
lim -~path -~ OC
1
• w(path).
Npath p a t h
(4.9)
According to (4.8) and (4.9) W(C) is computed when W1 is known. However even Wt can be computed with this method by using the loop equations. In fact they set a relation between W1 and the two-plaquette loop expectations. By computing these last functions with (4.8) and (4.9) one obtains a relation to determine W~. 4.2. DISCUSSION ON THE ALGORITHM (a) Owing to our choice of the a priori probabilities (i.e. satisfying the Markov property) the computer time for each step in the path is minimized, but the shrinking of a loop can be favoured only in a limited way: we require p _ > p + in order to inhibit the formation of very nested loops. Because of the expected low efficiency of this importance sampling we have used the formulation which requires path C ~ P rather than C ~ 0. However, in order to compute large Wilson loops, one probably must split the calculation into smaller steps by using the loop equations along the lines described in subsect. 3.2. (b) As clear from (4.8) the sign of w(path) is not definite, therefore in the sum (4.9) one expects cancellations. These cancellations are not important in the strongcoupling region where the paths are dominated by the ones in which at each step the transition is C ~ C (j, p). In the weak-couping region instead one expects cancellations since these paths are not so dominant. (c) Moreover, within this formulation, one also expects large fluctuations in the weak-coupling region although rare. In fact in this region one requires a small relaxation parameter so that the transition C ---, C contributes with a factor (1 - r)/po sizably larger than one. However when the numbers in (4.7) for the various transitions are around their average values the large factors coming for C--, C transitions are compensated by the factors coming from the other transitions (for
G. Marchesini / Loop dynamics
150
TABLE 2 Example of path in loop space for D = 2
6
Transition
1. C ~ C ' ( 3 , 1 )
3
Resulting loop 1
it 3['i:1 so -
2
2. C - , C
I I
3. C ---,C-(6,2)
4. C --, C-(5,1)
13~12 5. C ---, C - ( 4 , - 2)
5-
6. C ~ C2.1OCLO,2
7. C ---, C
Resulting loop 2
-4
:1 !: !
8. C --* C - ( 6 , - 1)
9. C--.C
]i
I
-5
8
7
I
I
I I
II
I
I
TABLE 2
(continued)
i0. C ~ C
/
/
11. C---,C
I I 5-
1/ II -4
II
13. C ~ C
]
/
I1
14. C---*C
/
/
II
15. C --* C (2, - 1)
plaquette 2-
16. C --, C - ( 4 , 1 )
1 2tl0
18. C ---, C - ( 2 , - 11)
6
!]
3
2
~
4
20. C ---"C
I I
21. C - . C
I I
22. C ~ C
II
23. C ---, C
I I
2,:1. C - - + C
I (2,-2)
9-5
1-
19. C ~ C - ( 7 , 2 )
25. C - - , C
7
11
12"
17. C --, C - ( 6 , - 1)
48
I
plaqucttc
G. Marchesini / Loop (l~,narnics
152
TABLE 3 Rate of successful path
D= 2 D= 4
C2
C4
90% 20%
66% 1%
,X~3
Q2b
F=0.11 o
0.1
;:oo
~,:~:o
:.,1__~_~2°-.°°--~=~-o~ o = : ° _ . ~ : .
....... ¢
¢¢
o
,-;=~2==:
oo
oe
o
0
.-.
I
.
.
.
.
.
...................... ,NPath
I
.
10
20
3()
40
50x10 ~
(a) 2
J
,.,~, = 2
0.~--
!
7" ~0.08
O# 03i ° o
u°
o
• o
0.2~
°
.. o
o
•
•
u°
°
o oo ' : :___ .;=:
Q*
•
°
oO
.......
•
=
i
i
10
_.L
I
20
30
o
o
oo
o
. . . . .
o
o°
I
40
50x103
(b) Fig. 6. Monte Carlo results for the calculation of W(C) for U ( N ) gauge theory in D dimensions for N ~ oc. These results refer to D = 2 and C given by a two-plaquettc loop (n c =6). Solid lines axe the exact results. The dots and the circle~ refer to results obtained with two different sets of random number generators. (a) Strong-coupling region. (b) The phase transition point. (c) Weak-coupling region. The crossed dots are obtained from the dots by removing the fluctuation at :'Vpath----35 000. The crossed circles are obtained from the circles by removing the two large fluctuations at ,¥path -- 30000 and at :Vpath= 40000.
G. Marchesini / Loop dynamicv
153
O O O O
.~=1
O
r =0.04
O O t
O
O
O =O
O
O
O O O O O
O
° o
O
*ix I °°.[o**
i
°
0.5~
O m
Q D
I
Q
elx
X X ~ X ~ ( X ) < X -~ X . x
~ X X
°°
0 0 0 0 0 0
9
I0 0 0 0
B
i
O
I
•
L
...........
10
I
I
20
30
......
_*
i
0I
t
40
..... j 50x103
fig.6c Fig. 6 (continued)
suitable relaxation parameters). The rare paths with a large fluctuation for these numbers, instead, may give a very large contribution. The effect of these fluctuations on the convergence wilt be analyzed later. 4.3. F I R S T N U M E R I C A L
RESULTS
As mentioned, the aim of the present calculations is: (i) to obtain a first indication about the efficiency of our importance sampling by evaluating the rate of successful paths; (ii) to study the convergence of the formalism in D = 2 and to explore the relative importance of the fluctuations. The reported results are obtained for the following probabilities p_ = 0.3,
p+ = 0.1,
Ps = 0 . 3 ,
P0 = 0 . 3 ,
( w i t h intersections),
Ps = 0,
P0 = 0.6,
(without intersections),
(4.10)
In table 2 we report an example of a path C ~ P in D = 2, where C is a loop which wraps around two plaquettes ( n c = 6 ) . According to (4.8) this path give the
(71.Marchesini / Loop dynamics
154
contribution 2r )10(
2r ]z(
2r
1
0 ,4)
In order to study the rate of successful paths C ~ P we have placed a cutoff in the numbers of actual transitions. In the notation of (4.7) we set N + + N + Ns'+ N,~" ~< 100.
(4.12)
We considered two initial loops: C 2 the loop around two plaquettes as in table 2 and C 4 the square loop around four plaquettes. We report in table 3 the rate of successful paths. This result shows that even for this simple form of importance sampling, the rates of successful paths are relatively high also for D = 4. Our present program executes over 5000 transitions in one second of CPU time on a CDC 7600 computer with an average n c < 40. We come now to discuss the convergence of the algorithm and its fluctuations. We report only the results for the D = 2 theory, where the solution is known [13]. Using the exact value for W1, W(C) is computed and then compared with its exact value. In fig. 6 we report the results for the calculation of W(C), where C is the two-plaquette loop (n c = 6). This is obtained from 105 paths generated with the a priori probabilities in (4.10). The cutoff in (4.12) is assumed and 90% of these paths are successful, i.e. contribute to (4.9). We show the results for three typical values of the couping constant: ~ = 3 in fig. 6a and h = 1 in fig. 6c for the strongand weak-coupling region respectively; and ~ = 2, the phase transition point, in fig. 6b. The relaxation parameter r is fixed at a value for which the contribution of a path in (4.8) is w(path)_< 1 when the rates of the various transitions in (4.7) are around their average values. To get a feeling about the occurrence of large fluctuations we report the data for two runs with Np~th up to 5 . 1 0 4 (tWO different seeds were chosen for random number generator): dots refer to run 1 and circles to run 2. In fig. 6a (strong-coupling region) the average approaches rapidly the exact value which is reached by Npath -- 10 000. For the run 1 the average is stable while for the run 2 fluctuations appear at Np~th= 30000 and Neath = 40000. In fig. 6b (~. at the critical value) the situation is similar except that the fluctuations in the run 2 at Npath ----30 000 and Np~th = 40000 are more pronounced. The two fluctuations are with opposite sign so that they compensate for Npath > 40 000. In the weak-coupling region these fluctuations become very large (fig. 6c). For the run 2 they are not yet washed out at Npath = 50000. As mentioned before such fluctuations are due to paths for which the rates of the different loop transitions are very far from their average. For instance the path responsible for the first large fluctuations of the run 2 (Npath=30000) has 28 steps with 25 loops with no intersections. The rates (4.7) for this path are N0' = 2,
No'' = 23,
U~ = 1,
N_ = 2,
Ns = 0,
(4.13)
G. Marchesini / Loop dynamics"
155
to be constructed with the average values for instance of No" = 15 and N_ = 8.5. As a result the large factors in (4.8) due to the C ~ C steps are not compensated for. These fluctuations of indefinite sign, although very rare, may undermine the convergence of the algorithm. To emphasize this in fig. 6c we have removed the paths with [w(path)[ > 1000.
(4.14)
For the run 1 there is a single path at Npath----35 000 which gives a contribution w(path) = 2600. The cross-dotted points are obtained by removing this path. For the run 2 the effect is more relevant. There are two paths which exceed the bound (4.14): one around Npath ~--30000 with w(path)-- - 2 1 600 and the other around Npath-40000 with w(path)-- +8800. When these two paths are removed the cross-circles are obtained which are remarkably close to the exact result. The picture we have described is rather typical also for other runs we have made: the convergence is rather stable in the strong-coupling region. In the weak-coupling region rare but large fluctuations obscure the convergence. When they are removed, even with a very large cut off, the exact result is approached. We have also performed this calculation for the square loop around four plaquettes and confirmed these features. 5. Conclusions In this paper we have proposed and tested a new numerical algorithm for pure gauge theories. The algorithm aims at solving the loop equations of the theory. The advantages with respect to more conventional Monte Carlo methods are clear: (a) No large computer memory is needed; the method is based on the motion of one or few loops in the space-time lattice. (b) A single path in the loop space gives a typical contribution in (4.8) with a simple analytical dependence on the parameters such as the coupling constant. This means that a single Monte Carlo run gives the full coupling constant dependence of a given W(C). (c) Finally the formalism is based only on physical quantities. The crucial point of the method is the importance sampling strate~" to increase the efficiency i.e. the rate of successful paths. We have shown that the simple one we have used for a first attempt gives a sufficient rate even at D = 4 at least for small loops. This choice however can be certainly improved by taking into account more information on the structure of the loop to be varied and on the path. We have fully discussed only the case N = 1 and N ~ oo. Along the lines of appendix A it is possible to generalize the formalism to other finite N. Some indications for the generalization have been pointed out in sect. 3 and appendix A and will be described in a future paper. It is a pleasure to thank M. Ciafaloni, F. Duimio and G. Veneziano for various discussions. I am also grateful to B. Webber for some advice on computer matters. I
156
G. Marchesini / Loopdynami~v
owe special thanks to E. Onofri for various discussions, advice and help especially in the last part of this work.
Appendix A In this appendix the formal solution of the loop equations of sect. 3 is deduced from the stochastic quantization [11,9] formalism. This may provide a useful technical framework to generalize the solution to U ( N ) theories with N finite and to other gauge theories. Moreover operators such as K in the text have here an explicit representation as differential operators. For a general loop C defined in (2.1) consider the loop function (D = 4) (p(C) = 1 Tr U(C),
(A.1)
u ( c ) = u , , ( < ) . . , g,,(x,,).
(A.2)
where
U~(x) are the usual U ( N ) gauge fields associated to the link x, x + e~. Consider now a five-dimensional electric field E~(x) associated to the link x, x + e~. For U ( N ) theory one has a = 1,2 . . . . . N 2,
(T~T ~ = ½N),
[e;(x),U~(x)]= r%(x),
(A.3)
and zero for non-corresponding links. A stochastic Fokker-Planck operator for this problem is
H=-~2 ~_~ e-S(V)E~,(x)eS(°lE~,(x).
(A.4)
x,lx.a
When for the action S ( U ) one takes the expression in (2.11) one finds
n r ( c ) = ,,c,~ (c) + 2 g
s,,r (%),~ (%)
j
nC
+X E E'['~(c+(J,~'))-'~(c-(J,"))], J=l
v
where the notation was introduced in sect. 2.
(A.5)
G. Marchesini / Loopdynamics
157
Of course for the U(1) theory in (A.5) one has ~(Cij)~(Cii ) = 9~(C). For N v~ 1 instead one has to consider the action of H on two-loop functions, and so on. Suppose that two loops Ct and C 2 overlap on a single rink. One has then . 2s1~
U(,p(C~)~(C~)) = (H~(C~)),p(C~) + ~(CI)(U~(C2)) + - ~ ( C ~
u C~), (A.6)
where C~ to C2 is a loop obtained by joining C~ and C 2 a t the overlapping link and st2 = _+1 according to whether the overlapping is parallel or antiparallel. The generalization to more complex or simple situations is obvious. For N --+ o¢ the last term in (A.6) drops and the action of H on multi-loop functions factorizes. The expression (A.4) shows that H has a non-negative spectrum and that the ground state wave function is q)(0)-- I corresponding to no loop. Introducing the abstract notation (see sect. 3)
~(c)=
~(o)=
(A.7)
one finds for the right 14,0 R) and left (4'}1 ground state of H I~R) = 10),
(A.8)
(ff~/I = (0J es(c'),
whose scalar product gives the partition function
z = ( ~,o~i~,o ~) = f[DU ]e s`~'.
(A.9)
Here DU represents the integration over all U(N) gauge fields U~(x). Under the hypothesis that the ground state (A.8) is unique one has the following expression for the ground state projection l i m e -u~
IxP~)(~P~}
(A.IO)
which allows the following representation for the Wilson loop function:
W(c)-l f [DU]es(v'l Tr U(C) =l(OfeS(U)[C) Z
= lira ( O l e - " ~ l C ) ,r __+ ~ C
(A.11) "
G. Marchesini / Loop dynamics
158
Similarly, for the two-loop Wilson function, one has
1 =-2fiDUleS(" '
w(cl,c2)
Tr V(C ) I Tr U(Cz)
= lim (01e-n~lC1,C2), '/" "--~ O(2
(A.12)
where (UICa, C2> = ( 1 / N ) T r U(C1)(1/N)Tr U(C2). It is clear from (A.7) that if one drops the last term of order ( 1 / N z) the action of e -n" factorizes and one proves then the factorization [12] W(C,,C2) = W(C1) W(C:)(1 + O ( 1 / N 2 ) ) .
(A.13)
It is interesting to observe, as was shown in ref. [9], that the loop equation recalled in sect. 2 can be readily deduced from the representation (A.11). In fact if the limit z ~ ~ exists one has 0 = lim (01e-mill C> = H ( W ( C ) ) .
(m.14)
The limit r ~ ~ in (A.11) can actually be done (formally) and W(C) can be represented in terms of a resolvent. To do th~s consider the identity (01e-talC> -- <01(T)zqC>,
(A.15)
where T = e -~n and e = "r/L. We shall keep e fixed while taking the limit r ~ ~ . Since TI0 > = [0> it is convenient to project out of T the state 10>:
T = 10><01 + T'.
(A.16)
Using (A.11) and (A.15) this gives W ( C ) = lim (01(T)zqC> = ]~ <0I(T')Jl C> L---, oc
=
(0l(1
j=0
-
T')-IlC>.
(A.17)
Note that this representation holds for any U ( N ) theory. From (A.16) one has T ' = e-~X(1 - 10>(01).
(A.18)
The value of e is actually immaterial as a consequence of the loop equations. As shown in sect. 3 these equations allow us to obtain a different representation for the resolvent in (A.17).
G. Marchesini / Loop dynamics
159
Appendix B The partition function of the single plaquette model is
Z=fdVeS' ,
(B.1)
where U is a matrix of the U(N) group and dU is the group invariant measure. The Wilson loops are ~ =(1TrUn
) = 21J drU e s w ~ l~ T r
U" .
(B.2)
These functions of ~. can be computed for any N and for N -o vc they undergo [13] a phase transition at h = 2. For instance for n = 1 { 1/~, ~b1= 1 - ¼ ~ ,
h> 2 2~<2.
(B.3)
When the formalism of sect. 3 is developed in this case, the Wilson loops are expressed as matrix elements of a resolvent. In this appendix we want to show that, when a relaxation parameter is introduced, the expansion of the resolvent converges also in the weak-coupling region. From this
,~.=~ r =0.05
'6 - - - -
_/
j
loo
~0 3
f
200
Fig. 7. Wilson loop for the singlc-plaquette model. Values of ~b.([) in (B.5) as function of I. The exact results are also indicated.
G. Marchesini / Loop dynamics
160
resolvent one can deduce the following expansion of the Gauss-Seidel type n-- 1
T,~(i+a)= (1 - r)~P U)' ~-~n (
-
1
--
r
E (.,.(i)a,(,) \~k 't'n
k
_~f)
,
1
(B.4)
~p(0i+a) = 1,
where r is the relaxation parameter and for instance ~(i=1)= 6n,0"In fig. 7 we report the results for X = 0.1. Even for this very small value of the coupling constant, the average /
7 converges quite fast to the exact values.
References [1] J.B. Kogut, Rev. MOd. Phys. 55 (1983) 775; M. Creutz, L. Jacobs and C. Rebbi, Phys. Reports 95 (1983) 201 [2] K.G. Wilson, Phys. Rev. D10 (1974) 2445; Phys. Reports 23 (1975) 331 [3] G.F. Dc Angelis, D. De Falco and F. Guerra, Nuovo Cim. Lett. 19 (1977) 55; F. Guerra, R. Marra and G. Immirzi, Nuovo Cim. Lett. 23 (19771 237 [4] Yu.M. Makeenko and A.A. Migdal, Phys. Lett. 88B (1979) 135 [5] A.M. Polyakov, Nucl. Phys. B164 (1979) 171 [6] D. Foerster, Phys. Lett. 87B (1979) 83 [7] T. Eguchi, Phys. Lett. 87B (1979) 91 [8] B. Sakita, Phys. Rev. D21 (1980) 1067: A. Jevicki and B. Sakita, Nucl. Phys. B165 (1980) 511 [9] G. Marehesini, Nucl. Phys. B191 (1981) 214 [10] (i.E. Forsythe and R.A. Leibler, Math. Tables Comp. 4 (1950) 127; F.J. Oswald, Matrix inversion by Monte Carlo methods, in Mathematical methods for digital computers vol. I (1960); J. Kuti, Phys. Rev. Lett. 49 (1982) 183; B. Berg and D. Foerster, Phys. Lett. 106B (1981) 323 [11] G-. Parisi and Y. Wu, Sci. Sinica 24 (1981) 483 [12] G. 't Hooft, Nucl. Phys. 72 (1974)461; G. Veneziano, Phys. Lett. B52 (1974) 220; E. Witten, Nucl. Phys. B160 (1979) 57 [13] D. Gross and E. Witten, Phys. Rev. D21 (1980) 446; G. Paffuti and P. Rossi, Phys. Lett. 92B (1980) 321