Simulations of odd flavors QCD can be performed in the framework of the
hybrid Monte Carlo algorithm where the inverse of the fermion matrix is
approximated by a polynomial. In this exploratory study we perform three
flavors QCD simulations. We make a comparison of the hybrid Monte Carlo
algorithm and the R-algorithm which also simulates odd flavors systems but has
step-size errors. We find that results from our hybrid Monte Carlo algorithm
are in agreement with those from the R-algorithm obtained at very small
step-size.Comment: Lattice 2000 (Algorithms), 5 pages, 8 figures, LaTe

Alexandrou

Boriçi

Bunk

de Forcrand

Duane

Frezzotti

Gottlieb

Iwasaki

Lippert

Lüscher

Montvay

Philippe de Forcrand

Tetsuya Takaishi

English

arXiv

Simulations of odd flavors QCD can be performed in the framework of the hybrid Monte Carlo algorithm where the inverse of the fermion matrix is approximated by a polynomial. In this exploratory study we perform three flavors QCD simulations. We make a comparison of the hybrid Monte Carlo algorithm and the R-algorithm which also simulates odd flavors systems but has step-size errors. We find that results from our hybrid Monte Carlo algorithm are in agreement with those from the R-algorithm obtained at very small step-size.Simulations of odd flavors QCD can be performed in the framework of the hybrid Monte Carlo algorithm where the inverse of the fermion matrix is approximated by a polynomial. In this exploratory study we perform three flavors QCD simulations. We make a comparison of the hybrid Monte Carlo algorithm and the R-algorithm which also simulates odd flavors systems but has step-size errors. We find that results from our hybrid Monte Carlo algorithm are in agreement with those from the R-algorithm obtained at very small step-size

Takaishi, Tetsuya

de Forcrand, Phillippe

CERN Document Server

