A previously introduced multi-boson technique for the simulation of QCD with
dynamical quarks is described and some results of first test runs on a
$6^3\times12$ lattice with Wilson quarks and gauge group SU(2) are reported.Comment: 7 pages, postscript file (166 KB

B. Bunk

B. Jegerlehner

Duane

Fabricius

Fried

Geradin

H. Simma

Jegerlehner

K. Jansen

Kennedy

Lüscher

M. Lüscher

Meyer

Press

R. Sommer

English

arXiv

Crossref

Nuclear Physics B - Proceedings Supplements

A new simulation algorithm for lattice QCD with dynamical quarks

A previously introduced multi-boson technique for the simulation of QCD with dynamical quarks is described and some results of first test runs on a 6^3\times12 lattice with Wilson quarks and gauge group SU(2) are reported.A previously introduced multi-boson technique for the simulation of QCD with dynamical quarks is described and some results of first test runs on a $6~3\times12$ lattice with Wilson quarks and gauge group SU(2) are reported

Bunk, B.

Jansen, K.

Jegerlehner, B.

Luscher, M.

Simma, H.

Sommer, R.

CERN Document Server

1A new simulation algorithm for lattice QCD with dynamical quarksB. Bunka, K. Jansenb, B. Jegerlehnerc, M. Luscherb, H. Simmaband R. SommerdaHumboldt Universitat, Institut fur PhysikInvalidenstrasse 110, D-10099 Berlin, GermanybDeutsches Elektronen-Synchrotron DESYNotkestrasse 85, D-22603 Hamburg, GermanycMax-Planck-Institut fur PhysikFohringerring 6, D-80805 Munchen, GermanydCERN, Theory DivisionCH-1211 Geneve 23, SwitzerlandA previously introduced multi-boson technique for the simulation of QCD with dynamical quarks is describedand some results of rst test runs on a 63 12 lattice with Wilson quarks and gauge group SU(2) are reported.1. INTRODUCTIONThe basic idea of the algorithm proposed inref.[1] is to map lattice QCD to a local bosonictheory and to simulate the latter using standardtechniques. The method is generally applicable,but attention is here restricted to the case of twodegenerate avours of Wilson quarks on a four-dimensional lattice with periodic boundary con-ditions in all directions. We shall rst recall thebasic steps leading to the bosonic theory and thenpresent some results from test runs on small lat-tices with gauge group SU(2).2. EQUIVALENT BOSONIC THEORYIn terms of the gauge covariant forward andbackward lattice derivatives, rand r, theWilson-Dirac operator D may be written asD =123X=0f(r+r)  rrg ; (1)where the lattice spacing has been set equal to 1for convenience. Our starting point is the eec-tive distributionPe[U ] / [det(D +m)]2e Sg[U ](2)Talk given by M.L. at the International Symposium onLattice Field Theory, Sept. 27{Oct. 1, 1994, Bielefeldof the gauge eld U which one obtains after inte-grating over the quark elds. Sg[U ] denotes theusual plaquette action and the bare quark massm is related to the Wilson hopping parameter Kthrough K = (8 + 2m) 1.It is easy to show that the lattice Dirac opera-tor satisesDy= 5D5; kD +mk  8 +m: (3)The operatorQ = 5(D +m)=[cM(8 +m)]; cM 1; (4)is hence hermitean and all its eigenvalues are be-tween  1 and 1. In the following it will be advan-tageous to work with Q rather than D + m andso we rewrite eq.(2) in the formPe[U ] / detQ2e Sg[U ]: (5)The constant cM 1 will be xed later when wediscuss the simulation of the bosonic theory.We now choose a sequence of polynomials P (s)of even degree n such thatlimn!1P (s) = 1=s for all 0 < s  1 (6)(an explicit example is given below). From theabove we then infer thatdetQ2= limn!1detP (Q2) 1: (7)2We may, furthermore, choose the polynomials sothat their roots z1; z2; : : : ; zncome in complexconjugate pairs with non-zero imaginary parts.As a consequence the root factorization may bewritten in the formP (Q2) = constant nYk=1(Q  k)2+ 2k; (8)wherek+ ik=pzk; k> 0: (9)After inserting this product in eq.(7) the deter-minant factorizes and each factordet(Q  k)2+ 2k 1(10)can be represented as a Gaussian integral oversome auxiliary bosonic eld k.In this way one establishes the identityPe[U ] = limn!11ZbZD[]D[y] e Sb[U;]; (11)whereSb[U; ] = Sg[U ] +PxPnk=1j(Q  k)k(x)j2+ 2kjk(x)j2(12)(Zbis a normalization constant independent ofthe gauge eld).We have thus shown that lattice QCD is a limitof a purely bosonic local theory. Note that theelds kare coupled to the gauge eld but notamong themselves. Their action is strictly pos-itive and the bosonic theory is, therefore, suit-able for numerical simulation. Since only a rela-tively small number of scalar elds can be storedin memory, it is however important to make surethat the limit n!1 is reached rapidly. In otherwords, the polynomial P (s) must be chosen withcare so as to obtain a good approximation of 1=salready for low values of n.3. CHOICE OF P (s)The polynomial dened below approximatesthe function 1=s with a uniform relative error inthe range "  s  1, where " is an adjustableparameter. It is somewhat simpler than the poly-nomial proposed in ref.[1], but so far seems to doequally well in practice.Let us introduce a scaled variable u and an an-gle  throughu = (s   ")=(1  "); cos  = 2u  1: (13)The Chebyshev polynomial Tr(u) of degree r isgiven byTr(u) = cos(r): (14)For specied n and " we now deneP (s) =1 +  Tn+1(u)=s; (15)where the constant  is chosen such that thesquare bracket vanishes at s = 0 (and so is di-visible by s).0 .2 .4 .6 .8 1-.01-.0050.005.01R(s)sFigure 1. Relative deviation R(s) of the polyno-mial (15) from 1=s for " = 0:05 and n = 16.In the range "  s  1 the relative t errorR(s) = [P (s)  1=s] s (16)of this polynomial satisesjR(s)j  21 p"1 +p"n+1; (17)i.e. the convergence is exponential with a rateroughly equal to 2p". When s < " the polyno-mial continues to converge with an s{dependentexponential rate which approaches zero in thelimit s! 0. Moreover, as can be seen from g. 1,R(s) is oscillating around 0 so that the errors may3be expected to partly cancel when the productover all eigenvalues of Q2is taken.The roots of P (s) are given byzk=12(1 + ") 12(1 + ") cos2kn+1  ip" sin2kn+1(18)(k = 1; 2; : : :; n). They lie on an ellipse aroundthe spectral interval 0  s  1 and satisfyImzk6= 0 (recall that n is assumed to be even).The transformation to the bosonic theory thusworks out in the way explained above.4. NUMERICAL SIMULATIONThe bosonic theory with action Sb[U; ] can besimulated straightforwardly using local heatbathand over-relaxation algorithms. A few technicalissues are addressed here and further details arediscussed in ref.[2].We rst consider the updating of the \matter"eld k. If all eld variables are kept xed exceptkat point x, the action assumes the quadraticformSb[U; ] = [k(x)]yAkk(x)+[bk]yk(x) + [k(x)]ybk+ constant;(19)where Akis an easily calculable constant positivematrix. The calculation of the (constant) vec-tor bkrequires the computation of [Qk](x) and[Q2k](x). To obtain the latter one also needsto calculate [Qk](y) at all neighbouring points yof x. The subroutine which applies Q to a giveneld at a given point must hence be called 10times. This part of the simulation thus tends tobe rather expensive.There is, however, a more ecient way to pro-ceed. One rst computes [Qk](z) at all latticepoints z and stores the result in an auxiliary ar-ray. After that one passes through the lattice insome order and updates the eld k. At any givenpoint x the vector bkis obtained easily by apply-ing Q to the auxiliary eld. The latter must becorrected at z = x and all neighbouring points assoon as k(x) has been replaced by its new value.This operation costs about as much as applyingQ once at a given site. In this way the total workis reduced to applying Q three times per latticepoint.Concerning the updating of the gauge eld, wenote that the action reduces toSb[U; ] = Re tr fU (x; )Fg+ constant; (20)if all eld variables are held xed with the excep-tion of the link matrix U (x; ). The force F isa sum of two contributions, the usual one fromthe pure gauge eld action and the other fromthe coupling to the \matter" elds. The widelyknown algorithmsmay thus be used to refresh thegauge eld. The computation of the force how-ever requires an eort proportional to n, i.e. agauge eld update sweep tends to be as expen-sive as the updating of all elds k.We nally remark that the autocorrelationtimes for heatbath updating of the elds kwithk near n=2 tend to be rather large if the parametercMin eq.(4) is set equal to 1 (as in ref.[1]). Theprobable cause for this is that the correspondingroots zkare rather close to the top s = 1 of thespectrum of Q2. The latter can be lowered bychoosing (say) cM= 1:1 and the problem thendisappears.5. LOW-LYING EIGENVALUES OF Q2The parameter " should be chosen such thatmost eigenvalues  of Q2are in the range "   1 where P () is a uniformly good approx-imation to 1= (if n is large enough). To checkwhether this condition is actually met in any par-ticular simulation, one must compute the lowesteigenvalue 0of Q2for a representative sample ofgauge eld congurations.There are various methods to compute the low-lying eigenvalues of large sparse matrices. A par-ticularly safe way to obtain 0is to minimize theRitz functional( ) = kQ k2=k k2(21)in the space of all (non-zero) quark elds  , us-ing a conjugate gradient algorithm [3{6]. Highereigenvalues may then be computed recursively byminimizing the Ritz functional in the space or-thogonal to the approximate eigenvectors alreadyfound.4It is our experience that this method works wellin the case of the lattice Dirac operator. In partic-ular, the degeneracy of the eigenvalues is correctlyobtained and the attainable relative accuracy ofthe computed eigenvalues is completely satisfac-tory on machines with only single precision arith-metic, provided one employs an accurate summa-tion method to compute scalar products of elds[7]. A rigorous bound on the nal numerical accu-racy may actually be deduced from the stoppingcriterion of the algorithm.6. TEST RESULTSWe now present some results from simulationsof the bosonic theory on a 63 12 lattice withSU(2) gauge group and bare parameters  = 2:12and K = 0:15. At this point the mass mofthe  meson is roughly equal to 1:28 (in latticeunits), while for the pion mass mone obtainsm= 0:93m. The spatial size L of the latticesatises mL ' 7 so that nite volume eects areexpected to be small.0 .0005 .001 .0015 .002 .002505001000150020002500dp(0)=d00Figure 2. Probability distribution of the lowesteigenvalue 0of Q2. The full line is from sim-ulations of the bosonic theory with n = 60 and" = 0:0015, while the dashed line represents theresult from the Hybrid Monte Carlo run.In our simulations of the bosonic theory a fullupdate cycle involved 1 heatbath plus 1 over-relaxation sweep for all matter elds kfollowedby 1 heatbath sweep for the gauge eld [8,9]. Foreach choice of n and " we have performed at least128'000 and sometimes up to 192'000 such cy-cles using the powerful APE computers at DESY-IfH. The high-quality random number generatorof refs.[10,11] was used in these simulations andcMwas set equal to 1:1 in all cases. For compar-ision we did also run a Hybrid Monte Carlo [12]program for the same lattice on the CRAY-YMPat HLRZ. Here we generated 2610 trajectories ofwhich 1610 were used for \measurement".The probability distribution of the lowest ei-genvalue 0of Q2is shown in g. 2. Note that thedistribution should be the same for all algorithmswhich simulate the eective gauge eld distribu-tion Pe[U ]. The fact that there is an only small(and probably not signicant) dierence betweenthe two curves in g. 2 thus is a rst indicationthat at this level of statistics one may be able todo with n = 60 and " = 0:0015.To further investigate this question we havemade runs with " = 0:0004; 0:0007;0:0015 andvalues of n ranging from 20 to 100. For " =0:0004 all eigenvalues  are always in the interval"    1, while for " = 0:0007 there are occa-sionally a few eigenvalues below ". It should beemphasized that it is not unreasonable to set " toa value as large as 0:0015, since P () continues tobe a good approximation to 1= for eigenvalues substantially smaller than " (cf. sect. 7).In g. 3 we present our results for the plaquetteexpectation value. The data are plotted againstthe parameter = max"s1jR(s)j (22)which indicates how well the polynomialP (s) ap-proximates 1=s in the interval containing all buta few eigenvalues of Q2. At xed " and increas-ing n,  is monotonically decreasing and eventu-ally vanishes exponentially [eq.(17)]. In particu-lar, the leftmost data points in g. 3 correspondto the larger values of n.The plot shows that there is a signicant depen-dence of the average plaquette on n and " when is greater than 5% or so. Below this level thedata are consistent with a constant value and also5Table 1Compilation of simulation results with  < 0:05n "  hplaquettei h0i mm60 0:0015 0:0177 0:5802(6) 0:00110(3) 1:173(16) 1:262(17)80 0:0007 0:0275 0:5806(8) 0:00107(5) 1:186(21) 1:257(22)80 0:0015 0:0038 0:5801(8) 0:00109(4) 1:176(15) 1:263(20)100 0:0004 0:0352 0:5805(9) 0:00115(3) 1:211(20) 1:300(20)100 0:0007 0:0095 0:5814(10) 0:00115(2) 1:204(11) 1:290(13)100 0:0015 0:0008 0:5799(7) 0:00114(4) 1:208(16) 1:290(15)HMC - - 0:5796(6) 0:00117(2) - -.001 .01 .1 1.57.575.58.585Figure 3. Plaquette expectation value plotted asa function of the t accuracy  [eq.(22)]. Opencircles, full circles and crosses correspond to dif-ferent values of " (0.0004, 0.0007, 0.0015). Theline is obtained by tting the data points with < 0:05 by a constant.with the result obtained from the Hybrid MonteCarlo (cf. table 1).Similiar conclusions are reached when consid-ering the expectation value of the lowest eigen-value 0of Q2(g. 4) and the masses of the and  mesons (gs. 5,6). The latter were denedthroughcoshm =C(T=2 + 1)C(T=2); (23)where C(t) is the appropriate correlation functionof local operators at time t and T denotes the time.001 .01 .1 1.0006.0008.001.0012.0014Figure 4. Same as g. 3, but showing the expec-tation value of the lowest eigenvalue 0of Q2.size of the lattice (i.e. T = 12 in the present case).These masses are, therefore, not quite the sameas the true masses dened through the eigenval-ues of the transfer matrix. To study the proper-ties of the new simulation algorithm the deni-tion (23) is, however, equally useful and perhapseven better suited, since extrapolation errors areavoided. The correlation functions C(t) were ob-tained in the conventional way, i.e. by computingquark propagators for a representative sample ofgauge eld congurations.7. CONCLUDING REMARKSThe simulations that we have performed so farshow that the new algorithm works out in the6.001 .01 .1 1.811.21.4Figure 5. Results for the pion massm(cf. g. 3)..001 .01 .1 1.811.21.4Figure 6. Results for the  meson mass m(cf. g. 3).sense that correct and accurate results can be ob-tained with a reasonable number of boson elds.Some further consolidation by studies on largerlattices and at smaller quark masses is of courseneeded. One would also like to better understandhow to choose the parameters n and " in anygiven instance. For the lattice considered aboveone may take n = 60 and " = 0:0015 to obtaincorrect results at the level of statistical precisiongiven in table 1.It is somewhat surprising that the system-atic eects (which are otherwise clearly visible ings. 3{6) are apparently negligible for t accura-cies  as large as a few percent. Perhaps this isdue to a coherent cancellation of errors when tak-ing the product over all eigenvalues of the quarkmatrix (as already discussed in sect. 3).Another interesting observation is that when is smaller than about 5% the \measured" quan-tities did not seem to depend on the values of "we have considered. This suggests that the maineect of the quark determinant comes from thehigh-energy modes and that correlations betweenthe low-lying eigenvalues in quark loops and theobserved quantities are small.Before further extensive tests are performed,we would like to install a number of improvementsin the present version of the simulation program.a. It is well-known that the condition numberof the quark matrix can be signicantly reducedby taking advantage of its special form with re-spect to the even and odd sublattices. Even-oddpreconditioning can be incorporated in a ratherstraightforward manner in the transformation tothe bosonic theory. As a result the same level ofaccuracy should be reached with only about halfof the boson elds required without precondition-ing.b. As discussed in ref.[2] unexpectedly long auto-correlation times are observed when the bosonictheory is simulated following the lines of sect. 4.The origin of this eect is now understood and anidea has been put forward of how to accelerate thesimulation algorithm [13].c. Let A[U ] be any observable depending on thegauge eld U and < A > its expectation valuewith respect to the true QCD distribution Pe[U ].One may then show that< A >=< XA >b= < X >b; (24)where < : : : >bdenotes an expectation value inthe bosonic theory with action Sb[U; ] andX[U ] = det1 +R(Q2): (25)Approximating < A > by < A >b(as we didin the test runs) thus amounts to neglecting thecorrelations between A and X. These are guar-anteed to be small when n is large, provided none7of the eigenvalues  of Q2are substantially below".In general we haveX[U ] =Y<"[1 + R()] [1 + O()] ; (26)and we may thus attempt to improve the con-vergence in the large n limit by including X inthe \measurement" process, using eqs.(24) and(26). This may be particularly worthwhile inthose instances, where there are just a few low-lying eigenvalues. In any case, one may in thisway obtain further insight into the importance ofthese eigenvalues for the quantities of interest.REFERENCES1. M. Luscher, Nucl. Phys. B418 (1994) 6372. B. Jegerlehner, Study of a new simulation al-gorithm for dynamical quarks on the APE-100 parallel computer, talk given at the In-ternational Symposium on Lattice Field The-ory, Sept. 27{Oct. 1, 1994, Bielefeld, to bepublished3. W. H. Press, S. A. Teukolsky, W. T. Vetter-ling and B. P. Flannery, Numerical recipes inFORTRAN, 2nd ed. (Cambridge UniversityPress, Cambridge, 1992)4. M. Geradin, Journal of Sound and Vibration19 (1971) 3195. I. Fried, Journal of Sound and Vibration 20(1972) 3336. A. Meyer, Modern algorithms for largesparse eigenvalue problems (Akademie Ver-lag, Berlin, 1987)7. M. Luscher and H. Simma, Accurate summa-tion methods for the APE-100 parallel com-puter, le /pub/quadrics/doc/sum.ps onftp.ifh.de8. K. Fabricius and O. Haan, Phys. Lett. 143B(1984) 4599. A. D. Kennedy and B. J. Pendleton, Phys.Lett. 156B (1985) 39310. M. Luscher, Comp. Phys. Comm. 79 (1994)10011. M. Luscher, A random number generator forthe APE-100 parallel computer,le /pub/quadrics/doc/ranfloat.ps onftp.ifh.de12. S. Duane, A. D. Kennedy, B. J. Pendletonand D. Roweth, Phys. Lett. B195 (1987) 21613. B. Jegerlehner, work in progress

A previously introduced multi-boson technique for the simulation of QCD with dynamical quarks is described and some results of first test runs on a $6^3\times12$ lattice with Wilson quarks and gauge group $SU(2)$ are reported

DESY

A New simulation algorithm for lattice QCD with dynamical quarks

A new simulation algorithm for lattice QCD with dynamical quarks

Abstract

Similar works

Full text

Available Versions

Crossref

CERN Document Server

DESY

DESY