Nature of the thermal phase transition with Wilson quarks

We describe a series of simulations of high temperature QCD with two flavors of Wilson quarks aimed at clarifying the nature of the high temperature phase found in current simulations. Most of our work is with four time slices, although we include some runs with six and eight time slices for comparison. In addition to the usual thermodynamic observables, we study the quark mass defined by the divergence of the axial vector current and the quark propagator in the Landau gauge. We find that the sharpness of the Nt = 4 thermal transition has a maximum around K=0.19 and 6/g2= 4.8.

PHYSICAL REVIEW D VOLUME 49, NUMBER 7 1 APRIL 1994 Nature of the thermal phase transition with Wilson quarks Claude Bernard Department of Physics, Washington University, St. Louis, Missouri 63130 Thomas A. DeGrand and A. Hasenfratz Physics Department, University of Colorado, Boulder, Colorado 80309 Carleton DeTar Physics Department, University of Utah, Salt Lake City, Utah 84112 Steven Gottlieb Department of Physics, Indiana University, Bloomington, Indiana 47405 and Department of Physics, Bldg. 510A, Brookhaven National Laboratory, Upton, New York 11973 Leo Karkkainen and D. Toussaint Department of Physics, University of Arizona, Tucson, Arizona 85721 R. L. Sugar Department of Physics, University of California, Santa Barbara, California 93106 (Received 27 October 1993) We describe a series of simulations of high temperature QCD with two flavors of Wilson quarks aimed at clarifying the nature of the high temperature phase found in current simulations. Most of our work is with four time slices, although we include some runs with six and eight time slices for comparison. In addition to the usual thermodynamic observables, we study the quark mass defined by the divergence of the axial vector current and the quark propagator in the Landau gauge. We find that the sharpness of the N, = 4 thermal transition has a maximum around k=0. 19 and 6/g2=4.8. PACS number(s): 12.38.Gc, 11.15.Ha I. INTRODUCTION Lattice simulations are an important source of infor­mation on the behavior of quantum chromodynamics at high temperature. Most work has been done with Kogut-Susskind quarks because of the exact remnant of chiral symmetry. Since the exact chiral symmetry of Kogut-Susskind quarks is a U(l) symmetry, there is some question about how well the results reproduce the real world with its SU(2) chiral symmetry. In the continuum limit the complete chiral symmetry is restored. However, in the continuum limit the results should be independent of the regularization used for the quarks. To test this it is important to study high temperature QCD with the other common form of lattice quarks, the Wilson quarks. The first simulations of high temperature QCD with two flavors of Wilson quarks revealed a potential problem-for the values of 6/g2 for which most low temperature simulations were done, 4.5 <6/g25 5.7, the high temperature transition occurs at a value of quark hopping parameter k for which the pion mass measured at zero temperature is quite large [1,2]. In other words, it is difficult to find a set of parameters for which the tem­perature is the critical temperature and the quark mass is small. Further work confirmed that the pion mass is large at the deconfinement transition for this range of 6/g2 [3,4], (A recent study has concluded that for four time slices the chiral limit is reached at a very small value of 6/g2ss3.9 [5].) Screening masses for color singlet sources show an ap­proach to parity doubling in the high temperature phase similar to what is seen with Kogut-Susskind quarks [2,3]. Also, measurements of the pion mass show a shallow minimum at the high temperature transition [6], Previous simulations with Wilson fermions have locat­ed Kt, the value of the hopping parameter at which the high temperature crossover or phase transition occurs, as a function of 6/g2 for N, =4 and 6. The critical value of the hopping parameter, kc, for which the pion mass van­ishes at zero temperature has been located with some­what less precision [1,2,6,7,3,4], Some measurements of hadron masses have been carried out on zero temperature lattices for values of k and 6/g2 close to the k, curve, al­lowing one to set a scale for the temperature, and to esti­mate kc in the vicinity of the thermal transition [3,4,8]. In more recent work at N,=z6 we have observed coex­istence of the low and high temperature phases over long simulation times, and we have extended these observa­tions in the present project. The change in the plaquette across the transition is much larger than for the high temperature transition with Kogut-Susskind quarks [4]. This unexplained behavior, as well as work by Hasenfratz and DeGrand on the effect of heavy quarks [9], has led us to extend our work. 0556-2821/94/49(7)/3574( 15)/$06.00 49 3574 © 1994 The American Physical Society 49 NATURE OF THE THERMAL PHASE TRANSITION WITH . 3575 This paper reports on a series of simulations with Wil­son quarks at high temperature, in which we have studied a number of indicators for the nature of the phases. Us­ing 83 X 4 lattices, we have extended earlier studies of the location of the thermal transition or crossover to #f=0.20, 0.21, and 0.22. In the range of 6/g2 and k that have been studied earlier, we have done extensive simula­tions on 8X8X20X4 lattices, with additional work on 123 X 6, 122 X 24 X 6, and 82X20X8 lattices. For one value of 6/g2 we made a series of runs on 6X6X20X4 lattices to make sure that the effects we see are not due to the spatial size of the lattice. In this article we will con­centrate on the results with four time slices. Simulations with Nt= 6 and 8 are still underway and will be described later. We see a number of inexplicable effects. At large P and small k the crossover from the confined phase to the high temperature phase is smooth. Beginning at (j8,/c) = (5. 1,0.16) and extending down to about (/?,«■) = (4.51,0.20) the crossover becomes abrupt, though probably not first order. A rapid crossover is seen in the plaquette, real part of the Polyakov loop, the entro­py, and the quark mass derived from the axial vector current. For f}<4. 5,k>0.20 the transition once again becomes very smooth. Section II discusses the quantities we measured, and Sec. Ill summarizes the simulations and the results. Con­clusions are in Sec. IV. II. MEASURED QUANTITIES In our simulations we have measured the expectation values of the Polyakov loop, the space-space and space­time plaquettes, the chiral condensate rf/ip, the entropy, screening masses for meson sources, the quark mass defined by the divergence of the axial vector current, and quark propagators in Landau gauge. The expectation value of the Polyakov loop, (P), is simply interpreted as exp(- Fq/T), where Fq is the free energy of a static test quark. With dynamical quarks (P) is always nonzero, but it increases dramatically at the high temperature transition. We also measured the space-space and space-time plaquettes (DiS) and {). In our normalization these are equal to three on a com­pletely ordered lattice. The energy, pressure, entropy, and ipip with Wilson quarks are obtained by differentiating the partition func­tion with respect to the temporal size, the spatial size, and the quark mass, respectively. Details are given in the Appendix. We study the entropy to lowest order in g and tpip, using the formulas .4N 1 <S*a4>aiV, 3 day 3 daz = _J 3 N?N, f \ a a, = -I-2U"r+4(cx-ca> (□*-□„>) n\n, * \ g I (2) and {sfa*}~^;=€fa*+pfaA -2 N*N, 2 Re -(2Tr-^: \ M *o-T2< (3) where sg and Sj- are the gluon and fermion entropies, re­spectively. We measured screening masses for meson sources with quantum numbers of the tr, a, p, and ax. These measure­ments are a standard hadron spectrum calculation, ex­cept that the propagation is in the z direction. We used a wall source covering the entire z = 0 slice of the lattice, with the gauge fixed to a spatial Coulomb gauge which maximizes the traces of the x, y, and t direction links. After blocking five to ten measurements together to mini­mize the autocorrelations, we fit all the propagators to a single exponential using the full covariance matrix of the propagator elements. A quark mass can be defined from the divergence of the axial current [10,11]. The basic relation is a current algebra relation V < $r5^(0)$y5y^U)) =2mq < $y5ip{0)ipy5ft(x) > . (4) If we sum over x, y, t slices, and measure distance in the z direction, this becomes dz x,y,t = 2mq 2 iifty ^(0)rpy ^(x )) • (5) x,y,t We define II(z) as the pion correlator with a point sink, mz)=(W(0)^y^(z)) (6) and A (z) also the axial vector current correlator: A(z)={ fV(0)tpy5yzif>(z)) , (7) <w>= 4k Nf N?N, (1) and where Jf'XO) is the wall source at z =0. At long distances both II(z) and A (z) will fall off as exp(-m„z). There­fore we perform a simultaneous fit to the two propagators on a lattice periodic in the z direction using three param­eters, C, mT, and mq:3576 CLAUDE BERNARD et al. 49 n(z) - C sinh(mir)[exp( - m1rz) + exp[mv.(Nz - z)]} , where A(z) = C2m Jexp( - m^z)-exp[mw{Nz - z)] j (8) (9) The factor of sinh(m^) in Eq. (8) comes from using the lattice difference /(z + 1)-/(z -1) for the derivative in Eq. (5). Note that II(z) is periodic in z while A(z) is an- tiperiodic. We use the pointlike axial vector current i/»(z )YsYllip(z) rather than a point split current. These are quark masses in lattice units; to convert to continuum quark masses requires a lattice-to-continuum renormal­ization. See Ref. [17] for a discussion of this point. The quark propagator in the Landau gauge was also measured. This propagator has been studied with Kogut-Susskind quarks in Ref. [12]. We chose a source constant in the y direction and a 8 function in x, z, and t with only the real part of the first Dirac component nonzero (in the Weyl basis we use). Because of the 8 function all possible momenta in x, z, and t directions were excited. To distinguish among the different momen­ta we performed a Fourier transform of the propagator in x and t directions (taking into account that it has to have odd frequencies in t direction). This gives the propaga­tion of the quark in the z direction as function of kx and k„ i.e., the dispersion relation of the screening propaga­tor. In order to keep the amount of generated data at a reasonable level, the propagator was saved only for on- axis momentum values of kx and kt. This enabled us to measure the on-axis dispersion relation of the quark screening mass, in particular the screening mass difference of the quark and light doublers. The form to which the spatial propagator is fitted is usually motivated by the form of the free propagator. We suppose that at large distances, each separate momentum component of the spatial propagator resem­bles the corresponding free quark form, but with its own renormalized quark mass, or in this case of Wilson fer­mions, with a renormalized k. In momentum space the free Wilson propagator is G(k )= - -2/c]£ cos(pp f* (10) -+4/r 2sin (Pfj With our choice of the source the first Dirac com­ponent of the propagator, G t, is real: fl- -iK^COSlp^ fi 1- 2/c]jr cos(p^) j2 + 4/c22sin2(pM) 1m J m Then, for nonzero z values, Lz G(pup2,z,p0) = 2 exp[i2irkz /Lz]G(pup2,k,p0) k = 1 (11) =L, X ( 1 - 6k+4kA )2-4k(B + 1) 8/c( 1-6k+4kA )2 cosh[wa(z - Lz/2)] sinh( ma ) sinh( maL7 /2) (12) A= 2 sin [pM/2] , /i-1,2,0 B= 2 sin2[pJ . fi= 1,2,0 and ma =2arcsinh 4k2B + { 1 - 8k+4jC/4 )2 8#c( 1 - 6k+4kA ) 1/2 (13) (14) (15) At zero momentum (on a low temperature lattice where the lowest Matsubara frequency is close to zero) this rela­tion turns into ma =ln 1~6 K 2k (16) The mass vanishes when as expected. Inverting this for k gives 1 2exp[ma] + 6 ' which, for small masses, reduces to the naive relation, 1 K = 2 ma +8 (17) (18) that one expects looking at the terms of the Lagrangian. For large lattices the lowest doubler mass becomes ma doubler = lim In k-*1/8 1-2* 2k = ln[3] = l. 09861 . (19) For Kogut-Susskind fermions [12] the free propagator turns out to be a sum of two terms, having parts with an alternating sign in z direction. For Wilson fermions, with our choice of the source, the propagator is a single ex­ponential, or hyperbolic cosine, on a finite lattice. Furthermore, the sign of Gt at k -0 changes at kc ={. Therefore, measuring the sign of the propagator can be used as an indicator of whether k is effectively greater or less than kc . One can infer from Eqs. (13)-(15) that the only effect of finite spatial lattice size is the discretization of the mo­menta. For a given momentum all lattice sizes give the same value of the screening mass. For a smaller lattice, the range of allowed momenta is more restricted, of course. To be specific let us look at what happens with our lat­tice size: 82X20X4. This is shown in Fig. 1. At kc the lowest momentum screening mass is at its minimum. If one increases k the screening masses start to converge to a single value close to one at k=0.152. At this point the dispersion relation is flat. The sign of the propagator with this source depends on the momentum. Generally, the value of k at which the sign changes increases with the momentum. For zero momentum it occurs at kc; for the smallest nonzero momentum in our lattice size it takes place at /c=(2 + v/2 -v^J/l^O. 1374. The amplitude for the doubler does not change sign in this k range.49 NATURE OF THE THERMAL PHASE TRANSITION WITH . 3577 K FIG. 1. The spatial screening mass at different spatial mo­menta for free Wilson fermions as a function of k for an 82 X 20 X 4 lattice. At kc = j the higher masses are for higher momenta. For our lattice size, inserting the appropriate momenta to Eq. (15) one obtains the following quark screening masses at kc = j-. k =(0,0,0,ir/4), ma= 0.6610, A: = (ir/4,0,0,ir/4), ma=0.8906, k = (ir/2,0,0,ir/4), ma = 1.1171, (20) /c=(3ir/4,0,0,ir/4), ma = 1.2149 , k = (ir,0,0,ir/4), ma = l.24ll. For purely temporal momenta the free field screening mass is k = (0,0,0,±7t/4), ma= 0.6610 , (21) A: =(0,0,0,±3ir/4), ma = 1.0711. Hence, the temporal doubler's screening mass is smaller than that of the spatial doubler. III. SIMULATIONS AND RESULTS Simulations were run on the Intel iPSC/860 and Para­gon, and the nCUBE-2 at the San Diego Supercomputer Center, on the Thinking Machines Corporation CM5 at the National Center for Superconducting Applications, and on a cluster of RS6000 workstations at the Universi­ty of Utah. We used the hybrid Monte Carlo algorithm with two flavors of dynamical quarks in all our simula­tions [13]. The parameters of our runs are listed in Tables I, II, and III. For the 82 X 20 X 4 runs we used trajectories with a length of one unit of simulation time and made measure­ments after every second trajectory. The step size for these runs ranged, in the normalization of Ref. [14], from 0.033 for the largest 6/g2 and smallest k to 0.02 at the other extreme. Acceptance rates for these runs range from 70% to 90%, with an average over all the runs of 87%. For computation of the fermion force in the updat­ing and the propagators in the measurements we used the conjugate gradient algorithm with even-odd incomplete lower upper (ILU) preconditioning [15]. The conjugate gradient residual, defined as |M^Mx -b\/\b \ where M is the preconditioned matrix, b is the source vector, and x is the solution vector, was 10-6. Runs were made at 6/g2 = 5.3, 5.1, 5.0, and 4.9 with Nt= 4. At 6/g2 = 5.3 we also made a series of runs with N, = 6. At 6/g2=5.1 we ran two points with N, = S and at 6/g2 = 5.0 two points with Nt =6. We also ran a series of simulations at 6/g2 = 5.1 on 62 X 20 X 4 lattices to verify that the spatial size of the lattice was not seriously affecting our results. TABLE I. Table of runs at fixed 6/g2 with varying k. "(h)" and "(c)" indicate hot and cold starts. JV, Nx,y 6/g2 K Traj. Ignore dt Accept 4 8 4.9 0.180 650 100 0.02 0.88 4 8 4.9 0.181 320 100 0.02 0.94 4 8 4.9 0.182 810 100 0.02 0.90 4 8 4.9 0.1825 824(c) 100 0.02 0.86 4 8 4.9 0.1825 780(A) 100 0.02 0.88 4 8 4.9 0.183 624(A) 100 0.02 0.90 4 8 4.9 0.183 810(c) 600 0.02 0.87 4 8 4.9 0.184 610 100 0.02 0.93 4 8 5.0 0.173 540 100 0.025 0.86 4 8 5.0 0.175 500 100 0.02 0.92 4 8 5.0 0.177 400 100 0.02 0.90 4 8 5.0 0.178 474 100 0.02 0.94 4 8 5.0 0.180 630 100 0.02 0.92 4 8 5.0 0.182 256 100 0.02 0.86 6 8 5.0 0.175 240 80 0.02 0.94 6 8 5.0 0.180 360 60 0.0167 0.87 8 8 5.0 0.175 350 100 0.0167 0.913578 CLAUDE BERNARD et al. 49 At 6/g2 = 5.3 a series of short runs on 63X4 lattices was made for very large k. For reference we show a phase diagram for the relevant range of k and 6/g2 in Fig. 2. Previous work showed that as k increased from 0.16 to 0.19 along the N, =4 high temperature crossover line the pion mass de­creased, suggesting a closer approach to the high temper­ature transition in the chiral limit [3]. More recent work by Iwasaki et al., beginning from the 6/g2=0 limit, sug­gested that a high temperature transition for zero quark mass might be found at k~ 0.225 [16]. We have done a series of runs on 83X4 lattices in which we varied 6/g2 at /c = 0.20, 0.21, and 0.22 to extend the previous work. As expected, the number of conjugate gradient iterations re­quired in the updating increases as k increases in this range, and the size of the possible updating time step de­creases. Thus these runs have limited statistics. In Fig. 3 we show the plaquette and Polyakov loop as a function of 6/g2 for the various values of k. Notice that the transi­tion appears to be sharpest at k «0.19, becoming smooth­er for larger and smaller k. Even in those cases where the transition is very abrupt, we do not see the sorts of meta­TABLE II. Table of runs at fixed 6/g2 with varying k. "(h )" and "(c)" indicate hot and cold starts. Most of the TV,=4 runs were on 82 X 20 X 4 lattices. The run indicated with a t at k=0.168,6/g2=5.3 was done on a 82X40X4 lattice. N, 6/g2 K Traj. Ignore dt Accept 4 8 5.1 0.165 1120 100 0.025 0.89 4 8 5.1 0.167 2660 100 0.025 0.89 4 8 5.1 0.169 460 100 0.025 0.87 4 8 5.1 0.170 700 100 0.025 0.90 4 8 5.1 0.171 980 100 0.025 0.89 4 8 5.1 0.172 3380 100 0.025 0.87 4 8 5.1 0.173 500 100 0.025 0.85 4 8 5.1 0.175 460 100 0.025 0.86 4 8 5.1 0.177 1800 100 0.025 0.88 4 8 5.1 0.179 660 100 0.025 0.88 4 6 5.1 0.169 860 200 0.0333 0.82 4 6 5.1 0.170 1000 100 0.0333 0.85 4 6 5.1 0.717 1360 100 0.0333 0.81 4 6 5.1 0.172 1120 100 0.0333 0.83 4 6 5.1 0.175 720 100 0.0333 0.81 8 8 5.1 0.167 512 100 0.025 0.82 8 8 5.1 0.173 279 100 0.02 0.87 8 8 5.1 0.177 440 100 0.02 0.66 4 8 5.3 0.155 2400 100 0.0333 0.78 4 8 5.3 0.157 660 100 0.0333 0.76 4 8 5.3 0.158 1239 100 0.0333 0.89 4 8 5.3 0.159 660 100 0.0333 0.89 4 8 5.3 0.160 1777 100 0.025 0.88 4 8 5.3 0.161 480 100 0.025 0.88 4 8 5.3 0.162 480 100 0.025 0.83 4 8 5.3 0.163 680 100 0.025 0.87 4 8 5.3 0.164 460 100 0.025 0.90 4 8 5.3 0.165 720 100 0.025 0.88 4 8 5.3 0.166 660 100 0.025 0.87 4 8 5.3 0.167 912 100 0.025 0.86 4 8 5.3 0.168 540 100 0.025 0.89 4f 8 5.3 0.168 840 100 0.025 0.80 4 8 5.3 0.169 380 100 0.025 0.86 4 8 5.3 0.170 380 100 0.025 0.79 4 8 5.3 0.172 440 100 0.025 0.87 6 12 5.3 0.155 320 100 0.0177 0.88 6 12 5.3 0.160 552 60 0.0177 0.91 6 12 5.3 0.165 666 216 0.0177 0.85 6 12 5.3 0.166 1403 400 0.0177 0.84 6 12 5.3 0.167 760 302 0.0177 0.84 6 12 5.3 0.168 603 200 0.0177 0.8549 NATURE OF THE THERMAL PHASE TRANSITION WITH . 3579 stability and tunneling characteristic of strong first order transitions. We do find cases where equilibration takes a long time. The worst case was in the run at 6/g2=4.9 and k=0. 1825. In this case we have plotted two points, from hot and cold starts. These points are marked by ar­rows in Fig. 3. However, these two runs eventually con­verged to similar values, lying in between the values in the early parts of the runs. The time history of the Po­lyakov loop in these two runs is shown in Fig. 4. We now examine the 82X20X4 runs in more detail. Figure 5 shows the real part of the Polyakov loop as a function of k for the different values of 6/g2. For 6/g2 = 5.3 we also include values for N, =6 to show how the transition point moves as Nt increases. For all of these values of 6/g2 we see the expected sharp increase in the Polyakov loop at a value of k less than kc, where kc is the value at which the squared pion mass vanishes on a zero temperature lattice. We estimate kc at these values of 6/g2 from published values of kc in Refs. [7] and [4] and a recent measurement at 6/g2 = 5.3 by the HEMCGC group: Ke(5.3)=0.16794 [17]. From a quad­ratic fit to these values, shown by a line in Fig. 2, we find Kc(6/g2)=0.1687(2) at 5.3, 0.1795(4) at 5.1, 0.1861(12) at 5.0, and 0.1941(40) at 4.9. Although not a physical quantity, the number of conju­gate gradient iterations used in solving M Mx=b indi­cates how singular M is on the average. This quantity has been used as a probe of the physics in Ref. [16]. In Fig. 6 we show the average number of conjugate gradient iterations used in an updating step, where a linear extra- TABLE III. Table of runs at fixed k with varying 6/g2. "(h)" "(c)" indicate hot and cold starts. The acceptance rate gives the average over all runs in the sample kept for measurement, whether or not dt was changing during the runs. Nt 6/g2 K Traj. Ignore dt Accept 4 8 4.75 0.19 578 100 0.014286 0.912(13) 4 8 4.755 0.19 1504 750 0.014286 0.960(6) 4 8 4.76 (c) 0.19 837 500 0.014286 0.948(11) 4 8 4.76 (h) 0.19 362 100 0.014286 0.962(12) 4 8 4.32 0.20 156 50 0.02 0.83(4) 4 8 4.36 0.20 172 50 0.02 0.83(3) 4 8 4.40 0.20 188 50 0.02 0.79(3) 4 8 4.44 0.20 152 50 0.02 0.77(4) 4 8 4.48 0.20 368 50 0.02 0.69(3) 4 8 4.50 0.20 244 100 0.014286 0.84(3) 4 8 4.52 0.20 841 150 0.014286 0.773(16) 4 8 4.54 0.20 566 150 0.014286 0.72(2) 4 8 4.56 0.20 1324 200 0.014286 0.941(7) 4 8 4.60 0.20 365 50 0.014286 0.937(14) 4 8 4.64 0.20 244 50 0.02 0.959(14) 4 8 4.10 0.21 74 50 0.01 0.79(8) 4 8 4.20 0.21 267 50 0.005 0.90(2) 4 8 4.26 0.21 586 100 0.005 0.85(2) 4 8 4.28 0.21 478 100 0.005-»0.0025 0.82(2) 4 8 4.30 0.21 454 100 0.005 0.904(16) 4 8 4.32 0.21 227 50 0.005 0.94(2) 4 8 4.34 0.21 259 50 0.007 143 0.943(15) 4 8 4.36 0.21 281 100 0.0025->0.007 143 0.960(14) 4 8 4.40 0.21 197 50 0.0025->0.01 0.95(2) 4 8 4.44 0.21 249 50 0.007 143-*-0.01 0.977(9) 4 8 4.50 0.21 120 50 0.05->-0.01 0.94(3) 4 8 3.80 0.22 59 15 0.001 0.93(4) 4 8 3.90 0.22 56 30 0.002->0.0004 0.69(9) 4 8 3.96 0.22 39 25 0.002->0.0005 0.36(13) 4 8 4.00 0.22 119 80 0.004->0.002 0.50(8) 4 8 4.04 0.22 161 50 0.004 0.71(4) 4 8 4.06 0.22 119 40 0.003 333 0.78(5) 4 8 4.10 0.22 234 50 0.005 0.86(3) 4 8 4.20 0.22 90 50 0.007 143 0.98(3) 4 8 4.30 0.22 58 50 0.005->0.007 143 1.00(0) 4 8 4.40 0.22 90 50 0.01 0.95(3) 4 8 4.50 0.22 122 50 0.01->0.02 0.96(2)3580 CLAUDE BERNARD et al. 49 K FIG. 2. Phase diagram showing estimates for the high tem­perature transition and kc. Circles represent the high tempera­ture transition or crossover for N, =4, squares the high temper­ature transition for N,= 6, and diamonds the zero temperature kc. Previous work included in this figure is from Refs. [2,6,7,3,4]. We show error bars where they are known. For series of runs done at fixed k the error bars are vertical, while for series done at fixed 6/g2 the bars are horizontal. Points coming from this work are shown in heavier symbols. The solid lines are fits to k, for TV, =4 and to kc used in interpolating and extrapolating. 6/g‘ 6/g‘ FIG. 3. The plaquette and Polyakov loop as a function of 6/g2 for various values of k. The diamonds are previous results of Ref. [3] for <c = 0.12, 0.14, 0.16, 0.17, 0.18, and 0.19. For /c=0.12 and 0.14 data from long runs as well as some data from short runs collected while generating hysteresis loops are shown. The octagons at k-0.20, 0.21, and 0.22, are new results from 83X4 lattices. The squares come from runs on 82X20X4 lattices. These runs were done at fixed values of 6/g2 with vary­ing k. They have been mapped onto this figure by fitting the 6/g,2,K, line (with a fit shown as a line in Fig. 2), and moving the points in the K,6/g2 plane parallel to this line. Specifically, we plot the points at 6/g2ffectiv<. = 6/gJun-9(6/g,2)/3(f,(/frun-<c,). The fit for k, at 6/g2=5.3, 5.1, 5.0, and 4.9 is 0.1579, 0.1713, 0.1772, and 0.1827, respectively. time FIG. 4. Time history of the real part of the Polyakov loop for runs with hot and cold starts at 6/g2=4.9 and k=0.1825. polation of the last two time steps was used to produce a starting guess for the solution vector. For 6/g2 = 5.3 and N, =4 there is very little effect on the number of itera­tions at Kt. As 6/g2 is decreased for N, =4 there is an in­creasingly sharp peak in the number of iterations at kc. Notice also the sharp peak in the Nt = 6 results for 6/g2 = 5.3. Figure 7 shows the average plaquette in these runs. Our normalization is such that the plaquette is three for a lattice of unit matrices. The plaquette also shows a sharp rise as the high temperature crossover is passed. Notice that for 6/g2 = 5.3 we have results for TV, =4 and 6 show­ing that this increase is in fact due to the time size of the lattice, or the temperature. _ The chiral condensate iptfi is less useful for Wilson K FIG. 5. Expectation value of the Polyakov loop as a function of k for the various values of 6/g2. Results are shown for N, =4 for 6/g2 = 5.3, 5.1, 5.0, and 4.9 (octagons). For 6/g2 = 5.3 we also show results for N,= 6 (diamonds). The crosses along the 6/g2 = 5.1 line are results on a 62 X 24 X 4 lattice at 6/g2 = 5.1, to show that the spatial size of the lattice is not greatly affecting the results. The dotted symbols extending the 6/g2=5.3 line are short runs on a 63X4 lattice, showing that the behavior is smooth out to very large k. The vertical lines mark the zero temperature kc for 6/g2=5.2, 5.1, and 5.0, respectively. [For 6/g2 = 4.9, kc( 7' = 0)«0.194.]49 NATURE OF THE THERMAL PHASE TRANSITION WITH . 3581 K FIG. 6. Conjugate gradient iterations for updating step, as a function of k for the various values of 6/g2. Again the dia­monds are N, -6 results at 6/g2=5.3. K FIG. 8. Expectation value of i/n/z as a function of k for the various values of 6/g2. The dotted symbols for 6/g2=5.3 are short runs on a 63 X 4 lattice extending to k far beyond the zero temperature kc shown by the vertical line. quarks than for Kogut-Susskind quarks, since it does not go to zero in the high temperature phase without difficult subtractions. Nevertheless we plot it in Fig. 8. There is a clear drop in ifnf> as the high temperature transition is crossed. This drop increases dramatically as 6/g2 de­creases. Perhaps the most physically relevant observable is the entropy. In Fig. 9 we plot TXs in units of a ~4. To give an idea of the normalization of this graph, for eight gluons in free field theory on an 8X8X20X4 lattice the gluon entropy would be 7sglue free=0.040a ~4, while for two flavors of free Wilson quarks at k=kc =0.125 the en­K FIG. 7. Expectation value of the plaquette as a function of k for the various values of 6/g2. Here we included values for larger N, to emphasize the effect of the temperature. The dot­ted symbols for 6/g2 = 5.3 are short runs on a 63X4 lattice ex­tending to k far beyond the zero temperature kc shown by the vertical line. tropy would be 7iquarkjree =0.125a ~4. The effects of the lattice spacing and spatial size are very large here; in the continuum with infinite spatial extent these numbers are 0.027 and 0.036, respectively. Strangely, when we divide the entropy into gauge and fermion parts as in Eqs. (2) and (3) we find that the gauge entropy is comparable to the fermion entropy instead of much smaller as would be the case with free fields on a lattice of this size. The quark mass defined by the divergence of the axial pion propagator is plotted in Fig. 10. When this quark mass was small we had great difficulty in getting good fits to the forms in Eqs. (8) and (9). This is expected, because when the quark mass is small the amplitude for the prop- K FIG. 9. Entropy (actually TXsa4) as a function of k for the various values of 6/g2. Again we show the 62X24X4 results for comparison.3582 CLAUDE BERNARD et al. 49 K FIG. 10. Quark mass from the axial current as a function of k for the various values of 6/g2. Points marked with question marks indicate runs where we were unable to get consistent fits as a function of distance. The plus signs on the m, =0 line are the zero temperature kc for 6/g2 = 5.3, 5.1, and 5.0. agator A(z) is very small. Additionally, there is a ten­dency for the effective quark mass, or the quark mass coming from a fit over a short distance range, to increase with distance from the source. In cases where we were unable to get a fit with a satisfactory x2 or where the quark mass was not convincingly independent of dis­tance, we plot the point with a question mark in Fig. 10. To pursue this further we ran one of the difficult points, 6/g2 = 5.3 and *c=0.168, on a 82 X 40 X 4 lattice, allowing us to measure the ratio out to a distance of 20. Figure 11 summarizes the results. In this figure we show the effective pion mass obtained from IKx) and A(z) by fitting to two successive distances, and the quark mass FIG. 11. Pion effective screening masses from n(z) (circles) and from A(z) (squares), and the effective quark screening mass from their ratio. The results are from an 82X40X4 lattice with 6/g2 = 5.3 and k=0.168. obtained from simultaneously fitting both propagators at the two successive distances (a one degree of freedom fit). Unfortunately, in all other cases the lattice was only 20 sites long and we have to draw conclusions from dis­tances less than ten. In Fig. 10 we see that when the TV, =4 lattice enters the high temperature regime the pointlike axial current quark mass no longer agrees with the low temperature lattices (N,= 6 and 8). The plusses at rnq- 0 in Fig. 10 are estimates for the zero tempera­ture kc. The axial current quark masses go through zero at k less than the zero temperature kc. When the axial current quark mass vanishes, the system is in the high temperature phase for j8> 5.0, while at /3=4.9 k, appears to coincide with the point where the axial current quark mass vanishes, within experimental uncertainty. Note, however, that the pion screening mass in the confinement phase is still nonzero at the transition point at /?=4.9. In Fig. 12 we show the squared pion screening masses in these runs. Again we see an increasingly sharp dip at Kt as 6/g2 decreases and k increases. The appearance of the cusp at /3=5.1 coincides with the beginning of the re­gion where the transition is abrupt. Screening masses for the n, p, a, and ay mesons are shown in Figs. 13, 14, and 15. In all cases we see the screening masses coming to­gether as the high temperature transition is crossed. However, we do not see any indication that the it-a or p~a{ splittings in the high temperature regime are de­creasing as 6/g2 decreases. Although the smaller pion masses in the cold regime suggest that chiral symmetry is being approached as we move toward smaller 6/g2 along the k, line, we do not see this trend in the high tempera­ture screening masses. Also notice that there are nonzero splittings between the parity partners at the points where the axial current quark mass is zero. Thus the vanishing of this quark mass is not an indicator for complete chiral K FIG. 12. Pion screening mass squared as a function of k for the various values of 6/g2. Again, the circles are for N, =4, the diamonds for N, =6, and the crosses for N, = 8. The bursts are zero temperature pion masses from the HEMCGC Collabora­tion at 6/g2 = 5.3.49 NATURE OF THE THERMAL PHASE TRANSITION WITH . 3583 K FIG. 13. Meson screening masses for 6/g2=5.30. The points connected by solid lines are for N, =4 and the points con­nected by dashed lines for N, = 6. symmetry restoration in the system. To investigate the contributions of the doublers to thermodynamic quantities such as the entropy we mea­sured the effective masses from the quark propagator in Landau gauge at a few values of k and 6/g2. We find that fitting the quark screening propagators is more difficult than fitting the meson propagators. In part this is because the quark propagators fluctuate more from configuration to configuration. There also seems to be a systematic trend toward larger effective quark masses at larger distances. With these caveats, the masses of the physical quark and the lightest doublers are given in Table IV. The fits were selected by choosing the largest fit range that gives an acceptable confidence level. The ranges and confidence levels are also given in Table IV. K FIG. 14. Meson screening masses for 6 /g 2= 5.10 at N,= 4. tc FIG. 15. Meson screening masses for 6/g2=4.90 at N, =4. The two points at k=0. 1825 are from cold and hot starts. IV. CONCLUSIONS The most naive expectation regarding the thermo­dynamics of two flavors of Wilson quarks at fixed N, is that there would be a line in the k,P plane at which a confinement-deconfinement transition occurs, that the transition would be smooth (crossover or second order), that the pion mass would smoothly decrease along that line, and that at some point, possibly corresponding to the point where the transition line crossed the zero tem­perature Kc-p line, the pion mass would go to zero. At that point one would have a finite temperature confinement-deconfinement or chirally restoring transi­tion analogous to that seen in staggered fermions. Simple arguments [9] would put this point around 0=5.0 at Nt= 4. These naive expectations are not borne out by the data. The chiral limit is reached at a very small (3 value if it is reached at all. However, near (3=5.0, N,=4 Wilson thermodynamics displays a number of features which have no analogs in staggered fermion systems. The tran­sition becomes very sharp, though not first order as far as we can tell. A cusp in the pion screening mass appears as one crosses from the confined to the deconfined phase. The axial vector quark mass becomes strongly N, depen­dent at this point and for small N, does not go to zero at its zero temperature k value (at fixed [3). The sharp tran­sition persists down to j3=4.5,/c=0.20 or so, at which point it once again becomes smooth. As far as we can tell, the zero temperature k=kc point plays no role in any N,= 4 effects we have observed. It is tempting to speculate that the crossover line in the K,6/g2 plane is close to some phase boundary where the transition is steepest. We are currently exploring this re­gion with N,= 6, where preliminary results indicate a change in the nature of the high temperature transition around this value of k. Indicators for the nature of the high temperature phase3584 CLAUDE BERNARD et al. 49 TABLE IV. The screening masses for the quark and the lightest doublers. AmasU) is the difference of the spatial (temporal) dou­bler screening mass to the quark screening mass at the lowest momenta. The sign is for G(k, =ir/4). The fits were done simultane­ously to all three propagators taking into account the cross correlations. The confidence level q and the range of each fit is also displayed.________________________________________________________________________________________________ N, K P Sign ma(0,ir/4) ma(TT,v/A) ma{ 0,3ff/4) Amas A ma, 9 Range 4 0.165 5.10 + 1.13(4) 2.3(4) 1.8(2) 1.1(4) 0.7(2) 0.47 3-10 4 0.167 5.10 + 1.13(4) 2.4(5) 1.8(3) 1.3(5) 0.6(3) 0.14 3-10 4 0.172 5.10 + 0.97(16) 1.5(3) 1.5(2) 0.5(3) 0.6(4) 0.43 4-10 4 0.177 5.10 - 1.05(9) 1.28(15) 1.33(11) 0.23(19) 0.27(15) 0.69 4-8 4 0.155 5.30 + 1.06(2) 1.63(13) 1.65(11) 0.57(13) 0.60(11) 0.68 3-9 4 0.160 5.30 + 0.89(4) 1.38(9) 1.22(5) 0.49(9) 0.34(7) 0.84 3-10 4 0.167 5.30 - 0.92(6) 1.51(9) 1.36(7) 0.59(11) 0.44(10) 0.57 3-9 give a somewhat mixed picture. It is clear from the meson screening masses and from that chiral symme­try is at least partially restored at high temperature. While the axial current quark mass goes to zero the tt - a and p - splittings in the screening masses remain nonzero. Quark propagators in the Landau gauge sug­gest a large constituent quark mass at the transition, at least for 6/g2 = 5.3 and 5.1. This is consistent with ear­lier work [3] where at the JVt=4 crossover point near these (j3,/c) values the zero temperature pion was found to be quite heavy. Notice that the series of runs at 6/g2 = 5.3 extends to k significantly larger than the zero temperature kc, and there is no noticeable effect on any of the measured quan­tities when this kc is crossed. (In fact, we have done short runs on 83X4 lattices for k as large as 0.19 at 6/g2 - 5.3 and seen no effects.) The sign of the propaga­tor of the zero momentum quark, shown in Table IV, is consistent with the sign of the axial current quark mass. Both of these quantities are behaving in the way one would expect in a free field theory at k>kc. ACKNOWLEDGMENTS These calculations were carried out on the iPSC/860, the Paragon and the nCUBE-2 at the San Diego Super­computer Center, on the CM5 at the National Center for Supercomputing Applications, on 15 IBM/RS6000 workstations in the Physics Department at the University of Utah, on an IBM/RS6000 cluster at the Utah Super­computing Institute and on our local workstations. We are grateful to the staffs of these centers for their help. We also thank Tony Anderson and Reshma Lai of Intel Scientific Computers for their help with the Paragon. We would like to thank Akira Ukawa and Frithjof Karsch for helpful discussions. Several of the authors have en­joyed the hospitality of the Institute for Nuclear Theory, the Institute for Theoretical Physics, and the UCSB Physics department, where parts of this work were done. This research was supported in part by Department of Energy Grants Nos. DE-2FG02-91ER-40628, DE- AC02-84ER-40125, DE-AC02-86ER-40253, DE-FG02- 85ER-40213, DE-FG03-90ER-40546, DE-FG02-91ER- 40661, and National Science Foundation Grants Nos. NSF-PHY90-08482, NSF-PHY93-09458, NSF-PHY91- 16964, NSF-PHY89-04035, and NSF-PHY91-01853. APPENDIX Expressions for the energy, pressure, and iptfi are found by differentiating the partition function with respect to l/T, volume, and quark mass, respectively. First, we write the action with adjustable lattice spacings in all directions. Introducing dimensionless parameters aM, we write the lattice spacing in the (x direction as afi=aafl. Clearly this is redundant, since we have five parameters a and ap to specify four lattice spacings, but it is con­venient and symmetric. In the conventional notation of Karsch, 1=0,/a,, where all the spatial a's are the same. When we are done taking derivatives, all the a will be set to one. The partition function is Z- f [dU]e s„+sf (Al) where the gauge action is 2 axava,at 2 2 4- 2 2 Dev <A2) x n> v 5 {tv where DMV is the plaquette in the fiv plane normalized to three for unit matrices. We allow a different gauge cou­pling gfLV in each plane. The fermion action is nf Trln MfM where M = l-#c2--0, (A3) (A4) where =< i+tv > ty * nytX+P+< i - )£/; u - a )8yiX (A5) The a,, in the coefficient of takes care of the dimen- fj. fl sional scaling of the first derivative. Notice that we have made a somewhat arbitrary choice in M when we scaled the irrelevant second derivative part with afl in the same way that we scaled the first derivative part. The in the coefficient must be adjusted to get correlation functions to be Euclidean invariant. Its role is similar to the49 NATURE OF THE THERMAL PHASE TRANSITION WITH . 3585 Karsch coefficients Ca and CT in the gauge action. Presumably eM has a power series expansion in g just as Ca and CT. Once again we have more parameters than we need: four e and k for four directions. This parame­trization is convenient because it includes the customary k and, later, kc as parameters. We can fix the ambiguity well enough for our purposes by requiring that e^= 1 when all the a are equal. In other words, if all direc­tions are scaled by the same factor the only thing that changes is k. Let kc be the value of k at which the pion mass and quark mass vanish, at least on an infinite lattice. Follow­ing free field theory, we introduce a quark mass 2 ma =k 1-kc 1 so that M - 1 ' + 2ma 2ma +k. *- 2 - (A6) (A7) Here k~ 1 will be a function of the couplings g^v and the scale factors aIn free field theory, = (e^-l in free field theory.) We find the energy, pressure, and ipip by differentiating the partition function: e= - 1 31nZ V 3/3 V,m const 1 31nZ 18V dm (A8) (A9) (A10) 0, V const Here V is the volume, V=aiY\iNiai=a2N^a3s. /3 is the inverse temperature, P=aNtat. Here the energy and pressure derivatives are taken with m constant, rather than with k constant. This is because kc depends on the afl, so that if we distort the lattice while holding k fixed, the quark mass, and every physical mass, will vary sharp­ly- 1 0/3 N,a da, and 1 dV 31V?a3a? 3as (All) (A 12) Alternatively, it may be easier to vary the volume by varying only one of the spatial lattice spacings 3 1 3 dV N^a3axay 3az (A 13) The gauge energy and pressure are standard [18-20]. Following Ref. [18], we define two derivatives: C = -2 MV 3a, (A 14) where k is not one of or v and C = a 9a, (A 15) where k is one of /z or v. Because stretching both the time and space directions is equivalent to changing the lattice spacing, Ca and CT are related to the 13 function: gab2=' ~2C„ -2Ct -2 ain(a) (A 16) The contributions to the energy and pressure from Sg and Sj- add. Doing the differentiation, and then setting the aM to one, the gluon energy is N, ?() T(nJ,-nB)+6cffnB+6c. rO*) • (A 17) Here and are the space-space and space-time pla- quettes, again normalized to three for a lattice of unit matrices. For the gluon pressure, we find ) -2Ca(DSJ+2n„) -2CT(2D„+D (A 18) We also consider the linear combination e+p, the en­tropy: ,4_L aN, = ega*+pga* ■+4(C T-C„) (A 19) The entropy is obviously zero at T=0. Just as the gauge couplings vary with the lattice spac­ings, k~ 1 and vary with the lattice spacings as we try to hold m fixed. There is an explicit dependence of k~ 1 on m plus a dependence of k~ 1 on g, where g is varying with the a^. Now it clearly does not matter which direc­tion we stretch the lattice, since k~x is defined on the infinite lattice, so 3 k. -l a k -i 3a, Therefore 3*71 3 ar 3 ac = 3- 3 k. -l V-i 3a, (A20) (A21) There are two independent derivatives of the eM, analo­gous to Ca and Cr. Define3586 CLAUDE BERNARD et a l 49 B = and de B"=3Z' ^ To compute ipip we can just set and eM to one at the beginning: (A22) a rptp (A23) iV, N?N, 3 Tr lnMfM dma Nr N?N, Tr M1M 3„ , .,+ sm -------M+M ------- dma dma Nf N?N, Tr 1 3 M t + 1 3 M Mf 3ma M dma (A24) The two parts are complex conjugates, so keep only one and take twice the real part. Using Eq. (A7), a 3rpip= -j-------- Re/Tr- N?Nt 2 \ M (Kr 1 + 2ma )2 2 ma +kc 1 ^ + 2 1 Kr ' + 2ma 1 AkN f / JVjty 2 Re^Tr -1 )■ (A25) where in the last step we used \/kc 1 +2ma =k. Looking at the derivation shows that the 1 in Eq. (A25) comes from differentiating the 1 /K~2 + 2ma outside the parentheses in Eq. (A7). Had we taken the fermion matrix to be M = this term would be absent. Since this latter form is closer to the usual continuum Lagrangian, we prefer =(Trir) _ i 4 kN, \pip=~-------zr1- Re( N\Nt 2 as our expression for ipip. Now for the fermion energy, differentiate Eq. (A4) and then set and efL to one: 1 3 Mu (A26) (A27) ej-a -1 *//. 7Vs3iV; 2 V -2 Nf r JV/JV, 2 1 2 kNf 2 Re Tr Re 3ar ■1 (kc x+2ma )2 - + 1 3a, kc ' + 2wa + be + dKc 1 Re(- N?N, 2 \ 3a, 1-- M Mt (A28) (A29) where in the last step we used 1 /kc 1 + 2ma =k. For the fermion pressure, pfa 1 ^/2ReTr^-^ 3 N*Nt 2 3a, r-Re/- 3JV/JV, 2 \ Af 1 (/c_ ' + 2ma )2 3/c,. 1 3 ar + - 1 '+2ma de + 3kc 1 > (A30) 3Ns3Nt 2 Re( 3 ac 1- MT (A31)49 NATURE OF THE THERMAL PHASE TRANSITION WITH . 3587 Just as for ipip, the 1 term in ey and Pf comes from differentiating the overall factor of k_1. It can be includ­ed or not, as desired. It will cancel when the finite parts of the energy and pressure are calculated by subtracting the zero temperature result from the nonzero tempera­ture result. However, the (1 /da,) term will not cancel out, since xl>ip is temperature dependent. In practice (Trl/M^) is fairly close to (Trl)=4X3, so numerically it may be best to leave the 1 in. Then we would use 1- M K0' t M (A32) to express the energy and pressure just in terms of the ex­pectation values of the spatial and temporal components of#f. Much of the difficulty cancels out if we look at the en­tropy, or sum of energy and pressure: 4 1 __ 4 i 4 SfU ^~€fa +Pf° aN, -2 "Nf / 1 -r-----Re(Tr--?[ 1 + (BT-Ba )] N*Nt 2 \ Mf T ° X >' The relation 3/c da, 1 d*-' - = 3- 3a( (A3 3) (A 34) resulted in all the derivative terms canceling. (Remember, by 3/3as we mean 3 /dax + d/day + 3 /da2 vary all the spatial a,- together.) As usual, the entropy is obviously zero at T = 0, where 0- Since no zero temperature subtraction is required for the entropy, the terms involving BT and Ba will be higher order in g2 than the "1" term, and we have neglected them in Eq. (3). Obviously, the big problem in getting the energy and pressure separately is to find 3*7'/3aM and BT and Ba. From Euclidean invariance, dK~l/dafi is independent of IiThe variable k~ 1 depends on the ap in two ways. First, there is an "explicit" dependence. From examining the fermion matrix, Eq. (A4), we see that if all the a^ are scaled together with g held fixed, kc is proportional to a. Thus 3k-, -l 3 a„ (A35) explicit Secondly, there is an "implicit" dependence of k~ on coming from the fact that kc depends on g-2, and we ad­just the g2v as we adjust the ap. Again, Euclidean invari­ance says 3*71 da,, 1 3/c„ implicit 4 3 ln( a) (A36) so only the (3 function appears: 3/c, 3a„ 1 3 K, -l implicit 4 3 ln(a) 1 a<cc 1 36/g2 4 36/g2 31n(a ) ' (A37) We could proceed by estimating 3x'71/36/g2 from our data at various values of 6/g2, or from correlations of the hadron propagators with the plaquette. Similarly, we could take a /3 function either from perturbation theory or from some set of lattice simulations. We will not solve this problem here, so we will only quote the entropy rath­er than the energy and pressure separately. [1] M. Fukugita, S. Ohta, and A. Ukawa, Phys. Rev. Lett. 57, 1974(1986). [2] A. Ukawa, in Lattice '88, Proceedings of the International Symposium, Batavia, Illinois, 1988, edited by A. S. Kron­feld and P. B. Mackenzie [Nucl. Phys. B (Proc. Suppl.) 9, 463 (1989)]. [3] HEMCGC Collaboration, K. M. Bitar et al., Phys. Rev. D 43, 2396(1991). [4] C. Bernard, T. A. DeGrand, C. DeTar, S. Gottlieb, A. Krasnitz, M. C. Ogilvie, R. L. Sugar, and D. Toussaint, Phys. Rev. D 46, 4741 (1992). [5] Y. Iwasaki, K. Kanaya, S. Sakai, and T. Yoshie, in Lattice '92, Proceedings of the International Symposium, Amster­dam, The Netherlands, 1992, edited by J. Smit and P. van Baal [Nucl. Phys. B (Proc. Suppl.) 30, 327 (1993)]; in Lat­tice '93, Proceedings of the International Symposium, Dal­las, Texas, 1993, edited by T. Draper et al. [Nucl. Phys. B (Proc. Suppl.) 34, 314 (1994)]. [6] K. Bitar, A. D. Kennedy, and P. Rossi, Phys. Lett. B 234, 333 (1990). [7] R. Gupta, A. Patel, C. Baillie, G. Guralnik, G. Kilcup, and S. Sharpe, Phys. Rev. D 40, 2072 (1989). [8] K. Bitar et al., in Lattice '92 [5], p. 401. [9] A. Hasenfratz and T. A. DeGrand, Phys. Rev. D 49, 466 (1994). [10] M. Bochicchio, L. Maiani, G. Martinelli, G. Rossi, and M. Testa, Nucl. Phys. B262, 331 (1985); L. Maiani and G. Martinelli, Phys. Lett. B 178, 265 (1986). [11] S. Itoh, Y. Iwasaki, Y. Oyanagi, and T. Yoshie, Nucl. Phys. B274, 33 (1986); Y. Iwasaki, L. Kanaya, S. Sakai, and T. Yoshie, Phys. Rev. Lett. 67, 1494 (1991). [12] G. Boyd, Sourendu Gupta, and F. Karsch, Nucl. Phys. B385, 481 (1992). [13] S. Duane, A. D. Kennedy, B. J. Pendleton, and D. Roweth, Phys. Lett. B 195, 216 (1987). [14] S. Gottlieb, W. Liu, D. Toussaint, and R. L. Sugar, Phys. Rev. D 35, 2531 (1987). [15] T. DeGrand, Comput. Phys. Commun. 52, 161 (1988); R. Gupta, A. Patel, C. Baillie, G. Guralnik, G. Kilcup, and S. Sharpe, Phys. Rev. D 40, 2072 (1989); T. DeGrand and P. Rossi, Comput. Phys. Commun. 60, 211 (1990). [16] Y. Iwasaki, K. Kanaya, S. Sakai, and T. Yoshie, in Lattice3588 CLAUDE BERNARD et al. 49 '91, Proceedings of the International Symposium, Tsuku- ba, Japan, 1991, edited by M. Fukugita et al. [Nucl. Phys. B (Proc. Suppl.) 26, 311 (1992)]; Phys. Rev. Lett. 69, 21 (1992). [17] K. M. Bitar et al., this issue, Phys. Rev. D 49, 3546 (1994). [18] F. Karsch, Nucl. Phys. B205, 285 (1982). [19] J. Engels, F. Karsch, I. Montvay, and H. Satz, Phys. Lett. 101B, 89 (1981); Nucl. Phys. B205, 545 (1982); T. Celik, J. Engels, and H. Satz, Phys. Lett. 129B, 323 (1983). [20] R. Trinchero, Nucl. Phys. B227, 61 (1983).