arXiv:hep-lat/0011003v1  1 Nov 20001Simulation of nf = 3 QCD by Hybrid Monte CarloTetsuya Takaishia∗ and Philippe de Forcrandb caHiroshima University of Economics, Hiroshima, 731-0192, JAPANbInst. fu¨r Theoretische Physik, ETH Ho¨nggerberg, CH-8093 Zu¨rich, SwitzerlandcCERN, Theory Division, CH-1211 Gene`ve 23, SwitzerlandSimulations of odd flavors QCD can be performed in the framework of the hybrid Monte Carlo algorithm wherethe inverse of the fermion matrix is approximated by a polynomial. In this exploratory study we perform threeflavors QCD simulations. We make a comparison of the hybrid Monte Carlo algorithm and the R-algorithm whichalso simulates odd flavors systems but has step-size errors. We find that results from our hybrid Monte Carloalgorithm are in agreement with those from the R-algorithm obtained at very small step-size.1. IntroductionRecent lattice QCD simulations include effectsof dynamical fermions. Due to the algorithmiclimitation of the standard Hybrid Monte Carlo(HMC) algorithm [1], those simulations are lim-ited to even numbers of degenerate flavors. In or-der to include dynamical effects correctly, simula-tions of QCD with three flavors (u,d,s quarks) aredesirable. Simulations with an odd number of fla-vors can be performed using the R-algorithm [2].This algorithm, however, is not exact: it causessystematic errors of order ∆τ2, where ∆τ is thestep-size of the Molecular Dynamics evolution. Acareful extrapolation to zero step-size is thereforeneeded to obtain exact results. Nevertheless, itis common practice to forego this extrapolationand to perform simulations with a single step-sizechosen small enough that the expected systematicerrors are smaller than the statistical ones. Wewant to point out that there is an alternative tothe R-algorithm, which gives arbitrarily accurateresults without any extrapolation[3].Lu¨scher proposed a local algorithm, the so-called ”Multiboson algorithm” [4], in which theinverse of the fermion matrix is approximatedby a suitable Chebyshev polynomial. Originallyhe proposed it for two flavors QCD. Boric¸i andde Forcrand [5] noticed that the determinant∗Presented by T. Takaishi at LAT00.of a fermion matrix can be written in a mani-festly positive way using a polynomial approxi-mation, so that one can simulate odd flavors QCDwith the multiboson method. Indeed, using thismethod, one flavor QCD was simulated success-fully [6]. The same polynomial approximationcan be applied for the HMC [7]. Actually, inthe development stage of Ref.[6], one flavor QCDwas also simulated by HMC and it was confirmedthat the two algorithmically different methods —multiboson and HMC — give the same plaque-tte value [8]. Here we give the formulation ofthe HMC algorithm with odd flavors and performnf = 3 QCD simulations. Then we compare ourresults with those of the R-algorithm.2. Formulation2.1. nf = 2The application of Lu¨scher’s idea [4] to nf = 2QCD HMC was first made by the authors ofRef.[7], and later by [9]. The lattice QCD par-tition function with nf = 2 degenerate quark fla-vors is given byZ =∫dU detD2 exp(−Sgauge), (1)where D is the fermion matrix and in this studywe use Wilson fermions. In the formulation of the2HMC algorithm, detD2 is treated asdetD2 ∼∫dφ†dφ exp(−φ†D†−1D−1φ), (2)where the γ5 hermiticity of the fermion matrix D,i.e. D = γ5D†γ5, is used.Introducing momenta P conjugate to the linkvariables U , the partition function is rewritten asZ =∫dUdP exp(−H), (3)where the Hamiltonian H is defined byH =12P 2 + Sgauge + φ†D†−1D−1φ. (4)This Hamiltonian is used for the Molecular Dy-namics (MD) simulation of the standard HMC al-gorithm. Eq.(4) has a computational difficulty inMD simulations since one must solve x = D−1φtype equations which in general take a largeamount of computational time for a large fermionmatrix and/or for a small quark mass.Following Lu¨scher [4], the inverse of D can beapproximated by a polynomial:1/D ≈ Pn(D) ≡n∏k=1(D − Zk), (5)where Zk are the roots of the polynomial Pn(D).We choose Zk = 1− exp(i 2pik/(n+ 1)).Replacing D−1 in eq.(4) by Pn(D) we obtainan approximate Hamiltonian,Hn =12P 2 + Sgauge + φ†Pn(D)†Pn(D)φ. (6)An advantage of using Hn is that no solver calcu-lation is required in the MD evolution. Instead,one needs n multiplications by the matrix D.Hn does introduce some systematic errors fromthe polynomial approximation. For the nf=evencase, however, these errors are easily correctedby using the exact Hamiltonian of eq.(4) at theMetropolis step[7]. This guarantees that configu-rations will be distributed according to the exactmeasure ∝ detD2, for any polynomial Pn(D).However, the domain of convergence of Pn(D)is bounded by a circle centered at (1, 0) whichgoes through the origin. If all eigenvalues of Dfall inside this domain, Pn(D) converges expo-nentially. Otherwise, Pn(D) does not converge,which may happen for some exceptional configu-rations. Our algorithm will tend to reject theseconfigurations at the Metropolis step, leading toextremely long autocorrelation times. This do-main of convergence can be changed by adoptinganother approximating polynomial. However, theorigin must be excluded. Together with connect-edness and conjugate symmetry of the spectrum,this implies that the real negative axis is alwaysexcluded from the domain of convergence for anypolynomial. Configurations with real negativeDirac eigenvalues will be rejected by our poly-nomial algorithm.2.2. nf = 1After invention of the multiboson algorithm,Boric¸i and de Forcrand [5] noticed that a singledetD can be treated in a manifestly positive wayand an nf = 1 multiboson simulation was per-formed to study thermodynamics of nf = 1 QCD[6].As before, the inverse of the fermion matrix D,using a polynomial of degree 2n, is approximatedas [4,5]1/D ≈2n∏k=1(D − Zk), (7)where Zk = 1 − exp(i 2pik/(2n + 1)). Noticingthat the Zk’s come in complex conjugate pairs,eq.(7) is rewritten as1/D ≈n∏k=1(D − Z¯k)(D − Zk). (8)Using the γ5 hermiticity of the fermion matrix,we find that det(D − Z¯k) = det(D − Zk)†. Thusthe determinant of D is written asdet(D) ∼ det(T †n(D)Tn(D))−1, (9)where Tn(D) ≡∏nk=1(D − Zk), and then we ob-taindet(D) ∼∫dφ†dφ exp(−φ†T †n(D)Tn(D)φ). (10)The term φ†T †n(D)Tn(D)φ is manifestly positive.Then we may define the Hamiltonian of nf = 13QCD asH =12P 2 + Sgauge + φ†T †n(D)Tn(D)φ. (11)With this Hamiltonian there is no difficulty toperform HMC algorithm. To improve efficiencyand accuracy, one may use a polynomial of lowerdegree n during the Molecular Dynamics trajec-tory, and a much higher degree m ≫ n for theMetropolis step [6,10]. The domain of conver-gence of the approximation eq.(10) is the sameas for nf = 2. Exceptional configurations forwhich eigenvalues fall outside this domain willlikewise be rejected at the Metropolis step. Afurther difficulty is that the sampled measureis ∝ det(T †m(D)Tm(D))−1, which for exceptionalconfigurations differs from the desired detD, in-creasingly so with m.0 50 100 150n  ( of  Pn(D) )4860048700488004890049000X nNf=2   84   Xn κ=0.160 κ=0.180  κ=0.200 κ=0.215 κ=0.2300 50 100 150n  ( of  Pn(D) )10−1010−810−610−410−21001021041061081010 |X n−X exact|Nf=2   84  |Xn−Xexact| κ=0.160 κ=0.180  κ=0.200 κ=0.215 κ=0.230Figure 1. (top): Xn versus degree n;(bottom):|Xn −Xexact| versus degree n.2.3. nf = 2 + 1The partition function of nf = 2 + 1 QCD isgiven byZ =∫dU det D˜2 detD exp(−Sgauge), (12)where the notations D˜ and D are introduced todistinguish the two different quark masses. Usingeq.(2) for det D˜2 and eq.(10) for detD,det D˜2 detD∼∫dφ˜†dφ˜dφ†dφ exp(−φ˜†D˜†−1D˜−1φ˜−HT ),(13)where HT ≡ φ†T †n(D)Tn(D)φ. We define nf =2 + 1 Hamiltonian byH =12P 2 + Sg + φ˜†D˜†−1D˜−1φ˜+HT . (14)Two remarks are in order: (i) as for nf = 1,one could use during the MD trajectory a poly-nomial of lower degree than for the Metropo-lis step; (ii) the two bosonic fields φ and φ˜could be replaced by a single one, with actionφ†T †n(D)D˜†−1D˜−1Tn(D)φ. For simplicity, in thisexploratory study we use two distinct bosonicfields and a single approximating polynomial.3. Convergence3.1. nf = 2To see the rate of convergence of Pn(D), wecalculate the quantity Xn = φ†P †n(D)Pn(D)φ.In the limit n → ∞, Xn goes to Xexact ≡φ†D†−1D−1φ. First, we choose Xexact = η†ηwhere η is a random gaussian vector. Then thevector φ is set to φ ≡ Dη. The accuracy ofXn is measured by the difference between Xnand Xexact. We use a random gauge configu-ration. Figure 1:(top) shows Xn versus the de-gree n for different quark masses. Here the sameη is used for each calculation of Xn. Xn con-verges to one value as n increases, but at highdegree n, Xn diverges, which can be understooddue to the rounding errors of our computer, wherecalculations are performed with 64-bit accuracy.Figure 1:(bottom) shows the accuracy of Xn by|Xn − Xexact|. Exponential convergence is seenfor each quark mass, but the rate of convergenceis slow for small quark masses.40 50 100 150n  ( of  Tn(D) )3000400050006000 X nNf=1    44    Xnκ=0.160κ=0.180κ=0.200κ=0.215κ=0.2300 50 100 150n  ( of  Tn(D) )10−510−410−310−210−1100101102103104 |X n−X max|Nf=1    44     |Xn−Xmax|κ=0.160κ=0.180κ=0.200κ=0.215κ=0.230Figure 2. (top): Xn versus degree n;(bottom):|Xn −Xmax| versus degree n.3.2. nf = 1We do the same analysis as for nf = 2, butfor nf = 1, the value of Xexact is not known. Sowe calculate the quantity Xn = φ†T †n(D)Tn(D)φ,where the vector φ is a gaussian random vector,and we use a random gauge configuration. Weassume that Xn goes to a certain value in thelimit of n → ∞. Figure 2:(top) shows Xn asa function of degree n. Xn seems to convergesto a certain value when the degree n increases.At high degree n, Xn diverges as in the case ofnf = 2.To see the rate of convergence, we calculate|Xn −Xmax| where Xmax is defined by Xmax =Xm, m ≫ n. Due to the rounding errors, wecan not take m very large. We take a maximum−20 0 20 40 60 80 100n (  of Tn(D) )0.4090.40950.410.41050.4110.41150.412Plaq.82x10x4  β=5.0  κ=0.130  Nf=3HMCR−alg.  [Iwasaki]−20 0 20 40 60 80 100n (  of Tn(D) )0.010.0110.0120.0130.0140.0150.0160.0170.0180.019Re(Pol)82x10x4  β=5.0  κ=0.130  Nf=3HMCR−alg.  [Iwasaki]Figure 3. (top): Plaquette of nf = 3 flavorQCD on an 82 × 10 × 4 lattice at β = 5.0 andat κ = 0.130 as a function of degree n; (bottom):Real part of Polyakov loop.number m where the rounding errors still do notappear. Figure 2:(bottom) shows |Xn−Xmax| asa function of degree n. The dips seen in the fig-ure are just due to the fact that at those pointsXn = Xmax ≡ Xm. The convergence seems to beexponential, but the rate of convergence is slowfor small quark masses as in the nf = 2 case.4. SimulationsSimulations of three flavors QCD are performedon an 82 × 10 × 4 lattice at β = 5.0 with κ =0.130 and 0.160. We measure the plaquette andPolyakov loop varying the degree n and comparethem with those from the R-algorithm obtained5with a step-size ∆τ = 0.01 [11]. Figures 3 and4:(top) show the plaquette as a function of n atκ = 0.130 and 0.160, respectively. Except forvery small n, the results from the HMC algorithmagree with those from the R-algorithm within sta-tistical errors. Results of the Polyakov loop areshown in Figures 3 and 4:(bottom). Except fora small discrepancy seen in Figure 3, the resultsfrom the HMC algorithm are in agreement withthose from the R-algorithm. Note that conver-gence is not monotonic in n.5. ConclusionsWe formulated an odd-flavor HMC algorithmusing a polynomial approximation. Simulationsof three flavors QCD were performed. We foundthat the plaquette values are consistent withthose from the R-algorithm at very small step-size. In principle the HMC algorithm is able tosimulate any flavors of QCD, with arbitrary ac-curacy and without extrapolation [as long as allDirac eigenvalues are not real negative]. Howeverthe rounding errors should be under control whenwe use a large lattice or/and small quark masseswhere one may need a polynomial of high degreen to achieve sufficient approximation.This work was partially supported by the Min-sitry of Education, Science, Sports and Culture,Grant-in-Aid, No.11740159.REFERENCES1. S. Duane, A. Kennedy, B. Pendleton andD. Roweth, Phys. Lett. B 195 (1987) 216.2. S. Gottlieb, W. Liu, D. Toussaint,R.L. Renken and R.L. Sugar, Phys. Rev. D35 (1987) 2531.3. For yet another approach, see T. Lip-pert, hep-lat/0007029, Lect.Notes Com-put.Sci. and Eng. 15 (2000) 166-177.4. M. Lu¨scher, Nucl. Phys. B 418 (1994) 637;B. Bunk, K. Jansen, B. Jegerlehner,M. Lu¨scher, H. Simma and R. Sommer, Nucl.Phys. B(Proc.Suppl.) 42 (1995) 49.5. A. Boric¸i and Ph. de Forcrand, Nucl. Phys. B454 (1995) 645.6. C. Alexandrou, A. Boric¸i, A. Feo, Ph. de−20 0 20 40 60 80 100n (  of Tn(D) )0.4300.4320.4340.4360.4380.4400.442Plaq.82x10x4  β=5.0  κ=0.160  Nf=3HMC  ( this work )R−algorithm [Iwasaki]−20 0 20 40 60 80 100n (  of Tn(D) )0.0470.0490.0510.0530.0550.057Re(Pol)82x10x4  β=5.0  κ=0.160  Nf=3HMC ( this work )R−algorithm [Iwasaki]Figure 4. (top): Plaquette of nf = 3 flavorQCD on an 82 × 10 × 4 lattice at β = 5.0 andat κ = 0.160 as a function of degree n; (bottom):Real part of Polyakov loop.Forcrand, A. Galli, F. Jegerlehner andT. Takaishi, Nucl. Phys. B(Proc.Suppl.) 53(1997) 435; Nucl. Phys. B(Proc.Suppl.) 63(1998) 406; Phys. Rev. D 60 (1999) 034504.7. Ph. de Forcrand and T. Takaishi Nucl. Phys.B(Proc.Suppl.) 53 (1997) 968.8. C. Alexandrou, A. Boric¸i, A. Feo, Ph. deForcrand, A. Galli, F. Jegerlehner andT. Takaishi, unpublished.9. R. Frezzotti and K. Jansen, Phys. Lett. B 402(1997) 328.10. I. Montvay, Nucl. Phys. B 466 (1996) 259.11. Y. Iwasaki, K. Kanaya, S. Kaya, S. Sakai andT. Yoshie, Phys. Rev. D 54 (1996) 7010.

Simulation of $n_{f}$ = 3 QCD by Hybrid Monte Carlo

Crossref

Simulation of nf = 3 QCD by Hybrid Monte Carlo

Simulation of n_f =3 QCD by Hybrid Monte Carlo

Simulation of n_f =3 QCD by Hybrid Monte Carlo

Abstract

Similar works

Full text

Available Versions

CERN Document Server

Crossref