The standard hybrid Monte Carlo algorithm is known to simulate even flavors
QCD only. Simulations of odd flavors QCD, however, can be also 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: 9 pages, 8 figures, Proceedings of the International Workshop on
  Nonperturbative Methods and Lattice QCD, Guangzhou, Chin

de Forcrand, Philippe

Takaishi, Tetsuya

English

arXiv

The standard hybrid Monte Carlo algorithm is known to simulate even flavors QCD only. Simulations of odd flavors QCD, however, can be also 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.The standard hybrid Monte Carlo algorithm is known to simulate even flavors QCD only. Simulations of odd flavors QCD, however, can be also 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

CERN Document Server

arXiv:hep-lat/0009024v2  28 Sep 2000Odd-flavor Simulations by the Hybrid Monte Carlo aTetsuya TakaishiHiroshima University of Economics, Hiroshima, 731-0192, JAPAN.E-mail:tt-taka@hue.ac.jpPhilippe de ForcrandCERN, Theory Division, CH-1211 Geneva 23, SwitzerlandETH-Ho¨nggerberg, CH-8093 Zu¨rich, SwitzerlandE-mail:forcrand@itp.phys.ethz.chThe standard hybrid Monte Carlo algorithm is known to simulate even flavorsQCD only. Simulations of odd flavors QCD, however, can be also performed in theframework of the hybrid Monte Carlo algorithm where the inverse of the fermionmatrix is approximated by a polynomial. In this exploratory study we performthree flavors QCD simulations. We make a comparison of the hybrid Monte Carloalgorithm and the R-algorithm which also simulates odd flavors systems but hasstep-size errors. We find that results from our hybrid Monte Carlo algorithm arein agreement with those from the R-algorithm obtained at very small step-size.1 IntroductionRecent lattice QCD simulations include effects of dynamical fermions. Thestandard established algorithm to simulate dynamical QCD is the HybridMonte Carlo (HMC) algorithm 1 but it is limited to simulations of an evennumber of degenerate flavors. It would be desirable to simulate lattice QCDwith three flavors since there exist three light quarks, i.e. u,d,s quarks, inthe real world. Simulations with an odd number of flavors have been per-formed using the R-algorithm 2. This algorithm, however, is not exact: itcauses systematic errors of order ∆τ2, where ∆τ is the step-size of the Molec-ular Dynamics evolution. A careful extrapolation to zero step-size is thereforeneeded to obtain exact results. Nevertheless, it is common practice to foregothis extrapolation and to perform simulations with a single step-size chosensmall enough that the expected systematic errors are smaller than the statisti-cal ones. We want to point out that there is an alternative to the R-algorithm,which gives arbitrarily accurate results without any extrapolation3.Some time ago, a local algorithm, the so-called ”Multiboson algorithm”,was proposed by Lu¨scher 4, in which the inverse of the fermion matrix is ap-proximated by a suitable Chebyshev polynomial. Originally he proposed it fortwo flavors QCD. Boric¸i and de Forcrand 5 noticed that the determinant of aaTalk presented by T.Takaishi1fermion matrix can be written in a manifestly positive way using a polynomialapproximation, so that one can simulate odd flavors QCD with the multibosonmethod. Indeed, using this method, one flavor QCD was simulated successfully6. The same polynomial approximation can be applied for the HMC algorithm7, which means that the HMC algorithm also has the possibility of simulatingodd flavors. Here we give the formulation for simulating odd flavors with theHMC algorithm and perform three flavors simulations. Then we compare ourresults with those of the R-algorithm.2 Formulation2.1 nf=2The lattice QCD partition function with nf degenerate quark flavors is givenbyZ =∫dU det(D)nf exp(−Sgauge), (1)where D is the fermion matrix and in this study we use Wilson fermions. Fornf = 2, the partition function isZ =∫dU det(D)2 exp(−Sgauge). (2)In the formulation of the HMC algorithm, det(D)2 is treated asdetD2 ∼∫dφ†dφ exp(−φ†D†−1D†φ), (3)where the γ5 hermiticity of the fermion matrix D, i.e. D = γ5D†γ5, is used.Introducing momenta P conjugate to the link variables U , the partitionfunction is rewritten asZ =∫dUdP exp(−H), (4)where the Hamiltonian H is defined byH =12P 2 + Sgauge + φ†D†−1D−1φ. (5)This Hamiltonian is used for the Molecular Dynamics (MD) simulation ofthe standard HMC algorithm. Eq.(5) has a computational difficulty in MDsimulations since one must solve x = D−1φ type equations which in general2take a large amount of computational time for a large fermion matrix and/orfor a small quark mass.Following Lu¨scher4, the inverse ofD can be approximated by a polynomial:1/D ≈ Pn(D) ≡n∏k=1(D − Zk) (6)where Zk are roots of the polynomial Pn(D):Zk = 1− exp(i2pikn+ 1). (7)The rate of convergence of the approximation depends on the quark mass (SeeSec.3).Replacing D−1 in eq.(5) by Pn(D) we obtain an approximate Hamiltonian,Hn =12P 2 + Sgauge + φ†Pn(D)†Pn(D)φ. (8)An advantage of using Hn is that no solver calculation is required in the MDevolution. Instead, one needs n multiplications by the matrix D. OriginallyHn was introduced to reduce computational work. Indeed, it was shown thatHn can provide some gain over the standard HMC algorithm7.Hn does introduce some systematic errors from the polynomial approxi-mation. For the nf=even case, however, these errors are easily corrected atthe Metropolis step by using the exact Hamiltonian of eq.(5)7.The domain of convergence of Pn(D) is bounded by a circle centered at(1,0) which goes through the origin. If all eigenvalues of D fall inside thisdomain, Pn(D) converges exponentially. Otherwise, Pn(D) does not converge,which may happen for some exceptional configurations. Our algorithm rejectsthese configurations at the Metropolis step. This domain of convergence can bechanged by adopting another approximating polynomial. However, the originmust be excluded. Together with connectedness and conjugate symmetry of thespectrum, this implies that the real negative axis is always excluded from thedomain of convergence for any polynomial. Configurations with real negativeDirac eigenvalues will be rejected by our polynomial algorithm.2.2 nf = 1In this case, we have to consider detD. detD can not be expressed in a mani-festly positive manner using the same treatment of eq.(3). Thus the standardHMC algorithm can not handle nf = 1 or nf=odd simulations.3The multiboson algorithm was originally developed for a simulation ofnf=2 QCD4. After invention of the multiboson algorithm, Boric¸i and de For-crand 5 noticed that a single detD can be treated in a manifestly positive wayand an nf = 1 multiboson simulation was performed to study thermodynamicsof nf = 1 QCD6.As before, the inverse of the fermion matrix D, using a polynomial ofdegree 2n, is approximated as 4,51/D ≈2n∏k=1(D − Zk), (9)where Zk = 1− exp(i 2pik/(2n+ 1)). Noticing that the Zk’s come in complexconjugate pairs, eq.(9) is rewritten as1/D ≈n∏k=1(D − Z¯k)(D − Zk). (10)Using the γ5 hermiticity of the fermion matrix, we find that det(D − Z¯k) =det(D − Zk)†. Thus the determinant of D is written asdet(D) ∼ det(T †n(D)Tn(D))−1, (11)where Tn(D) ≡∏nk=1(D−Zk). Using an integral form of the determinant, weobtaindet(D) ∼∫dφ†dφ exp(−φ†T †n(D)Tn(D)φ). (12)The term φ†T †n(D)Tn(D)φ is manifestly positive. Then we may define theHamiltonian of nf = 1 QCD asH =12P 2 + Sgauge + φ†T †n(D)Tn(D)φ. (13)With this Hamiltonian there is no difficulty to perform HMC algorithm. Thedomain of convergence of the approximation eq.(12) is the same as for nf = 2.Exceptional configurations for which eigenvalues fall outside this domain willbe rejected at the Metropolis step.2.3 nf = 2 + 1The partition function of nf = 2 + 1 QCD is given byZ =∫dU det D˜2 detD exp(−Sgauge), (14)40 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: (left): Xn versus degree n. (right):|Xn −Xexact| versus degree n.where the notations D˜ and D are introduced to distinguish the two differentquark masses. Using eq.(3) for det D˜2 and eq.(12) for detD,det D˜2 detD ∼∫dφ˜†dφ˜dφ†dφ exp(−φ˜†D˜†−1D˜−1φ˜− φ†T †n(D)Tn(D)φ). (15)We define nf = 2 + 1 Hamiltonian byH =12P 2 + Sg + φ˜†D˜†−1D˜−1φ˜+ φ†T †n(D)Tn(D)φ. (16)Two remarks are in order: (i) the degree n of the approximating poly-nomial may be different during the Molecular Dynamics trajectory and forthe Metropolis acceptance test; the former can be made arbitrarily small andtuned for maximum efficiency, while the latter should be taken very large toenforce the correct measure; (ii) the two bosonic fields φ and φ˜ could be re-placed by a single one, with action φ†T †n(D)D˜†−1D˜−1Tn(D)φ. For simplicity,in this exploratory study we use two distinct bosonic fields and a single degreefor the approximating polynomial.3 Convergence3.1 nf = 2In order to see the rate of convergence of Pn(D), we calculate the quan-tity Xn = φ†P †n(D)Pn(D)φ. In the limit n → ∞, Xn goes to Xexact ≡50 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: (left): Xn versus degree n. (right):|Xn −Xmax| versus degree n.φ†D†−1D−1φ. First, we choose Xexact = η†η where η is a random gaus-sian vector. Then the vector φ is set to φ ≡ Dη. The accuracy of Xn ismeasured by the difference between Xn and Xexact. We use a random gaugeconfiguration for this analysis. Figure 1:(left) shows Xn versus the degree nfor different quark masses. Here the same η is used for each calculation of Xn.Xn converges to one value as n increases, but at high degree n, Xn diverges,which can be understood due to the rounding errors of our computer, wherecalculations are performed with 64-bit accuracy.Figure 1:(right) shows the accuracy of Xn by |Xn −Xexact|. Exponentialconvergence is seen for each quark mass, but the rate of convergence is slowfor small quark masses.3.2 nf = 1We do the same analysis as for nf = 2, but for nf = 1, the value ofXexact is notknown. So we 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 the limit of n→∞.Figure 2:(left) shows Xn as a function of degree n. Xn seems to convergesto a certain value when the degree n increases. At high degree n, Xn divergesas in the case of nf = 2.To see the rate of convergence, we calculate |Xn −Xmax| where Xmax isdefined by Xmax = Xm, m >> n. Due to the rounding errors, we can not6−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: (left): Plaquette of nf = 3 flavor QCD on an 82 × 10 × 4 lattice at β = 5.0 andat κ = 0.130 as a function of degree n. (right): Real part of Polyakov loop.take very large m. We take a maximum number m where the rounding errorsstill do not appear. Figure 2:(right) shows |Xn−Xmax| as a function of degreen. The dips seen in the figure are just due to the fact that at those pointsXn = Xmax ≡ Xm. The convergence seems to be exponential, but the rate ofconvergence is slow for small quark masses as in the nf = 2 case.4 SimulationsWe perform simulations of three flavors QCD on an 82 × 10 × 4 lattice atβ = 5.0 with κ = 0.130 and 0.160. We measure the plaquette and Polyakovloop varying the degree n and compare them with those from the R-algorithmobtained with a step-size ∆τ = 0.018. Figures 3 and 4:(left) show the plaquetteas a function of n at κ = 0.130 and 0.160, respectively. Except for very smalln, the results from the HMC algorithm agree with those from the R-algorithmwithin statistical errors. Results of the Polyakov loop are shown in figures 3and 4:(right). Except for a small discrepancy seen in figure 3, the results fromthe HMC algorithm are in agreement with those from the R-algorithm. Notethat convergence is not monotonic in n.5 ConclusionsWe formulated an odd-flavor HMC algorithm using a polynomial approxima-tion. Simulations of three flavors QCD were performed. We found that theplaquette values are consistent with those from the R-algorithm at very small7−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: (left): Plaquette of nf = 3 flavor QCD on an 82 × 10 × 4 lattice at β = 5.0 andat κ = 0.160 as a function of degree n. (right): Real part of Polyakov loop.step-size. In principle the HMC algorithm is able to simulate any flavors ofQCD, with arbitrary accuracy and without extrapolation [as long as all Diraceigenvalues are not real negative]. However the rounding errors should be un-der control when we use a large lattice or/and small quark masses where onemay need a polynomial of high degree n to achieve sufficient approximation.AcknowledgementsThis work is partially supported by the Grant-in-Aid for Scientific Researchby Monbusho, Japan(No.11740159).References1. S.Duane, A.Kennedy, B.Pendleton and D.Roweth, Phys. Lett. B 195,216 (1987)2. S. Gottlieb, W.Liu, D.Toussaint, R.L.Renken and R.L.Sugar, Phys. Rev.D 35, 2531 (1987).3. For yet another approach, see T.Lippert, hep-lat/0007029, Lect.NotesComput.Sci. and Eng. 15 (2000) 166-177.4. M. Lu¨scher, Nucl. Phys. B 418, 637 (1994);B.Bunk, K.Jansen, B.Jegerlehner, M.Lu¨scher, H.Simma and R.Sommer,Nucl. Phys. B(Proc.Suppl.) 42, 49 (1995)5. A. Boric¸i and Ph. de Forcrand, Nucl. Phys. B 454, 645 (1995).86. C. Alexandrou, A.Boric¸i, A.Feo, Ph. de Forcrand, A.Galli, F.Jegerlehnerand T.Takaishi, Phys. Rev. D 60, 034504 (1999).7. Ph. de Forcrand and T. Takaishi Nucl. Phys. B(Proc.Suppl.) 53, 968(1997)8. Y. Iwasaki, K.Kanaya, S.Kaya, S.Sakai and T.Yoshie, Phys. Rev. D 54,7010 (1996).9

Odd-flavor Simulations by the Hybrid Monte Carlo

Abstract

Similar works

Full text

Available Versions

CERN Document Server

Crossref