- Email: [email protected]
Numerical simulations of quantum chromodynamics
Volume 124B, number 1,2
PHYSICS LETTERS
21 April 1983
NUMERICAL SIMULATIONS OF QUANTUM CHROMODYNAMICS * H.W. HAMBER The Institute for Advanced Study', Princeton, NJ 08540, USA E. MARINARI CEN Saclay, Gif-Sur Yvette, France G. PARISI Universitd di Roma [1, Tor Vergata, ltaly and C. REBB1 Brookhaven National Laboratory, Upton, N Y 11973, USA
Received 25 January 1983
A previously proposed method to perform Monte Carlo simulations with dynamical fermions is applied to lattice QCD. The main features of the computational technique are reviewed and results for some selected quantities, such as the average plaquette action and the quark condensate (~-t~) are presented. Both the dependence on the quark mass and on the number of flavors are considered. An estimate of the scale parameter in the presence of fermions is given.
Monte Carlo (MC) simulations of systems with fermions pose special problems due to the anticommuting nature o f the fermionic fields. In ref. [1] a method was proposed to simulate a fermionic system by means of a parallel set of bosonic variables, the socalled pseudofermions, and in ref. [2] the method was successfully applied to a numerical analysis o f the Schwinger model. In this letter we report the first results o f the application of the method to a realistic SU(3) model o f interacting quarks and gluons. The possibility, opened by our technique, of numerically simulating the complete theory of quantum chromodynamics makes the interest of even preliminary results obvious. We recall the most relevant features of the lattice formulation o f a fermionic system and of the m e t h o d Supported by the US Department of Energy under contract No. DE-AC02-76CH00016.
proposed in ref. [l ]. We shall adopt Susskind's discretization of Dirac's lagrangian, in the representation of ref. [3 ]. Computational convenience and the advantage of a surviving continuous chiral symmetry for massless fermions [4,5] motivate our choice. A single component ~x of the fermionic field is associated to each lattice site x and the fermionic action is 1 /2 /2! SF = ~ ~rx(~xU x~x+~ - ~ x + ~ U ~ x ) X,/2
+m ~
~x~x.
(1)
X
In this equation we have set the lattice spacing a equal to 1 (as we shall do in most o f the following), UUx is the gauge variable associated with the link emerging from x in the positive/l direction,/St is a displacement vector of 1 lattice unit in t h e / l direction, and PUx is a factor equal to 1 for/a = 1, to ( - 1 ) xl f o r / l = 2, to 99
( - 1 ) xz+x2 for/~ = 3, to ( - 1 ) xl+x2+x3 for #L= 4. In Susskind's formulation these factors effectively play the role of 7 matrices, inducing an analogous algebra in Fourier space. We shall write SF=~
XX t
~x(Dxx,+mOxx,)Ox,=~(D+m)$ ,
(2)
the operator D(U) being implicitly defined by eq. (1). Color indices will be suppressed throughout most of this paper. S F, together with the standard Wilson action for the gauge field S G = 13 ~
(1
--
~-RetrUV?x
21 April 1983
PHYSICS LETTERS
Volume 124B, number 1,2
" ut -v
- u x")
U x+OU x+~U
X,#
under the "3,5'' transformation ~x -" (-1) xl+x2+x3+x4 X ~x guarantees positivity of det(D+m) and det(D+m): det(D~ + m). We shall make this explicit by using a measure factor where the square root of the product [det(D t + rn)] [det(D + m)] = det[(Dt + m)(D + m)] appears. Furthermore, the action of eq. (1) is known to describe four independent flavors in the low-frequency limit. If one introduced n identical fields ¢(i)x, i = 1 .... n, interacting only through their common coupling to U, the number of flavors would be nf = 4n and the determinant would be raised to the power n = ] n f in the measure. Thus we shall use as measure factor exp(-SG)(det[(D? + m)(D
def. = 13 ~
E# v ,
(3)
+ rn)]}nf/8 ,
or, equivalently
X,,~
forms the total action S = S G + S F. The vacuum expectation value of a quantum observable 0 (U, ~-, 4) is given by
(O)=z-l f [-[ dU# [1 (d*xd~x) Oexp(-S), (4) X,Ia
X
with
Z=fx~,. dUff ~x (dff/xd~x) exp(-S).
(5)
Common to several proposals [1,6,7] to perform MC simulations of fermionic systems is a preliminary step carrying out the gaussian integration over the fermionic variables. This reduces the computation of the vacuum expectation values to an average over the bosonic gauge variables only, but with a measure given by exp [-SG(U) ] det[D(U) + m]:
(0> =z - l f I-I dUUx(O)Uexp(-SG)det(D+m), (6) X,l~
Z
= f H dU"x exp(-SG) det(D + m ) ,
(7)
X,#
where (O)u represents the expectation value of O (as determined by the gaussian integration) in the presence of a fixed background gauge field Ux~. For our purposes it is convenient to rewrite the measure in eqs. (6) and (7) in a slightly different form. First we notice that the property D + m ~ D t + m = -/9 + m 100
exp [-Serf(U)]
(8)
with Seff(U)= SG(U)
-~nfTrln[(Dt+m)(D+m)] . '
(9)
In this last equation Tr stands for the full trace of the operator, whereas traces over color indices only are denoted by tr. Although the measure factor of eqs. (8) and (9) defines a fully consistent quantum theory only for nf multiple of 4, we will use it also with nf = 2, trying to compensate in this fashion for the unwanted multiplication of flavors introduced by the Susskind formulation. The difficulty in performing numerical simulations derives from the non-local nature of Seff, which makes the evaluation of the change in action induced by the upgrading U-~ U' computationally very involved. The "quenched approximation", which has been successfully used in the analysis of the quark model spectrum [8-14] avoids the problem altogether by setting Serf(U ) = S G (U). The method proposed in ref. [1] assumes first that the variation of U in the upgrading procedure, 6U = U' - U, is taken sufficiently small so that terms of order higher than the first in 6U can be safely ignored in the evaluation of 6Seff. 6Seff is thus approximated by 6Seff = 5SG - ~nf ~ (tr J u ; S U U x + h.c.) , x,~.
where the "current" J"x is given by
(10)
Volume 124B, number 1,2
PHYSICS LETTERS Table 1 Results for ~ = 5.5, nf = 2.
-1
J~xee' = ~I'Ux [(D + m)x+;zcxe' - (D + + m)x+lOcxe,]
(11)
(color indices are here made explicit). Then the method offers a recipe to calculate JUx efficiently. The basic idea is that the current induced by a background gauge field is the same whether the propagating fields are bosonic or fermionic (provided o f course they are coupled through the same lagrangian); it is only the sign by which the term with JUx enters in 6Sef f which distinguishes between the two cases. JUx can therefore be computed b y a MC simulation over a set o f bosonic variables. We introduce bosonic fields q~x, ~bx (the pseudofermions) and, for convenience, also fields Xx, Y,x deffmed by X = (D + m) ¢ .
(12)
We denote by double brackets averages taken over the pseudofermions with measure exp [-~-(D ) + m)(D + m) ~b] = exp (-XX) • From the fundamental relation -1
[(Dt + m)(D + m)]xc,x'c' = <(~x'c'e&c})
(13)
(14)
it is straightforward to prove the equation ~ - <(~xc'Xx+~>>) JUxce, =1~ r~u (((~e,~x+~>>
21 April 1983
(15)
Thus the matrix elements o f the current can be expressed as averages o f suitable produces o f the and X fields. These are estimated b y Monte Carlo simulation. The algorithm becomes computationally convenient if all the required J~x can be evaluated before upgrading the gauge variables. Again, this is justified if 8 U is taken sufficiently small. Thus, the whole procedure, spelled out in detail, works as follows. N Monte Carlo steps are performed over the pseudofermionic variables ¢ and ~-; the currents are computed by averaging over these N iterations; then one Monte Carlo step is performed over all the gauge variables U, using the previously estimated currents to calculate the variations 8Sef f. In our computations we have considered an 84 lattice with periodic boundary conditions. We have performed 300 MC iterations on the pure gauge system from an ordered start at ~ = 5.5, obtaining a configuration which we denote b y C1, and then 500 fur-
Conf.
m
N
(~-ff)
(E}
(J- U)
C1 50 100 50 100 C2 50 100 50 100
0.2 0.2 0.2 0.2 0.2 0.1 0.1 0.1 0.1 0.1
quenched 20 20 40 40 quenched 20 20 40 40
1.091 1,033 1,006 1.062 1.062 0.943 0.798 0.770 0.853 0.776
0.5030 0.4904 0.4792 0.4949 0.4908 0.5030 0.4808 0.4797 0.4888 0.4814
0.6964 0.6975 0.6978 0.6966 0.6981 0.7270 0.7273 0.7300 0.7284 0.7272
ther iterations at/3 = 5.8 obtaining a configuration C 2. C1 and C2 have been the starting configurations for sequences of MC iterations including now also the dynamical reactions from the fermions. Precisely, we have run 100 MC iterations at/3 = 5.5 starting from C1, with m = 0.2, nf = 2 , N = 20 and 4 0 , 1 0 0 MC iterations at fl = 5.8 starting from C2, with m = 0.2, nf = 2 , N = 20 and 40, and nf = 4 , N = 20. Similar runs with m = 0.1 were also performed, but starting from the 50th configuration obtained in the corresponding runs with m = 0.2. The 50th and 100th configurations of all runs were recorded on tape for further analysis. A few more technical details are discussed later. Our results are summarized in tables 1 - 3 . Displayed are the average value o f ~ f (defined as (1/Nsites) × Y'x ~x fix), the average value o f the plaquette internal energy E (defined as (1/Nplaquettes) X ~-,x,u
Table 2 Results for/3 = 5.8, nf = 2. Conf.
m
N
(~-4' >
(E>
(J- U}
C2 50 100 50 100 C2 50 100 50 100
0.2 0.2 0.2 0.2 0.2 0.1 0.1 0.I 0.1 0.1
quenched 20 20 40 40 quenched 20 20 40 40
0.830 0.769 0.762 0.784 0.788 0.583 0.498 0.517 0.522 0.505
0.4328 0.4188 0.4176 0.4180 0.4220 0.4328 0.4172 0.4184 0.4214 0.4202
0.7086 0,7094 0.7113 0.7116 0.7106 0.7358 0.7380 0.7380 0.7358 0.7381
101
Volume 124B, number 1,2
PHYSICS LETTERS
Table 3 Results for t3 = 5.8,nf = 4. Conf.
m
N
(~)
(E)
(J.U)
C2 50 100 C2 50 100
0.2 0.2 0.2 0.1 0.1 0.1
quenched 20 20 quenched 20 20
0.830 0.748 0.750 0.583 0.470 0.461
0.4328 0.4107 0.4095 0.4328 0.4099 0.4092
0.7986 0.7102 0.7135 0.7358 0.7380 0.7390
between J~x and UUx (defined as (1]Nlinks) X Y~x,u Re tr JUtxx UU). The quantities in tables 1 - 3 represent exact averages over the configurations considered, with numerical errors o f the order of the last significant digit. But one must remember that all observables fluctuate from configuration to configuration. We have limited our analysis to a few sets of configurations only because o f computational constraints. F r o m less precise evaluations of the observables in the course o f the simulations, however, we estimate the fluctuations o f @-~) to be of the order 1.5% and those o f (E) to be of the order of 0.5%. Also in comparing the values of (E) and ( J - U) it is useful to remember that 6Seff, in the upgrade of the gauge variables UUx, is given b y the sum of the variations of the internal energies o f 6 placjuettes, multiplied by/3, • 1 minus ~nf times Re tr JU)c6U~x. Thus the variations of the quantities E and J . U are weighted quite differently in the variation o f S e f f. The results show a clear decrease of @-~) and (E) when the dynamical reaction from the fermions is taken into account. As one might have expected the inclusion of the fermionic determinant in the measure shifts the statistical equilibrium toward more ordered configurations, lowering the internal energy and the
21 April 1983
average value o f ff ft. The effect becomes larger the lower the mass and the greater the number of flavors, again in agreement with theoretical expectations. While the variations in b o t h ( f ~ ) and (E) induced by the inclusion of the fermionic reaction is larger than the fluctuations from configuration to configuration, the differences among results obtained with N equal to 20 or 40 and with 50 or 100 MC iterations are compatible with statistical fluctuations and no clear systematic effect is discernible. If anything, the decrease in ( ~ ) and (E) seems slightly larger with N = 20 than with N = 40, although the difference becomes less marked as the number o f MC iterations increases. The average value o f J " U appears stable and not very sensitive to the parameters of the computation. Averaging the results obtained with N = 20 and 40 and 50 and 100 MC iterations we find the values reported in table 4. Extrapolating the results for @-~) to m = 0 we find (f~)Q=0.795,
(~)2f=0.557
(~)Q
(~-~)2f = 0 . 2 4 4 ,
= 0.336,
at/3=5.5,
(~-~)4f = 0.181
at/3 = 5 . 8 .
(17)
From these numbers we may attempt an estimate of the variation of the lattice scale parameter induced b y the fermionic dynamical reaction although the results, based as they are on a very limited sample of configurations, m a y be affected b y serious statistical errors. The relation between a, A and/3 is given by the formula a = A 1 ([87r2/(33 _ 2nf)]
fl)"459-57nf)/(33-2nf)2
× exp [--47r2/3/(33 --2nf)] [1 +O(/3-1)] .
Table 4 Summary of tile results for (~-qJ) and (E).
102
(16)
~3
nf
m
(~)
(~qJ)quenched
A(~-~ ) (%)
(E)
(/:')quenched
A(E) (%)
5.5 5.5 5.8 5.8 5.8 5.8
2 2 2 2 4 4
0.2 0.1 0.2 0.1 02 0.1
1.041 0.799 0.776 0.510 0.749 0.465
1.091 0.943 0.830 0.583 0.830 0.583
-4.6 -15.3 -6.5 -12.5 -9.8 -20.2
0.4888 0.4826 0.4191 0.4193 0.4101 0.4095
0.5030 0.5030 0.4328 0.4328 0.4328 0.4328
-2.8 -4.1 -3.2 -3.1 -5.2 -5.4
(18)
PHYSICS LETTERS
Volume 124B, number ! ,2
@-~) should scale as a 3 times an anomalous dimension factor, which to the order it has been computed, does not depend on nf. Thus the ratio between the quantities (~-~) computed with or without fermions should be just given by the third power of the corresponding ratios of lattice spacings. From eqs. (16), (17) and (18) we then obtain Azf/A Q = 0.465 A2f/AQ=0.460,
A4f/AQ=0.145
at/3 = 5.5,
(19)
at/3---5.8.
(20)
A decrease in A with the inclusion of inner fermionic loops is expected, and the consistency between the results at/3 = 5.5 and/3 = 5.8 is gratifying. However, the ratios A2f/AQ and A4f/AQ are much smaller than one would expect from a perturbative analysis of the relationship between Alattic e and a more conventional scale such as Affg [15], which would give A2f/AQ = 0 . 6 5 7 ,
A4f/AQ = 0 . 3 7 7 .
(21)
The discrepancy may mean that the 0(/3-1) terms neglected in eq. (18) are still important at the values of/3 considered; that the predictions from the quenched approximation and from the complete theory for the physical value of the condensate (~-~) are different (or equivalently that the two schemes lead to different estimates for A ~ ) ; that our method fails to account completely for the fermionic determinant in the measure, either because of systematic effects or because the simulation has not yet converged to equilibrium [much smaller values for ( ~ ) w o u l d be required for consistency with the ratios in eq. (21)]; or, of course, any combination of the above. The convergence toward equilibrium has been tested by doing longer runs on smaller lattices (64) on the VAX computers at BNL and the University of Rome. In some instances we have observed further, but contained, decreases in ( ~ ) and (k'). Also, some dependence has been noted on the width of the step used in the bosonic upgrading. Details will be given in a separate publication. In any event, further theoretical and analytical work will be needed to discriminate among the various possibilities cited above. Finally, we present a few more details about our computation together with considerations on the approximations involved. Our procedure becomes exact in the limits where: (i) 6U, in the upgrading, tends to zero;
21 April 1983
(ii) N tends to infinity. About point (i), the upgrading of the gauge field was done following the Metropolis procedure and taking for U' the product of U times a matrix R chosen at random among a set of 200 SU(3) matrices R k stored in memory. Apart from the requirement that R k and R~-1 should belong to the set, the matricesR k were constructed exponentiating the generators i~ k • k, where the components of k are Gell-Mann's matrices and the vectors ~k were chosen stochastically in the neighborhood of the origin with a distribution approximately gaussian and (~2) = 9/300. Eight upgradings of the same link were attempted before proceeding to the next one: the average number of accepted changes turned out to be on the order of 5. Altogether the change 6 U in the upgrading was generally contained. We have performed tests of the dependence of the results on the average magnitude of 6 U and will report about our conclusions in a separate paper, About point (ii), the upgrading of the pseudofermions was also done with the Metropolis algorithm and three attempted upgradings per variable. Actually, given the gaussian form of the pseudofermionic measure, a heat bath procedure is more efficient. We have written a code that does the upgrading with the heat bath method, but, since we had already done a substantial number of MC iterations with our former program, for consistency we have used the Metropolis algorithm throughout the runs. We have performed tests which showed that our upgrading procedure gives results consistent with the heat bath algorithm. Insofar as the dependence on N is concerned, 20 or 40 MC steps are not sufficient to achieve a good convergence of any specific J~x to its correct value: the fluctuations may still be of the order of 30%. However, on the average, the value of the currents is well approximated. A simulation done on a definite configuration produced the following sequence of results for (J- U): ( ] . U) = 0.7423 after 1 MC step, = 0.7352 after 5 steps, = 0.7313 after 10 steps, = 0.7351 after 20 steps, = 0.7361 after 40 steps, versus a correct value (J • U) = 0.7358. It can be shown that errors in the evaluation of 6S have only a small effect on the results on a Monte Carlo simulation, provided they are randomly distributed and average to zero [16]. Thus it is likely that good approximate results may be achieved also with N relatively small. This is supported by the consistency of the data with N = 20 and N = 40. 103
Volume 124B, number 1,2
PHYSICS LETTERS
The simulations we have reported about have been performed with the CDC 7600 computer at BNL. Our code requires approximately 18 + 3N s per MC iteration. Thus, even with N = 40, the computational requirement for a simulation including fermionic dynamical effects is only a few times greater than for the corresponding simulation of the gauge system alone. Beyond these computations, extensive simulations on a 64 lattice have been done on the VAX computers at the University of Rome and BNL and on the CDC 7600 computer at Saclay. We plan to present the results of these simulations and of a more ref'med analysis of the data, together with further theoretical considerations in a separate publication. One of us (HWH) wishes to thank the Data Acquisition Group of CBA at BNL for the use of their VAX 11/780. The work of HWH has been supported by US Department of Energy under grant No. DE-AC0276ER02220.
References [1 ] F. Fucito, E. Marinari, G. Parisi and C. Rebbi, Nucl. Phys. B180 (1981) 360.
104
21 April 1983
[2] E. Marinari, G. Parisi and C. Rebbi, Nucl. Phys. B190 (1981) 734. [3] N. Kawamoto and J. Smit, Nucl. Phys. B192 (1981) 10. [4] C. Rebbi, Brookhaven preprint, to be published in the Proceedings of the 19th Orbis Scientiae Meeting (1982). [5 ] H. Kluberg-Stern, A. Morel, O. Napoly and B. Peterson, Saclay preprint (1982). [6] D.J. Scalapino and R.L. Sugar, Phys. Rev. Lett. 46 (1981) 519. [7] J. Kuti, Phys. Rev. Lett. 49 (1982) 183. [8] H. Hamber and G. Parisi, Phys. Rev. Lett. 47 (1981) 1795;Phys. Rev. D (1982), to be published. [9] E. Marinari, G. Parisi and C. Rebbi, Phys. Rev. Lett. 47 (1981) 1798. [10] H. Hamber, E. Marinari, G. Parisi and C. Rebbi, Phys. Lett.108B (1982) 314. [11] D. Weingarten, Phys. kett. B109 (1982) 57; 1UHET preprint (1982). [12] F. Fucito, G. MartineUi,C. Omero, G. Parisi, R. Petronzio and F. Rapuano, CERN preprint TH-3288 (1982). [13] H. Hamber, Institute for Advanced Studies preprint (October 1982). [14] C. Bernard, T. Draper and K. Olynyk, UCLA preprint (1982). [15] H. Sharatchandra, H. Tun and P. Weisz, Nucl. Phys. B192 (1981) 205. [16] C. Lang, private communication.
PHYSICS LETTERS
21 April 1983
NUMERICAL SIMULATIONS OF QUANTUM CHROMODYNAMICS * H.W. HAMBER The Institute for Advanced Study', Princeton, NJ 08540, USA E. MARINARI CEN Saclay, Gif-Sur Yvette, France G. PARISI Universitd di Roma [1, Tor Vergata, ltaly and C. REBB1 Brookhaven National Laboratory, Upton, N Y 11973, USA
Received 25 January 1983
A previously proposed method to perform Monte Carlo simulations with dynamical fermions is applied to lattice QCD. The main features of the computational technique are reviewed and results for some selected quantities, such as the average plaquette action and the quark condensate (~-t~) are presented. Both the dependence on the quark mass and on the number of flavors are considered. An estimate of the scale parameter in the presence of fermions is given.
Monte Carlo (MC) simulations of systems with fermions pose special problems due to the anticommuting nature o f the fermionic fields. In ref. [1] a method was proposed to simulate a fermionic system by means of a parallel set of bosonic variables, the socalled pseudofermions, and in ref. [2] the method was successfully applied to a numerical analysis o f the Schwinger model. In this letter we report the first results o f the application of the method to a realistic SU(3) model o f interacting quarks and gluons. The possibility, opened by our technique, of numerically simulating the complete theory of quantum chromodynamics makes the interest of even preliminary results obvious. We recall the most relevant features of the lattice formulation o f a fermionic system and of the m e t h o d Supported by the US Department of Energy under contract No. DE-AC02-76CH00016.
proposed in ref. [l ]. We shall adopt Susskind's discretization of Dirac's lagrangian, in the representation of ref. [3 ]. Computational convenience and the advantage of a surviving continuous chiral symmetry for massless fermions [4,5] motivate our choice. A single component ~x of the fermionic field is associated to each lattice site x and the fermionic action is 1 /2 /2! SF = ~ ~rx(~xU x~x+~ - ~ x + ~ U ~ x ) X,/2
+m ~
~x~x.
(1)
X
In this equation we have set the lattice spacing a equal to 1 (as we shall do in most o f the following), UUx is the gauge variable associated with the link emerging from x in the positive/l direction,/St is a displacement vector of 1 lattice unit in t h e / l direction, and PUx is a factor equal to 1 for/a = 1, to ( - 1 ) xl f o r / l = 2, to 99
( - 1 ) xz+x2 for/~ = 3, to ( - 1 ) xl+x2+x3 for #L= 4. In Susskind's formulation these factors effectively play the role of 7 matrices, inducing an analogous algebra in Fourier space. We shall write SF=~
XX t
~x(Dxx,+mOxx,)Ox,=~(D+m)$ ,
(2)
the operator D(U) being implicitly defined by eq. (1). Color indices will be suppressed throughout most of this paper. S F, together with the standard Wilson action for the gauge field S G = 13 ~
(1
--
~-RetrUV?x
21 April 1983
PHYSICS LETTERS
Volume 124B, number 1,2
" ut -v
- u x")
U x+OU x+~U
X,#
under the "3,5'' transformation ~x -" (-1) xl+x2+x3+x4 X ~x guarantees positivity of det(D+m) and det(D+m): det(D~ + m). We shall make this explicit by using a measure factor where the square root of the product [det(D t + rn)] [det(D + m)] = det[(Dt + m)(D + m)] appears. Furthermore, the action of eq. (1) is known to describe four independent flavors in the low-frequency limit. If one introduced n identical fields ¢(i)x, i = 1 .... n, interacting only through their common coupling to U, the number of flavors would be nf = 4n and the determinant would be raised to the power n = ] n f in the measure. Thus we shall use as measure factor exp(-SG)(det[(D? + m)(D
def. = 13 ~
E# v ,
(3)
+ rn)]}nf/8 ,
or, equivalently
X,,~
forms the total action S = S G + S F. The vacuum expectation value of a quantum observable 0 (U, ~-, 4) is given by
(O)=z-l f [-[ dU# [1 (d*xd~x) Oexp(-S), (4) X,Ia
X
with
Z=fx~,. dUff ~x (dff/xd~x) exp(-S).
(5)
Common to several proposals [1,6,7] to perform MC simulations of fermionic systems is a preliminary step carrying out the gaussian integration over the fermionic variables. This reduces the computation of the vacuum expectation values to an average over the bosonic gauge variables only, but with a measure given by exp [-SG(U) ] det[D(U) + m]:
(0> =z - l f I-I dUUx(O)Uexp(-SG)det(D+m), (6) X,l~
Z
= f H dU"x exp(-SG) det(D + m ) ,
(7)
X,#
where (O)u represents the expectation value of O (as determined by the gaussian integration) in the presence of a fixed background gauge field Ux~. For our purposes it is convenient to rewrite the measure in eqs. (6) and (7) in a slightly different form. First we notice that the property D + m ~ D t + m = -/9 + m 100
exp [-Serf(U)]
(8)
with Seff(U)= SG(U)
-~nfTrln[(Dt+m)(D+m)] . '
(9)
In this last equation Tr stands for the full trace of the operator, whereas traces over color indices only are denoted by tr. Although the measure factor of eqs. (8) and (9) defines a fully consistent quantum theory only for nf multiple of 4, we will use it also with nf = 2, trying to compensate in this fashion for the unwanted multiplication of flavors introduced by the Susskind formulation. The difficulty in performing numerical simulations derives from the non-local nature of Seff, which makes the evaluation of the change in action induced by the upgrading U-~ U' computationally very involved. The "quenched approximation", which has been successfully used in the analysis of the quark model spectrum [8-14] avoids the problem altogether by setting Serf(U ) = S G (U). The method proposed in ref. [1] assumes first that the variation of U in the upgrading procedure, 6U = U' - U, is taken sufficiently small so that terms of order higher than the first in 6U can be safely ignored in the evaluation of 6Seff. 6Seff is thus approximated by 6Seff = 5SG - ~nf ~ (tr J u ; S U U x + h.c.) , x,~.
where the "current" J"x is given by
(10)
Volume 124B, number 1,2
PHYSICS LETTERS Table 1 Results for ~ = 5.5, nf = 2.
-1
J~xee' = ~I'Ux [(D + m)x+;zcxe' - (D + + m)x+lOcxe,]
(11)
(color indices are here made explicit). Then the method offers a recipe to calculate JUx efficiently. The basic idea is that the current induced by a background gauge field is the same whether the propagating fields are bosonic or fermionic (provided o f course they are coupled through the same lagrangian); it is only the sign by which the term with JUx enters in 6Sef f which distinguishes between the two cases. JUx can therefore be computed b y a MC simulation over a set o f bosonic variables. We introduce bosonic fields q~x, ~bx (the pseudofermions) and, for convenience, also fields Xx, Y,x deffmed by X = (D + m) ¢ .
(12)
We denote by double brackets averages taken over the pseudofermions with measure exp [-~-(D ) + m)(D + m) ~b] = exp (-XX) • From the fundamental relation -1
[(Dt + m)(D + m)]xc,x'c' = <(~x'c'e&c})
(13)
(14)
it is straightforward to prove the equation ~ - <(~xc'Xx+~>>) JUxce, =1~ r~u (((~e,~x+~>>
21 April 1983
(15)
Thus the matrix elements o f the current can be expressed as averages o f suitable produces o f the and X fields. These are estimated b y Monte Carlo simulation. The algorithm becomes computationally convenient if all the required J~x can be evaluated before upgrading the gauge variables. Again, this is justified if 8 U is taken sufficiently small. Thus, the whole procedure, spelled out in detail, works as follows. N Monte Carlo steps are performed over the pseudofermionic variables ¢ and ~-; the currents are computed by averaging over these N iterations; then one Monte Carlo step is performed over all the gauge variables U, using the previously estimated currents to calculate the variations 8Sef f. In our computations we have considered an 84 lattice with periodic boundary conditions. We have performed 300 MC iterations on the pure gauge system from an ordered start at ~ = 5.5, obtaining a configuration which we denote b y C1, and then 500 fur-
Conf.
m
N
(~-ff)
(E}
(J- U)
C1 50 100 50 100 C2 50 100 50 100
0.2 0.2 0.2 0.2 0.2 0.1 0.1 0.1 0.1 0.1
quenched 20 20 40 40 quenched 20 20 40 40
1.091 1,033 1,006 1.062 1.062 0.943 0.798 0.770 0.853 0.776
0.5030 0.4904 0.4792 0.4949 0.4908 0.5030 0.4808 0.4797 0.4888 0.4814
0.6964 0.6975 0.6978 0.6966 0.6981 0.7270 0.7273 0.7300 0.7284 0.7272
ther iterations at/3 = 5.8 obtaining a configuration C 2. C1 and C2 have been the starting configurations for sequences of MC iterations including now also the dynamical reactions from the fermions. Precisely, we have run 100 MC iterations at/3 = 5.5 starting from C1, with m = 0.2, nf = 2 , N = 20 and 4 0 , 1 0 0 MC iterations at fl = 5.8 starting from C2, with m = 0.2, nf = 2 , N = 20 and 40, and nf = 4 , N = 20. Similar runs with m = 0.1 were also performed, but starting from the 50th configuration obtained in the corresponding runs with m = 0.2. The 50th and 100th configurations of all runs were recorded on tape for further analysis. A few more technical details are discussed later. Our results are summarized in tables 1 - 3 . Displayed are the average value o f ~ f (defined as (1/Nsites) × Y'x ~x fix), the average value o f the plaquette internal energy E (defined as (1/Nplaquettes) X ~-,x,u
Table 2 Results for/3 = 5.8, nf = 2. Conf.
m
N
(~-4' >
(E>
(J- U}
C2 50 100 50 100 C2 50 100 50 100
0.2 0.2 0.2 0.2 0.2 0.1 0.1 0.I 0.1 0.1
quenched 20 20 40 40 quenched 20 20 40 40
0.830 0.769 0.762 0.784 0.788 0.583 0.498 0.517 0.522 0.505
0.4328 0.4188 0.4176 0.4180 0.4220 0.4328 0.4172 0.4184 0.4214 0.4202
0.7086 0,7094 0.7113 0.7116 0.7106 0.7358 0.7380 0.7380 0.7358 0.7381
101
Volume 124B, number 1,2
PHYSICS LETTERS
Table 3 Results for t3 = 5.8,nf = 4. Conf.
m
N
(~)
(E)
(J.U)
C2 50 100 C2 50 100
0.2 0.2 0.2 0.1 0.1 0.1
quenched 20 20 quenched 20 20
0.830 0.748 0.750 0.583 0.470 0.461
0.4328 0.4107 0.4095 0.4328 0.4099 0.4092
0.7986 0.7102 0.7135 0.7358 0.7380 0.7390
between J~x and UUx (defined as (1]Nlinks) X Y~x,u Re tr JUtxx UU). The quantities in tables 1 - 3 represent exact averages over the configurations considered, with numerical errors o f the order of the last significant digit. But one must remember that all observables fluctuate from configuration to configuration. We have limited our analysis to a few sets of configurations only because o f computational constraints. F r o m less precise evaluations of the observables in the course o f the simulations, however, we estimate the fluctuations o f @-~) to be of the order 1.5% and those o f (E) to be of the order of 0.5%. Also in comparing the values of (E) and ( J - U) it is useful to remember that 6Seff, in the upgrade of the gauge variables UUx, is given b y the sum of the variations of the internal energies o f 6 placjuettes, multiplied by/3, • 1 minus ~nf times Re tr JU)c6U~x. Thus the variations of the quantities E and J . U are weighted quite differently in the variation o f S e f f. The results show a clear decrease of @-~) and (E) when the dynamical reaction from the fermions is taken into account. As one might have expected the inclusion of the fermionic determinant in the measure shifts the statistical equilibrium toward more ordered configurations, lowering the internal energy and the
21 April 1983
average value o f ff ft. The effect becomes larger the lower the mass and the greater the number of flavors, again in agreement with theoretical expectations. While the variations in b o t h ( f ~ ) and (E) induced by the inclusion of the fermionic reaction is larger than the fluctuations from configuration to configuration, the differences among results obtained with N equal to 20 or 40 and with 50 or 100 MC iterations are compatible with statistical fluctuations and no clear systematic effect is discernible. If anything, the decrease in ( ~ ) and (E) seems slightly larger with N = 20 than with N = 40, although the difference becomes less marked as the number o f MC iterations increases. The average value o f J " U appears stable and not very sensitive to the parameters of the computation. Averaging the results obtained with N = 20 and 40 and 50 and 100 MC iterations we find the values reported in table 4. Extrapolating the results for @-~) to m = 0 we find (f~)Q=0.795,
(~)2f=0.557
(~)Q
(~-~)2f = 0 . 2 4 4 ,
= 0.336,
at/3=5.5,
(~-~)4f = 0.181
at/3 = 5 . 8 .
(17)
From these numbers we may attempt an estimate of the variation of the lattice scale parameter induced b y the fermionic dynamical reaction although the results, based as they are on a very limited sample of configurations, m a y be affected b y serious statistical errors. The relation between a, A and/3 is given by the formula a = A 1 ([87r2/(33 _ 2nf)]
fl)"459-57nf)/(33-2nf)2
× exp [--47r2/3/(33 --2nf)] [1 +O(/3-1)] .
Table 4 Summary of tile results for (~-qJ) and (E).
102
(16)
~3
nf
m
(~)
(~qJ)quenched
A(~-~ ) (%)
(E)
(/:')quenched
A(E) (%)
5.5 5.5 5.8 5.8 5.8 5.8
2 2 2 2 4 4
0.2 0.1 0.2 0.1 02 0.1
1.041 0.799 0.776 0.510 0.749 0.465
1.091 0.943 0.830 0.583 0.830 0.583
-4.6 -15.3 -6.5 -12.5 -9.8 -20.2
0.4888 0.4826 0.4191 0.4193 0.4101 0.4095
0.5030 0.5030 0.4328 0.4328 0.4328 0.4328
-2.8 -4.1 -3.2 -3.1 -5.2 -5.4
(18)
PHYSICS LETTERS
Volume 124B, number ! ,2
@-~) should scale as a 3 times an anomalous dimension factor, which to the order it has been computed, does not depend on nf. Thus the ratio between the quantities (~-~) computed with or without fermions should be just given by the third power of the corresponding ratios of lattice spacings. From eqs. (16), (17) and (18) we then obtain Azf/A Q = 0.465 A2f/AQ=0.460,
A4f/AQ=0.145
at/3 = 5.5,
(19)
at/3---5.8.
(20)
A decrease in A with the inclusion of inner fermionic loops is expected, and the consistency between the results at/3 = 5.5 and/3 = 5.8 is gratifying. However, the ratios A2f/AQ and A4f/AQ are much smaller than one would expect from a perturbative analysis of the relationship between Alattic e and a more conventional scale such as Affg [15], which would give A2f/AQ = 0 . 6 5 7 ,
A4f/AQ = 0 . 3 7 7 .
(21)
The discrepancy may mean that the 0(/3-1) terms neglected in eq. (18) are still important at the values of/3 considered; that the predictions from the quenched approximation and from the complete theory for the physical value of the condensate (~-~) are different (or equivalently that the two schemes lead to different estimates for A ~ ) ; that our method fails to account completely for the fermionic determinant in the measure, either because of systematic effects or because the simulation has not yet converged to equilibrium [much smaller values for ( ~ ) w o u l d be required for consistency with the ratios in eq. (21)]; or, of course, any combination of the above. The convergence toward equilibrium has been tested by doing longer runs on smaller lattices (64) on the VAX computers at BNL and the University of Rome. In some instances we have observed further, but contained, decreases in ( ~ ) and (k'). Also, some dependence has been noted on the width of the step used in the bosonic upgrading. Details will be given in a separate publication. In any event, further theoretical and analytical work will be needed to discriminate among the various possibilities cited above. Finally, we present a few more details about our computation together with considerations on the approximations involved. Our procedure becomes exact in the limits where: (i) 6U, in the upgrading, tends to zero;
21 April 1983
(ii) N tends to infinity. About point (i), the upgrading of the gauge field was done following the Metropolis procedure and taking for U' the product of U times a matrix R chosen at random among a set of 200 SU(3) matrices R k stored in memory. Apart from the requirement that R k and R~-1 should belong to the set, the matricesR k were constructed exponentiating the generators i~ k • k, where the components of k are Gell-Mann's matrices and the vectors ~k were chosen stochastically in the neighborhood of the origin with a distribution approximately gaussian and (~2) = 9/300. Eight upgradings of the same link were attempted before proceeding to the next one: the average number of accepted changes turned out to be on the order of 5. Altogether the change 6 U in the upgrading was generally contained. We have performed tests of the dependence of the results on the average magnitude of 6 U and will report about our conclusions in a separate paper, About point (ii), the upgrading of the pseudofermions was also done with the Metropolis algorithm and three attempted upgradings per variable. Actually, given the gaussian form of the pseudofermionic measure, a heat bath procedure is more efficient. We have written a code that does the upgrading with the heat bath method, but, since we had already done a substantial number of MC iterations with our former program, for consistency we have used the Metropolis algorithm throughout the runs. We have performed tests which showed that our upgrading procedure gives results consistent with the heat bath algorithm. Insofar as the dependence on N is concerned, 20 or 40 MC steps are not sufficient to achieve a good convergence of any specific J~x to its correct value: the fluctuations may still be of the order of 30%. However, on the average, the value of the currents is well approximated. A simulation done on a definite configuration produced the following sequence of results for (J- U): ( ] . U) = 0.7423 after 1 MC step, = 0.7352 after 5 steps, = 0.7313 after 10 steps, = 0.7351 after 20 steps, = 0.7361 after 40 steps, versus a correct value (J • U) = 0.7358. It can be shown that errors in the evaluation of 6S have only a small effect on the results on a Monte Carlo simulation, provided they are randomly distributed and average to zero [16]. Thus it is likely that good approximate results may be achieved also with N relatively small. This is supported by the consistency of the data with N = 20 and N = 40. 103
Volume 124B, number 1,2
PHYSICS LETTERS
The simulations we have reported about have been performed with the CDC 7600 computer at BNL. Our code requires approximately 18 + 3N s per MC iteration. Thus, even with N = 40, the computational requirement for a simulation including fermionic dynamical effects is only a few times greater than for the corresponding simulation of the gauge system alone. Beyond these computations, extensive simulations on a 64 lattice have been done on the VAX computers at the University of Rome and BNL and on the CDC 7600 computer at Saclay. We plan to present the results of these simulations and of a more ref'med analysis of the data, together with further theoretical considerations in a separate publication. One of us (HWH) wishes to thank the Data Acquisition Group of CBA at BNL for the use of their VAX 11/780. The work of HWH has been supported by US Department of Energy under grant No. DE-AC0276ER02220.
References [1 ] F. Fucito, E. Marinari, G. Parisi and C. Rebbi, Nucl. Phys. B180 (1981) 360.
104
21 April 1983
[2] E. Marinari, G. Parisi and C. Rebbi, Nucl. Phys. B190 (1981) 734. [3] N. Kawamoto and J. Smit, Nucl. Phys. B192 (1981) 10. [4] C. Rebbi, Brookhaven preprint, to be published in the Proceedings of the 19th Orbis Scientiae Meeting (1982). [5 ] H. Kluberg-Stern, A. Morel, O. Napoly and B. Peterson, Saclay preprint (1982). [6] D.J. Scalapino and R.L. Sugar, Phys. Rev. Lett. 46 (1981) 519. [7] J. Kuti, Phys. Rev. Lett. 49 (1982) 183. [8] H. Hamber and G. Parisi, Phys. Rev. Lett. 47 (1981) 1795;Phys. Rev. D (1982), to be published. [9] E. Marinari, G. Parisi and C. Rebbi, Phys. Rev. Lett. 47 (1981) 1798. [10] H. Hamber, E. Marinari, G. Parisi and C. Rebbi, Phys. Lett.108B (1982) 314. [11] D. Weingarten, Phys. kett. B109 (1982) 57; 1UHET preprint (1982). [12] F. Fucito, G. MartineUi,C. Omero, G. Parisi, R. Petronzio and F. Rapuano, CERN preprint TH-3288 (1982). [13] H. Hamber, Institute for Advanced Studies preprint (October 1982). [14] C. Bernard, T. Draper and K. Olynyk, UCLA preprint (1982). [15] H. Sharatchandra, H. Tun and P. Weisz, Nucl. Phys. B192 (1981) 205. [16] C. Lang, private communication.