2+1 Flavor QCD simulated in the epsilon-regime in different topological sectors
P. Hasenfratz, D. Hierl, V. Maillart, F. Niedermayer, A. Schäfer, C. Weiermann, M. Weingart
aa r X i v : . [ h e p - l a t ] J un ǫ -regime indifferent topological sectors P. Hasenfratz a , D. Hierl b , V. Maillart a , F. Niedermayer a , A. Sch¨afer b ,C. Weiermann a and M. Weingart aa Institute for Theoretical PhysicsUniversity of BernSidlerstrasse 5, CH-3012 Bern, Switzerland b Institute for Theoretical Physics, University of Regensburg, D-93040Regensburg, Germany
Abstract
We generated configurations with the parametrized fixed-point Diracoperator D FP on a (1 . box at a lattice spacing a = 0 .
13 fm. Wecompare the distributions of the three lowest k = 1 , , ν = 0 , , k and ν . After including the finite size correction from one-loopchiral perturbation theory we obtained for the chiral condensate in theMS scheme [Σ(2 GeV)] / = 0 . pontaneous chiral symmetry breaking and the related existence of lightGoldstone bosons is a basic feature of QCD. Chiral Perturbation Theory (ChPT)provides a systematic description of this physics in terms of a set of low en-ergy constants which encode the related non-perturbative features of QCD. Themethod of low energy effective Lagrangians simplifies the calculations signifi-cantly [1, 2, 3, 4] and over the years ChPT became a refined powerful technique.Gasser and Leutwyler recognized very early, decades before the numerical cal-culations could attack such problems, that these constants can also be fixedusing physical quantities which, presumably, will never be measured in realexperiments. They can be studied, however, in lattice QCD.The ǫ - regime [5, 6, 7, 8] describes physics close to the chiral limit in a boxwhose size is larger than the QCD scale. On the other hand, the size of thebox relative to the Goldstone boson correlation length must be small. Underthese conditions the Goldstone bosons, as opposed to other excitations, feelthe effect of boundaries strongly. ChPT provides a powerful and systematicway to calculate the finite size corrections. A nice additional feature of the ǫ -regime is that Random Matrix Theory (RMT) [9] makes precise predictions formicroscopic observables. RMT relates in particular the distribution of low-lyingeigenvalues of the Dirac operator in different topological sectors to the chiralcondensate.The pioneering numerical works [10, 11, 12, 13, 14, 15, 16] in the ǫ -regimein quenched QCD suggested that this regime can be an excellent tool to studylow-energy physics in QCD. The special problems of the ǫ -regime called fornew numerical procedures [14, 15, 16, 17] which were first tested also in thisapproximation. Combining these numerical developments with the renormal-ization properties of the spectral density the work [18] present a state of theart analysis for the quark condensate and the first Leutwyler-Smilga sum rulein quenched QCD. The quenched approximation, however, tends to be singularin the chiral limit and it is expected that quenching is even more problematichere than in other cases [19].Unfortunately, full QCD simulations are expensive. In addition, in the ǫ − regime the Dirac operator should have excellent chiral properties. The stan-dard choice is the overlap Dirac operator [20] with a hybrid Monte Carlo algo-rithm.The results in [21] reflect already the basic physics features of the ǫ -regime.The distribution of the low-lying eigenmodes could be fitted quite well to theRMT predictions and a reasonable value for the chiral condensate was obtained.These results are promising given the fact that the box was small (1 . m q ≥
40 MeV) and the lattice was coarse( a = 0 .
16 fm). The results in a larger box of 1 . m q ≈
20 MeV obtained in [22] were consistent with those of [21].2he first serious simulation results in the ǫ -regime have been presented by theJLQCD group recently [23]. In that work N f = 2 QCD was simulated with over-lap fermions (created with the Wilson kernel) on a lattice of size 1 . × . with a resolution a = 0 .
11 fm obtained from r . Different bare quark masses( am q ) were considered in the range 0 . , . . . , . , . ≈ ǫ -regime. Hybrid Monte Carlo with overlap fermions has problems when the topo-logical charge changes. To avoid this an action was used which prevented topol-ogy change and the whole run stayed in the Q = 0 sector. The authors observedan overall good agreement with RMT. The fermion condensate, which is theonly parameter to be fitted, was found to be Σ(2 GeV) / = 0 . D † GW + D GW = D † GW RD GW , (1)where R is a local operator and is trivial in Dirac space. The parametrized FPaction has many gauge paths and involves a special smearing with projectionto the gauge group SU(3). As a consequence, hybrid Monte Carlo algorithmscan not be used. In the work [25] we advised a partially global update withthree nested accept/reject steps which can reach small quark masses even oncoarse lattices. Several steps of this algorithm were developed in [26] usingsome earlier suggestions [27, 28, 29, 30]. All the three steps are preconditioned.Pieces of the quark determinant are switched on gradually in the order of theircomputational expenses. The largest and most fluctuating part of the deter-minant is coming from the ultraviolet modes. We reduce these fluctuations bycalculating the trace of D n FP , n = 1 , ..., ∼
100 lowest-lyingmodes are calculated and subtracted. The determinant of the reduced and sub-tracted D FP is calculated stochastically [33, 34, 35, 36]. For the determinantbreakup we generalized the mass shifting method of [31]. For the strange quarkwe performed a root operation. Since the strange quark mass is not very small,we used a polynomial expansion to approximate this root operator. Recentlywe replaced the polynomial expansion by a rational approximation [37] whichbrought a ≈
30% performance gain in this part.Due to the partial global updating procedure the algorithm, beyond a cer-tain volume, scales with V which constrains the size of lattices which can beconsidered. Due to the large number of subtracted low-lying modes, small quarkmass is not a barrier.We generated ∼ lattice us-ing the partially global algorithm discussed above. The distance between everysecond configurations in this chain is similar to that of two gauge configurationsseparated by a typical Metropolis gauge update sweep. We considered every3enth configurations from the Markov chain and performed the measurementson the remaining ∼
400 configurations. When calculating statistical errors weformed 20 bins and used the jackknife method.We fixed the lattice spacing from the Sommer parameter [38] r = 0 .
49 fmand found a = 0 . . .The Dirac operator D FP has no exact chiral symmetry due to parametriza-tion errors. As a consequence, the quark masses have an additive mass renormal-ization. The degenerate u, d and the s quark masses in the code are M ud = 0 . M s = 0 . M = 0 . m ud = 0 . m s = 0 . M A W I uu us ssM =0.0147(3) M Figure 1: The AWI mass for the three combinations of quark masses. Thelinear extrapolation to M AWI = 0 gives an additive mass renormalization M =0 . m ud were not useful either.RMT predicts the probability distribution p νk ( ξ νk ) for the k -th low-lying eigen-value ( k = 1 , , . . . ) of the Dirac operator in the topological sector ν . Denotingthe corresponding eigenvalues of the (continuum) Dirac operator by iα νk , thevariable ξ νk is related to the bare chiral condensate Σ as ξ νk = α νk Σ V . Here Vis the volume of the box. Fig. 2 shows the prediction of RMT for the cumulativedistributions R ξ νk dξp νk ( ξ ) for ν = 0 and k = 1 , ,
3. The distributions depend4n µ i = m i Σ V where m i are the quark masses. Decreasing the m ud mass bymore than a factor 10 (at fixed m s ) the cumulative distributions practicallyremain unchanged. On the other hand sending the strange quark mass to in-finity (at fixed m ud ) we land on a N f = 2 flavor theory with a visibly differentdistribution. The strange quark has a (modest) effect on our observables. ξ f =3, µ = µ =1.43, µ =12.3N f =2, µ = µ =1.43 N f =3, µ = µ =0.10, µ =12.3 Figure 2: Random Matrix Theory prediction for the cumulative distributionof ξ = ξ νk = α νk Σ V , where iα νk is the k-th eigenvalue of the continuum Diracoperator in a gauge background with topological charge ν , and µ i = m i Σ V .Here the ν = 0 results are shown, but the picture is similar for ν = 1 , R = 1 lies on the circle | λ − | = 1. In our casewhere 2 R is different from 1, it is convenient to introduce the rescaled oper-ator ˆ D GW = √ RD GW √ R , for which the 2 R factor is eliminated in the GWrelation, and the spectrum lies on the circle. (In fact, our operator 2 R for thelow-lying modes is effectively a constant close to 1 within a few percent.) Torelate the eigenvalues on the GW circle to those appearing in the RMT (or ingeneral, in continuum expressions) it is natural to use the stereographic projec-tion iα = λ − λ/ . (2)As Fig. 3 shows, the low-lying eigenvalues of D FP ( m ) (which we obtain duringour simulation) are close but not exactly on the GW circle, hence we project5hem first horizontally onto the GW circle, and make the stereographic pro-jection in the next step. This is justified by the observation that making asystematic GW improvement of D FP towards D GW the imaginary part of theeigenvalues stays practically constant – they move horizontally to the GW cir-cle. The values of α obtained by this procedure are then compared with theRMT predictions. δ r x24 Figure 3: The histogram of the deviation from the circle, δr = 1 − | λ − | , forthe eigenvalues of D FP on 8 ×
24 and 12 lattices for the complex eigenvalueswith | Im λ | < .
1, using the present mass parameters.Having the probability distributions from RMT we can calculate the ratios h ξ νk i / h ξ ν ′ k ′ i , where the factor Σ V cancels. These predictions are comparedwith the measured ratios in Fig. 4, where all the ratios in the topological sectors ν = 1 , , D FP has noexact chiral symmetry and has modes on the real axis. Occasionally, the realeigenvalue λ might be even far from zero in which case the topological interpre-tation is uncertain. We note here that using an exact Ginsparg-Wilson operator(with kernel D FP ) for measuring the topological charge is not a good solutionto this problem. The overlap just projects most of the real eigenvalues to thepoint λ = 0. Following the intuitive picture that topology is related to extendedobjects we investigated the correlation between the inverse participation ratio(IPR) of the eigenvector having a real eigenvalue λ and the value of λ . Here6 ν k/01 ν k/02 ν k/11 ν k/12 ν k/21 ν k/2202031112 132122230312 13222302 0312 1321 2223031322230203 12132223 031323 Figure 4: Ratios of expectation values h α νk i obtained from simulations com-pared to the results of RMT, where νk is indexing the k -th lowest eigenvalue inthe topological sector ν . The different symbols refer to the denominator while νk of the numerator is indicated at the data points. For example, the highestratio with value ≈ h ξ i / h ξ i .IPR is given by P x P i =1 | ψ ( λ ) i ( x ) | for a normalized eigenvector ψ ( λ ) . The IPRis inversely proportional to the size of the effective support of the eigenfunction.As Fig. 5 shows, there is a strong correlation indeed. As λ is moving awayfrom zero, the IPR increases indicating that the wave function becomes moreand more localized. Fig. 6 demonstrates that the real eigenvalues λ ( D FP ) arestrongly concentrated in a small region close to zero. The histogram of λ ( D FP )has a large, narrow peak followed by a long tail.Fig. 7 shows the histogram of IPR for the real eigenmodes of D FP . Thereis a strong, narrow peak which corresponds to extended real modes and a longtail corresponding to wave functions with smaller support.For the questions we study in this paper it is not necessary to use a procedurewhich leaves the relative weights of the different topological sectors intact. (Dueto the long autocorrelation time of the topological charge it would be anyhowdifficult to control these weights.) We introduce therefore a cut above whichthe real eigenvalues are discarded. This cut must be sufficiently small to allowfor a good identification of the topological charge, but large enough to allow fora reasonable number of configurations in each sectors we want to investigate.The cut should be smaller than the spectral gap (for the spectrum see Fig. 8),7 λ I P R Figure 5: The inverse participation ratio P x P i =1 | ψ ( λ ) i ( x ) | for real eigenvalues λ ( D FP ) vs. the eigenvalue. Eigenvectors with larger real part have smallerextension in agreement with the expectations.otherwise the effect of topology will be strongly reduced. We introduced a cut ofsize 0.03 which is approximately equal to the minimal gap in the ν = 1 sector.We accepted a real eigenvalue λ as a sign of topology if λ < .
03. This cutsplits our sample of 414 configurations into 328, 51 and 35 configurations forthe ν = 0 , , ν = 0 topological sector obtained using the overlap improved fixed-pointDirac operator D GW vs. D FP in the measurement. The distributions are verysimilar showing that the two operators are close to each other. The k = 1distribution from D GW has a tail at small eigenvalues (i.e. has a smaller gap than D FP ) demonstrating a systematic error one obtains when using different Diracoperators for the generation of the configurations and for the measurements:the small eigenvalues are less effectively suppressed.In Figs. 10, 11, 12 the cumulative distributions of the D FP operator arecompared with RMT for the ν = 0 , , µ i = m i Σ V and in α Σ V . We obtained the result for this bare quantity Σ / =0 . a . 8 λ Figure 6: Histogram of λ ( D FP ) for the real eigenvalues. In determining thetopological charge ν we considered only the real eigenvalues with λ < . R . For this reason the quark mass enters in a simple additiveway in the form ˆ D + m (1 − ˆ D/ D = √ RD FP √ R . Effectively theoperator 2 R behaves in the infrared like a constant close to 1. Its expectationvalue for the lowest ∼
100 eigenvectors is 1 .
05 within 1%. Using the spectrumof ˆ D (which is in fact identical to that of 2 RD FP ) the matching with RMT givesthe slightly changed result Σ / = 0 . N ) model one should interpret Σ asthe absolute value of the magnetization in the finite volume V . This differs fromthe value Σ ∞ defined in the infinite volume by a finite-size correction which canbe calculated in ChPT. In the presence of the magnetic field h the orientationof magnetization is controlled by the Boltzmann factor exp( hM cos θ ) where M is the total magnetization. This gives h M k i = M Y ′ N ( hM ) Y N ( hM ) , (3)where Y N ( z ) is related to the modified Bessel functions. Comparing this ex-pression with the result of ChPT [7] we get for N = 4 (corresponding to two9 Figure 7: Histogram of inverse participation ratio for the real eigenmodes of D FP .flavors) Σ = ρ Σ ∞ , where ρ = (cid:18) β F L (cid:19) . (4)Here the shape coefficient β takes the value 0.14046 for a symmetric box [7].The one-loop finite-size correction to the order parameter (magnetization in theexample above) has been calculated for the SU( N f ) × SU( N f ) symmetry groupin [5, 6], for the O( N ) group in [41] and up to the two-loop level in [7].Neglecting the strange quark contribution in the finite size effects we use thetwo-flavor result to correct the measured bare Σ / = 0 . / ∞ = Σ / / .
119 = 0 . . MS (2 GeV)] / = 0 . . (5)10 x24 12 Figure 8: The spectrum of D FP for 50 configurations on 8 ×
24 and 12 lattices,at the same bare mass parameters and gauge coupling. α FP D GW Figure 9: The cumulative distributions in the ν = 0 topological sector for thethree lowest complex eigenvalues as measured with D FP and with the overlapoperator D GW using the fixed-point operator as kernel.11 ξ Σ =0.291(3) GeV ν =0 Figure 10: The cumulative distribution of ξ νk = α νk Σ V for ν = 0, k = 1 , , α is obtained from λ ( D FP ) by stereographic projection described in thetext, and is rescaled using Σ / = 0 . µ = µ = m ud Σ V = 1 .
43 and µ = m s Σ V = 12 . ξ Σ =0.291(3) GeV ν =1 Figure 11: The cumulative distribution of ξ νk = α νk Σ V for ν = 1, k = 1 , , ξ Σ =0.291(3) GeV ν =2 Figure 12: The cumulative distribution of ξ νk = α νk Σ V for ν = 2, k = 1 , , cknowledgements We thank Gilberto Colangelo, Stephan D¨urr, J¨urgGasser, Anna Hasenfratz and the members of the BGR Collaboration for valu-able discussions. We also thank the LRZ in Munich and CSCS in Manno forsupport. The analysis was done on the PC clusters of ITP in Bern. This workwas supported by the Schweizerischer Nationalfonds.
APPENDIX
The bare lattice scalar density is not universal, it depends on the lattice actionand it needs renormalization. The conventional way is to express the result inMS scheme at some given scale, e.g. µ = 2 GeV.To relate the bare lattice result to MS scheme one generally uses the Rome-Southhampton RI/MOM method [42]. In this method one introduces an inter-mediate scheme RI (or RI ′ ) where some Green’s functions of the actual operator(calculated in a fixed gauge) are equated to their corresponding Born terms ata given scale p = µ .In the first step one relates the bare lattice operator non-perturbatively toits renormalized counterpart in the RI scheme O RIR ( µ ) = Z RI , lat O ( µ, a ) O latbare ( a ) . (A-1)The matching scale µ is restricted by two conditions: It should be sufficientlylarge to avoid non-perturbative effects and small enough to avoid large cut-offeffects. For the scalar density, the non-perturbative matching was done in [43],using a technique which removes a large part of O( a ) cut-off effects.In the second step one relates the renormalized operators in the RI and MSschemes in perturbation theory O MSR ( µ ) = Z MS , RI ( µ ) O RIR ( µ ) . (A-2)Obviously, the value of µ must lie in the perturbative regime. The factor Z MS , RI ( µ ) has been calculated up to NNLO [44] and NNNLO [45]. Combin-ing (A-1) and (A-2), we obtain the matching factor connecting the bare latticeand the renormalized MS results Z MS , lat O ( µ, a ) = Z MS , RI ( µ ) Z RI , lat O ( µ, a ) . (A-3)As the expansion parameter in the perturbative series is the running coupling α ( µ ) of QCD, we need to determine it for the scale µ where the matching is These schemes differ in their definition of the quark field renormalization factor Z q . α is described by the differential equation dα ( µ ) d ln µ = β ( α ) , (A-4)where the 4-loop β -function in MS scheme is given in [48]. As initial conditionwe use α (2 GeV) = 0 . .The procedure described above is carriedout at several matching scales µ ≥ Z MS , lats ( µ , a ) for the scalar density at µ = 2 GeV, Z MS , lats ( µ , a ) = Z MS , lats ( µ, a ) exp − Z α ( µ ) α ( µ ) dα γ s ( α ) β ( α ) ! . (A-5)Here, due to the relation Z s = Z − m , the anomalous dimension γ s is related tothe mass anomalous dimension by γ s = − γ m . The latter has been calculated tofour loops in [46, 47].In Fig. 13, we plot our result Z MS , lats ( µ , a ) at µ = 2 GeV vs. the matchingscale µ for the intermediate schemes RI and RI ′ . If the non-perturbative, cut-off and higher order perturbative effects were negligible, the two curves wouldcoincide and would be independent of the choice of µ . Taking the average valuesfor 2 GeV ≤ µ ≤ Z MS , lats ( µ = 2 GeV , a ) = 0 . This result is obtained by using the PDG value α (M Z ) = 0 . ± .
002 [50] and runningit from 5-flavors down to 3-flavors, across the m b and m c thresholds µ [GeV ] Z s M S , l a t ( µ , a ) with RI’with RI Figure 13: The renormalization factor Z MS , lats ( µ , a ) for the scalar density,connecting the MS and lattice schemes, at µ = 2 GeV, obtained using RI andRI ′ schemes at the intermediate step vs. the matching scale µ .16 eferences [1] S. Weinberg, Physica A (1979) 327.[2] J. Gasser and H. Leutwyler, Phys. Lett. B (1983) 321, 325.[3] J. Gasser and H. Leutwyler, Ann. Phys. (N.Y.). (1984) 142.[4] J. Gasser and H. Leutwyler, Nucl. Phys. B (1985) 465.[5] J. Gasser and H. Leutwyler, Phys. Lett. B (1987) 83; (1987) 477.[6] J. Gasser and H. Leutwyler, Nucl. Phys. B (1988) 763.[7] P. Hasenfratz and H. Leutwyler, Nucl. Phys. B (1990) 241.[8] F. C. Hansen, Nucl. Phys. B (1990) 685.[9] E. V. Shuryak and J. J. M. Verbaarschot, Nucl. Phys. A (1993) 306[arXiv:hep-th/9212088],J. J. M. Verbaarschot and T. Wettig, Ann. Rev. Nucl. Part. Sci (2000)343 [arXiv:hep-ph/0003017],P. H. Damgaard and S. M. Nishigaki, Phys. Rev. (2001) 045012[arXiv:hep-th/0006111][10] R. G. Edwards, U. M. Heller, J. E. Kiskis and R. Narayanan,Phys. Rev. Lett. (1999) 4188.[11] W. Bietenholz, K. Jansen and S. Shcheredin, JHEP (2003) 033.[12] W. Bietenholz, Th. Chiarappa, K. Jansen Kei-ichi Nagai and S. Shcheredin,JHEP (2004) 023.[13] L. Giusti, M. L¨uscher, P. Weisz and H. Wittig, JHEP (2003) 23[arXiv:hep-lat/0309189].[14] L. Giusti, P. Hernandez, M. Laine, P. Weisz and H. Wittig, JHEP (2004) 013 [arXiv:hep-lat/0402002].[15] R. G. Edwards, Nucl. Phys. Proc. Suppl. (2002) 38[arXiv:hep-lat/0111009].[16] T. DeGrand and S. Schaefer, Comput. Phys. Commun. (2004) 185[arXiv:hep-lat/0401011].[17] L. Giusti, C. Hoelbling, M. L¨uscher and H. Wittig, Comput. Phys. Com-mun. (2003) 31 [arXiv:hep-lat/0212012].[18] L. Giusti and S. Necco, JHEP (2007) 090 [arXiv:hep-lat/0702013].[19] P. H. Damgaard, Nucl. Phys. Proc. Suppl. (2004) 47[arXiv:hep-lat/0310037]. 1720] H. Neuberger, Phys. Lett. B (1998) 141 [arXiv:hep-lat/9707022]; (1998) 353 [arXiv:hep-lat/9801031].[21] T. DeGrand and S. Schaefer, Phys. Rev. D (2005) 034507[arXiv:hep-lat/0412005]; (2005) 054503 [arXiv:hep-lat/0506021]; PoSLAT2005 (2006) 140 [arXiv:hep-lat/0508025].[22] Th. DeGrand, Zh. Liu, S. Schaefer, Phys. Rev. D (2006) 094504,Erratum-ibid. D (2006) 099904 [arXiv:hep-lat/0608019];[23] H. Fukaya et al. [JLQCD Collaboration], Phys. Rev. Lett. , 172001 (2007)[arXiv:hep-lat/0702003].[24] P. Hasenfratz, S. Hauswirth, K. Holland, T. Jorg, F. Niedermayer andU. Wenger, Int. J. Mod. Phys. C (2001) 691 [arXiv:hep-lat/0003013],P. Hasenfratz, S. Hauswirth, T. Jorg, F. Niedermayer and K. Holland,Nucl. Phys.B (2002) 280 [arXiv:hep-lat/0205010],C. Gattringer et al. Bern-Graz-Regensburg Collaboration, Nucl.Phys.B (2004) 3 [arXiv:hep-lat/0307013],P. Hasenfratz, K. J. Juge and F. Niedermayer, Bern-Graz-Regensburg Col-laboration, JHEP 0412:030,2004 [arXiv:hep-lat/0411034].[25] A. Hasenfratz, P. Hasenfratz and F. Niedermayer, Phys.Rev. D72 (2005)114508 [arXiv:hep-lat/0506024].[26] A. Hasenfratz and F. Knechtli, Comput. Phys. Commun. (2002) 81[arXiv:hep-lat/0203010].[27] M. Grady, Phys. Rev. D (1985) 1496 [arXiv:hep-lat/0203010].[28] M. Creutz, Algorithms for Simulating Fermions (1992), in Quantum Fieldson the Computer, World Scientific Publishing.[29] A. Alexandru and A. Hasenfratz, Phys. Rev. D (2002) 114506[arXiv:hep-lat/0203026]; (2002) 094502 [arXiv:hep-lat/0207014].[30] F. Knechtli and U. Wolff, (Alpha Coll.) Nucl. Phys. B (2003) 3[arXiv:hep-lat/0303001].[31] M. Hasenbusch, Phys. Lett. B (2001) 177 [arXiv:hep-lat/0107019].[32] P. de Forcrand, Nucl. Phys.proc. Suppl. (1999)[arXiv:hep-lat/9809145].[33] A. D. Kennedy and J. Kuti, Phys. Rev. Lett. (1985) 2473.[34] A. D. Kennedy, J. Kuti, S. Meyer and B. J. Pendleton, Phys. Rev. D (1988) 627.[35] B. Joo, I. Horvath and K, F. Liu, Phys. Rev. D (2003) 074505[arXiv:hep-lat/0112033]. 1836] I. Montvay and E. Scholz, Phys. Lett. B (2005) 73[arXiv:hep-lat/0506006].[37] M. Weingart, ’Stochastic Estimator of the s Quark Determinant in FullQCD Simulation’, Diploma work, Uni. Bern, 2007.[38] R. Sommer, Nucl. Phys. B
839 (1994) [arXiv:hep-lat/9310022].[39] A. Hasenfratz, P. Hasenfratz, D. Hierl, F. Niedermayer and A. Sch¨afer, PoSLAT2006 (2006) 178 [arXiv:hep-lat/0610096].[40] P. Hasenfratz,V. Laliena and F. Niedermayer, Phys. Lett. B (1998)125 [arXiv:hep-lat/9801021].[41] H. Neuberger, Phys. Rev. Lett, (1988) 880; Nucl. Phys. B [FS22](1988) 180.[42] G. Martinelli, C. Pittori, C. T. Sachrajda, M. Testa and A. Vladikas, Nucl.Phys. B (1995) 81 [arXiv:hep-lat/9411010].[43] V. Maillart and F. Niedermayer, arXiv:0807.0030 [hep-lat].[44] E. Franco and V. Lubicz, Nucl. Phys. B (1998) 641[arXiv:hep-ph/9803491].[45] K. G. Chetyrkin and A. Retey, Nucl. Phys. B (2000) 3[arXiv:hep-ph/9910332].[46] K. G. Chetyrkin, Phys. Lett. B (1997) 161 [arXiv:hep-ph/9703278].[47] J. A. M. Vermaseren, S. A. Larin and T. van Ritbergen, Phys. Lett. B (1997) 327 [arXiv:hep-ph/9703284].[48] T. van Ritbergen, J. A. M. Vermaseren and S. A. Larin, Phys. Lett. B (1997) 379 [arXiv:hep-ph/9701390].[49] Y. Aoki et al. , Phys. Rev. D (2008) 054510, arXiv:0712.1061 [hep-lat].[50] W. M. Yao et al. [Particle Data Group], J. Phys. G33