Chiral condensate from the twisted mass Dirac operator spectrum
aa r X i v : . [ h e p - l a t ] O c t Prepared for submission to JHEP
DESY 13-038HU-EP-13/08SFB-CPP-13-17
Chiral condensate from the twisted mass Diracoperator spectrum
Krzysztof Cichy, a,b
Elena Garcia-Ramos, a,c
Karl Jansen a,d a NIC, DESY, Platanenallee 6, 15738 Zeuthen, Germany b Adam Mickiewicz University, Faculty of Physics, Umultowska 85, 61-614 Poznan, Poland c Humboldt Universität zu Berlin, Newtonstr. 15, 12489 Berlin, Germany d Department of Physics, University of Cyprus, P.O. Box 20537, 1678 Nicosia, Cyprus
E-mail: [email protected] , [email protected] , [email protected] Abstract:
We present the results of our computation of the dimensionless chiral conden-sate r Σ / with N f = 2 and N f = 2 + 1 + 1 flavours of maximally twisted mass fermions.The condensate is determined from the Dirac operator spectrum, applying the spectralprojector method proposed by Giusti and Lüscher. We use 3 lattice spacings and severalquark masses at each lattice spacing to perform the chiral and continuum extrapolations.We study the effect of the dynamical strange and charm quarks by comparing our resultsfor N f = 2 and N f = 2 + 1 + 1 dynamical flavours. ontents N f = 2 N f = 2 + 1 + 1 One of the most important phenomena of QCD is the spontaneous breaking of chiral sym-metry. This purely non-perturbative phenomenon was subject to many analyses in LatticeQCD, using different fermion discretizations and methods. In particular, the chiral con-densate – the order parameter of spontaneous chiral symmetry breaking – can be extractedfrom chiral perturbation theory fits of the quark-mass dependence of light pseudoscalar me-son observables [1–9], the topological susceptibility [8, 10–13] or the pion electromagneticform factor [14]. Other methods include using the ǫ -regime expansion and/or chiral ran-dom matrix theory [7, 13, 15–23], Wilson chiral perturbation theory fits of the integratedspectral density [24–26] and a calculation directly from the quark propagator [27, 28]. Asummary of recent results is provided in Ref. [29].The chiral condensate is related to the spectral density of the Dirac operator via theBanks-Casher relation [30]: lim λ → lim m → lim V →∞ ρ ( λ, m ) = Σ π , (1.1)where λ is the modulus of the eigenvalue, m the quark mass, ρ ( λ, m ) the spectral density, V the volume and Σ is the chiral condensate in the infinite-volume and in the chiral limit.Clearly, the triple limit on the left-hand side of the above equation makes it impractical toevaluate the chiral condensate on the lattice.– 1 –owever, recently a method has been proposed [31] to effectively make use of the Banks-Casher relation and explore the chiral properties of QCD on the lattice, in particular tocompute the chiral condensate. The method has also other applications, e.g. it allows tocompute the topological susceptibility or renormalization constants. Briefly, the methodconsists in stochastically evaluating the mode number , i.e. the number of eigenmodes ofthe Dirac operator below some spectral threshold value and using the dependence of thisnumber of eigenmodes on the threshold value to calculate the observable of interest. In thefollowing, we will refer to this method as spectral projectors . One of its essential advantagesfor computing the mode number is the fact that it is very effective in terms of computationalcost – the required computational effort grows linearly with the lattice volume instead ofquadratically, as is the case for a direct computation of eigenmodes and counting theirnumber below the spectral threshold.In this paper, we report our results for the chiral condensate with N f = 2 and N f =2 + 1 + 1 Wilson twisted mass fermions at maximal twist. Preliminary results of ourcomputations for the N f = 2 + 1 + 1 case were presented in Ref. [32]. The paper isorganized as follows. In the second section we provide a short description of the spectralprojector method. In section 3, we describe our lattice setup. Section 4 presents our resultsfor the chiral condensate both in the N f = 2 and the N f = 2 + 1 + 1 case. In section 5,we summarize and compare with other determinations of the chiral condensate found inliterature. An appendix presents our tests of the method. Many interesting properties of the chiral regime of QCD can be understood from the be-haviour of quantities related to the low-lying spectrum of the Dirac operator. One of suchspectral quantities, essential in the determination of the chiral condensate, is the modenumber , i.e. the number of eigenvectors of the massive Hermitian Dirac operator D † D witheigenvalue magnitude smaller than some threshold value M . We will denote this modenumber by ν ( M, µ ) , where µ is the quark mass.Here we provide a short description of the spectral projector method for the computa-tion of ν ( M, µ ) . For a more complete exposition, we refer to the original work of Ref. [31].In this section we assume that we work on an Euclidean lattice, but we will not specify theparticular form of the lattice Dirac operator.If P M is the orthogonal projector to the subspace of fermion fields spanned by the lowestlying eigenmodes of the massive Hermitian Dirac operator D † D with eigenvalues below somethreshold value M , the mode number ν ( M, µ ) can be represented stochastically by: ν ( M, µ ) = h Tr P M i = * N N X j =1 ( η j , P M η j ) + , (2.1)where η , . . . , η N are pseudo-fermion fields added to the theory.The orthogonal projector P M can be approximated by a rational function of D † D : P M ≈ h ( X ) , X = 1 − M ∗ D † D + M ∗ , (2.2)– 2 –here M ∗ is a mass parameter related to the spectral threshold value M . The function: h ( x ) = 12 (cid:0) − xP ( x ) (cid:1) (2.3)is an approximation to the step function θ ( − x ) in the range − ≤ x ≤ , where P ( y ) isin our case the Chebyshev polynomial (of some adjustable degree n ) that minimizes thedeviation: δ = max ǫ ≤ y ≤ | − √ yP ( y ) | (2.4)for some ǫ > . Computing the approximation to the spectral projector P M requires solvingthe following equation an appropriate number of times: ( D † D + M ∗ ) ψ = η (2.5)for a given source field η . Solving this equation is the main computational cost in thecalculation of the mode number. In particular, the computational cost scales linearly V with the volume.One can show [31] that the mode number is a renormalization group invariant, i.e.: ν R ( M R , µ R ) = ν ( M, µ ) , (2.6)where the subscript R denotes renormalized quantities. Note that the spectral thresholdparameter M renormalizes in the same way as the light quark mass (i.e. M R = Z − P M forWilson twisted mass fermions).Finally, we give here the relation between the mode number and the mass-dependentrenormalized chiral condensate: [31] Σ R = π V s − (cid:18) µ R M R (cid:19) ∂∂M R ν R ( M R , µ R ) , (2.7)which is defined to match the chiral condensate to leading order of chiral perturbationtheory. In this section, we will specify the lattice Dirac operator that is used for our work, i.e. theWilson twisted mass Dirac operator. For our computations of the chiral condensate, wehave used gauge field configurations generated by the European Twisted Mass Collaboration(ETMC) with N f = 2 [6, 33, 34] and N f = 2 + 1 + 1 [35–37] dynamical quarks.The gauge action is: S G [ U ] = β X x (cid:16) b X µ,ν =1 Re Tr (cid:0) − P × x ; µ,ν (cid:1) + b X µ = ν Re Tr (cid:0) − P × x ; µ,ν (cid:1)(cid:17) , (3.1) As shown in Ref. [31], the ratio
M/M ∗ depends on the chosen approximation to the projector. For thechoice of Ref. [31], which we apply also in our case, M/M ∗ = 0 . . – 3 –ith β = 6 /g , g the bare coupling and P × , P × are the plaquette and rectangularWilson loops, respectively. For the N f = 2 case, we use the tree-level Symanzik improvedaction [38], i.e. we set b = − , with the normalization condition b = 1 − b . In the caseof N f = 2 + 1 + 1 , we use the Iwasaki action [39, 40], i.e. b = − . .The Wilson twisted mass fermion action for the light, up and down quarks for boththe N f = 2 and N f = 2 + 1 + 1 cases, is given in the so-called twisted basis by: [41–44] S l [ ψ, ¯ ψ, U ] = a X x ¯ χ l ( x ) (cid:0) D W + m ,l + iµ l γ τ (cid:1) χ l ( x ) , (3.2)where m ,l and µ l denote, respectively, the bare untwisted and twisted light quark masses(for shortness, whenever there is no risk of confusion, from now on we will use the symbol µ to denote µ l ). The renormalized light quark mass is given by µ R = Z − P µ . The matrix τ acts in flavour space and χ l = ( χ u , χ d ) is a two-component vector in flavour space, relatedto the one in the physical basis by a chiral rotation. The standard massless Wilson-Diracoperator D W reads: D W = 12 (cid:0) γ µ ( ∇ µ + ∇ ∗ µ ) − a ∇ ∗ µ ∇ µ (cid:1) , (3.3)where ∇ µ and ∇ ∗ µ are the forward and backward covariant derivatives.The twisted mass action for the heavy doublet is given by: [43, 45] S h [ ψ, ¯ ψ, U ] = a X x ¯ χ h ( x ) (cid:0) D W + m ,h + iµ σ γ τ + µ δ τ (cid:1) χ h ( x ) , (3.4)where m ,h denotes the bare untwisted heavy quark mass, µ σ the bare twisted mass withthe twist along the τ direction and µ δ the mass splitting along the τ direction, introducedto make the strange and charm quark masses non-degenerate. The mass parameters µ σ and µ δ are related to the physical renormalized strange m sR and charm m cR quark masses by m s,cR = Z − P ( µ σ ∓ ( Z P /Z S ) µ δ ) . The heavy quark doublet in the twisted basis χ h = ( χ c , χ s ) is again related to the one in the physical basis by a chiral rotation.The twisted mass formulation allows for an automatic O ( a ) improvement of physicalobservables, provided the hopping parameter κ = (8 + 2 am ) − , where m ≡ m ,l = m ,h can be chosen, is tuned to maximal twist by setting it to its critical value, at which thePCAC quark mass vanishes [42, 46–50].The details of the ensembles considered for this work are presented in Tab. 1 for N f = 2 and Tab. 2 for N f = 2 + 1 + 1 . For both cases, they include 3 lattice spacings(from a ≈ . to a ≈ . fm) and up to 5 quark masses at a given lattice spacing.The renormalized light quark masses µ R are in the range from around 15 to 50 MeV. Thevalues of the renormalization constant Z P for different ensembles [53–55], used to convertbare light quark masses µ and bare spectral threshold parameters M to their renormalizedvalues in the MS scheme (at the scale of 2 GeV), are given in Tab. 3, where we also showthe values of r /a (in the chiral limit), used to express our results for the condensate as adimensionless product r Σ / . Our physical lattice extents L for extracting physical results For N f = 2 + 1 + 1 , the mass-independent renormalization constant Z P is extracted as a chiral limit ofa dedicated computation with 4 mass-degenerate flavours – see Refs. [56, 57] for details. – 4 –nsemble β lattice aµ µ R [MeV] κ c L [fm]b . × . × . × . × . × . × . × . × . × . × . × . × Table 1 : Parameters of the N f = 2 gauge ensembles [6, 33, 34]. We show the inverse barecoupling β , lattice size ( L/a ) × ( T /a ) , bare twisted light quark mass aµ , renormalizedquark mass µ R in MeV, critical value of the hopping parameter at which the PCAC massvanishes and physical extent of the lattice L in fm.Ensemble β lattice aµ l µ l,R [MeV] κ c L [fm]A30.32 1.90 × × × × × × × × × × × × × × × Table 2 : Parameters of the N f = 2 + 1 + 1 gauge ensembles [35–37]. We show the inversebare coupling β , lattice size ( L/a ) × ( T /a ) , bare twisted light quark mass µ l , renormalizedquark mass µ l,R in MeV, critical value of the hopping parameter at which the PCAC massvanishes and physical extent of the lattice L in fm.range from 2 to 3 fm (in the temporal direction, we always have T = 2 L ). To check forthe size of finite volume effects, we included different lattice sizes for β = 3 . , aµ = 0 . – 5 – f β a [fm] Z P ( MS , r /a Table 3 : The values of the lattice spacing a [37, 51], r /a [35, 51, 52] and the renormal-ization constant Z P in the MS scheme at the scale of 2 GeV [53–55], for different values of β and N f = 2 and N f = 2 + 1 + 1 .( N f = 2 ) and β = 1 . , aµ = 0 . ( N f = 2 + 1 + 1 ). In this section, we show our results of the calculation of the chiral condensate. First, weillustrate the procedure of extraction of the chiral condensate and discuss the influence ofthe various errors that enter the computation. Then, we analyze finite volume effects inour simulations. Finally, we move on to our chiral and continuum extrapolations.
We show here how to extract the mass-dependent chiral condensate according to Eq. (2.7),illustrating the procedure for ensemble B40.32. Using the spectral projector method, wecomputed the dependence of the mode number on the renormalized spectral thresholdparameter M R for 5 values of M R , from around 2.5 times the renormalized quark mass toaround 120 MeV. Shortly above the latter value one starts to see deviations from the linearregime of ν R ( M R , µ R ) vs. M R (see Appendix A).Fig. 1 shows the dependence of the mode number on the renormalized spectral thresh-old parameter M R . The solid line is a linear fit to all 5 points. The slope of this line ∂ν R ( M R , µ R ) /∂M R determines the value of the mass-dependent chiral condensate accord-ing to Eq. (2.7). The error of this slope includes two sources: the error of the slope ofthe bare mode number as function of the bare threshold parameter M , ∂ν ( M, µ ) /∂M and the error of Z P needed to convert from bare to renormalized quantities. Although ∂ν R ( M R , µ R ) /∂M R appears to be constant as a function of M R within errors, we will takeits value to be the middle point of the chosen fitting interval, see below for details of thefitting intervals considered. Finally, Eq. (2.7) yields, after taking the cubic root : a Σ / = 0 . , where the first error is the one of the slope ∂ν ( M, µ ) /∂M and the other one comes from Z P = 0 . and is dominated by systematic effects – hence, we take it as a systematicerror of our computation. – 6 – ν R M R β =3.90, L/a=32, a µ =0.0040fit Figure 1 : Dependence of the mode number on the renormalized spectral threshold pa-rameter M R for ensemble B40.32. The solid line is a linear fit to all 5 points.fit range lowest aM R highest aM R χ / dof r Σ / Table 4 : Values of r Σ / extracted from different fitting ranges. Every fit includes atleast 3 values of M R . The fit labeled “1 – 5” is the full fit. We estimate the error from thechoice of the fitting range by comparing the value from the full fit with the ones from fits“1 – 4” and “2 – 5”. The error given is statistical only.The value of a Σ / can be further converted to a dimensionless product r Σ / (whichwill be the final result of this paper, after taking the chiral and continuum limits) or to aphysical value in MeV. For the former, we use the value in the chiral limit r /a = 5 . .Since the error of this value is again mostly systematic, we quote it as another systematicerror of r Σ / : r Σ / = 0 . , The values of Σ that we give (also in our plots) are always for the renormalized condensate. – 7 –here the two errors are as above and the third one comes from r /a . For a conversion toMeV, one needs to choose a value of the lattice spacing in physical units. There are severalsuch estimates for ETMC 2-flavour ensembles, giving for β = 3 . values including 0.079 fm[6], 0.085 fm [51] and 0.089 fm [58]. Taking this spread into account, the relative error on thelattice spacing is around 7%, which leads to a similar relative uncertainty in the value of thechiral condensate Σ / in physical units, which amounts to about 20 MeV. This is roughlyan order of magnitude more than other errors entering our computation. Even being lessconservative and using for the error the value quoted in Ref. [51] – 0.085(2) stat (1) syst fm –the error that it yields is still of the order of 10 MeV. Therefore, we decided to give ourfinal results as the dimensionless product r Σ / and we chose not to quote any value inMeV for it until a significantly improved determination of the lattice spacing is available.Another source of the error is the choice of the fitting range and hence the value of M R that enters the square root in Eq. (2.7). Of course, physical results should not dependon this choice, provided that the whole fitting interval lies in the linear regime of themode number vs. M R dependence. Hence, varying the fitting range serves two purposes:establishing whether non-linear effects are already present and checking that the choice of M R in Eq. (2.7) does not influence the final result. The values of r Σ / resulting fromdifferent fitting ranges in M R are shown in Tab. 4. For all fits, χ / d . o . f . is below 1. Thecompatibility of all results and the good values of χ / d . o . f . imply that for this ensemblethe choice of the fitting range and M R does not affect the results in a substantial way . Toquantify this error, we considered the 4-point fits with the lowest or highest value of M R excluded. Excluding the first or last point leads to a similar change of the result and hencewe took the larger of the two as our conservative error from the choice of the fitting range.Note, however, that Tab. 4 implies a systematic tendency towards increasing of Σ whenthe fitting range moves towards higher values of M R . This indicates an onset of non-linearbehaviour for values of M R only slightly above the ones we considered.Finally, our estimate of r Σ / for ensemble B40.32, including all sources of error, is: r Σ / = 0 . stat (38) fit range (39) Z P (53) r /a . We note that the total error is dominated by systematic errors. This means that increasingstatistics would not essentially change our total error. It should be considered an importantadvantage of the method of spectral projectors that rather moderate statistics (in our casearound 230 independent gauge field configurations for this ensemble) leads to a practicallynegligible statistical error. Let us also mention that the quoted statistical error takesautocorrelations fully into account. We performed an analysis of autocorrelations usingtwo methods and found that in general the autocorrelations are small, even at our smallestlattice spacings. For the details of our autocorrelation analysis, we refer to Appendix B.
One of the main sources of systematic effects in Lattice QCD simulations are finite volumeeffects (FVE). In Ref. [31], theoretical arguments were provided that FVE should be small The ensemble B40.32 is somewhat special in this aspect. As we show below, in general, the fitting rangeuncertainty is the most important source of error in our analysis. – 8 – ν / V / - L/a β =3.9, a µ =0.004aM R =0.0225aM R =0.0300aM R =0.0375aM R =0.0450aM R =0.0525 r Σ / ν / V / - L/a β =1.9, a µ =0.004aM R =0.0252aM R =0.0336aM R =0.0420aM R =0.0504 r Σ / Figure 2 : The main plots show the volume dependence of the mode number density ν/V for different values of the renormalized spectral threshold M R . The horizontal bands showthe result at the largest volume. The insets show the volume dependence of the chiralcondensate r Σ / (the error of each point includes the statistical error and the systematicone originating from the choice of the fitting interval). (Left) N f = 2 , β = 3 . , aµ = 0 . ,(right) N f = 2 + 1 + 1 , β = 1 . , aµ = 0 . .for the chiral condensate computed from the mode number – with exponentially smalldifference between the finite volume and infinite volume results of O (exp( − M Λ L/ , where M = 2ΛΣ /F , Λ = p M − µ and F is the pion decay constant in the chiral limit. Sincein practice the mass-dependent chiral condensate is extracted at Λ ≫ µ , the mass M Λ ismuch higher than the pion mass, which typically governs FVE. Hence, one expects that forthe computation of the chiral condensate from the mode number, FVE will be rather small.FVE for the mode number itself were computed in SU(2) chiral perturbation theory [24].The resulting formula leads to a prediction that FVE from lattices with L ≥ fm should besmall, O ( . for M R ≈ O (60 − MeV and renormalized quark masses of O (10 − MeV (with larger FVE at smaller M R ). Indeed, in practice, it was shown in Ref. [31] thatfor L ≥ fm the results deviate from their infinite volume values by less than 1%.To show that it is also the case in our setup, we performed the computation of themode number and the chiral condensate for: • N f = 2 : β = 3 . , aµ = 0 . , lattice extents: L/a =16 , 20, 24 and 32, with corresponding physical values of 1.4, 1.7, 2.0 and 2.7 fm,respectively, • N f = 2 + 1 + 1 : β = 1 . , aµ = 0 . , lattice extents: L/a = 20 , 24 and 32, with corresponding physical values of 1.7, 2.1 and 2.8 fm,respectively.In Fig. 2, we show the volume dependence of the mode number density ν/V for 4-5different values of the renormalized spectral threshold M R . The mode number density canbe computed very precisely. It hence provides a strong test of finite size effects.– 9 –he left plot shows our data for N f = 2 . The results for L/a = 20 and especially
L/a = 16 are systematically lower than
L/a = 32 , signaling large FVE. However, the modenumber density for
L/a = 24 is compatible with the one for
L/a = 32 for 3 intermediatevalues of M R , while it differs by 2-3 σ for the lowest and highest M R , thus changing the slopeof the mode number vs. M R dependence and the extracted chiral condensate (see the insetof Fig. 2 (left)). This change of slope is statistically significant, but it is still a relativelysmall, 1-1.5% effect. Taking into account the uncertainty from the choice of the fittingrange, the final results for the chiral condensate are compatible for all cases – including theones for small volumes, indicating that even if the mode number density goes systematicallydown, the slope of the whole ν ( M R , µ R ) dependence is less affected. We have also tried adescription of FVE in the framework of the formula derived in Ref. [24]. We conclude thatit provides a reasonable agreement with actual lattice data for L/a & , while FVE forsmaller volumes are somewhat underestimated (by a factor of O (2) at L/a = 16 , comparedto the actually observed FVE).In the right plot, we show analogous data for N f = 2 + 1 + 1 . Similarly, we observesignificant finite size effects in the mode number density (and also in the chiral condensate)for L/a = 20 , while
L/a = 24 and
L/a = 32 are always compatible.This allows us to conclude that finite size effects are indeed very small when one reachesa linear lattice extent of around 2 fm (
L/a = 24 at both β = 3 . ( N f = 2 ) and β = 1 . ( N f = 2 + 1 + 1 )). Therefore, we used in our analysis all available ETMC ensembles witha linear lattice extent of at least 2 fm . N f = 2 We now show our results for the 2-flavour case. For each value of β , we have 2-4 sea quarkmasses, according to Tab. 1. For each ensemble, we perform computations of the modenumber at 5 values of the renormalized spectral threshold M R , from around 50 to 120 MeV.We follow the procedure outlined in Sec. 4.1, i.e. we extract the mass-dependent condensatefrom the slope ν R ( M R , µ R ) as function of M R for each ensemble.The results for r Σ / for all considered ensembles are gathered in Tab. 5. These resultsare then used to extrapolate to the chiral limit for each value of β . The chiral correctionsto the mass-dependent condensate were calculated at the next-to-leading order of chiralperturbation theory in Ref. [31]. The obtained formula suggests that the mass-dependentcondensate is equal to the chiral condensate in the chiral limit up to terms linear in µ R and higher order effects. In particular, there are no corrections proportional to µ R ln µ R .Moreover, the size of these chiral corrections is small, as illustrated explicitly in Ref. [31]– inserting the values of low energy constants, it was shown that regardless of the value of M R at which the condensate is extracted, the curvature is very mild. Hence, in practice alinear extrapolation of the mass-dependent condensate to the chiral limit is fully justifiedand we follow this conclusion in our chiral extrapolations. As a check, we tried fits of theNLO formula (inserting values of the low energy constants used in Ref. [31]) and we found The exception to this rule is ensemble d65.32 with L ≈ . fm. We decided to use this ensemble,because without it there would only be one quark mass at β = 4 . and no chiral extrapolation could beperformed. – 10 –nsemble aµ r Σ / b . . . . . . . . . Table 5 : Results for r Σ / for all considered N f = 2 ensembles. The given errors are:statistical, resulting from Z P , resulting from r /a , respectively. We also show results in thechiral limit, where we also give the systematic error from the choice of the fitting interval(4th error). See text for more details.that the differences with respect to the linear extrapolation are negligible compared to ourerrors.Our extrapolations for all three values of β are shown in the main plot of Fig. 3 (weplot r Σ vs. r µ R to allow comparisons between different values of β ). The plotted errorsare only statistical, since in extrapolations at fixed β , the relative errors from Z P and r /a are the same for all quark mass values (we use chirally extrapolated values of Z P and r /a )– we give them in Tab. 5. To estimate the systematic error originating from the choice ofthe fitting range in ν R ( M R , µ R ) vs. M R fits, we repeated all chiral extrapolations for twotailored fitting ranges – excluding the first value of M R (to account for effects of comingtoo close to the renormalized quark mass) or the last value thereof (to account for possibledeviations from the linear behaviour for too high values of M R ).The chiral limit values, with all sources of error, are also shown in Tab. 5. In general,the total error originates in practice only from the choice of the fitting range and the latterincreases when approaching the continuum limit. The reasons for this behaviour include thefact that the number of quark masses that we use decreases for smaller lattice spacings andat β = 4 . the slope of the quark mass dependence of the chiral condensate is apparentlylarger than at coarser lattice spacings , making the final chiral limit value more susceptibleto changes in the fitting interval.Finally, we can use the chirally extrapolated values of the condensate to perform anextrapolation to the continuum limit. We start by discussing the O ( a ) -improvement of thechiral condensate. For on-shell quantities, O ( a ) -improvement amounts to the quantity being Note that this slope may be affected by the smaller volume of ensemble d65.32 – hence, it may be aresidual FVE and not an indication of the dependence of the slope on the lattice spacing. In such case, ourfitting range error at β = 4 . implicitly reflects this FVE. – 11 – r Σ r µ R β =3.90 β =4.05 β =4.20 r Σ / (a/r ) N f =2 cont.limit = 0.689(16) Figure 3 : Main plot: chiral extrapolations of the chiral condensate r Σ for N f = 2 ensembles and β = 3 . , 4.05, 4.2. The lines are extrapolations to the chiral limit, linear inthe quark mass. The values in the chiral limit for β = 3 . and 4.2 are slightly shifted forbetter presentation. The errors are statistical only. Inset: continuum extrapolation of thechirally extrapolated chiral condensate r Σ / vs. ( a/r ) . The errors include: statisticalerrors, errors from Z P and errors from r /a .even under the R parity transformation: ψ → iγ τ ψ , ¯ ψ → i ¯ ψγ τ [42]. Let us considerthe spectral sums [31]: σ k ( µ, m q ) = h Tr n ( D † m D m + µ ) − k o i , where k ≥ for reasonsexplained in the given reference. The spectral sums are related to the mode number [31]and the improvement (or lack thereof) of the spectral sums implies the improvement of thechiral condensate. Representing the spectral sums as a density chain correlation function(for k = 3 ): σ ( µ, m ) = − a X x ...x h P +12 ( x ) P − ( x ) P +34 ( x ) P − ( x ) P +56 ( x ) P − ( x ) i , (4.1)it is straightforward to show that the object on the right-hand side is even under R transformation, since the number of densities is even. However, one also needs to considercontact terms arising from Eq. (4.1), i.e. terms in the sum with x i = x j for some i = j .It can be demonstrated [59] that such terms give rise only to O ( am ) terms in the modenumber – hence they vanish at maximal twist. In this way, the contact terms do not spoilautomatic O ( a ) -improvement of the chiral condensate.– 12 –ence, our continuum limit extrapolation is performed linearly in a , using results atthree lattice spacings, with fixed fitting range of the mode number vs. spectral thresholddependence, corresponding to M R ≈ MeV (entering the square root in Eq. (2.7)) forall values of β . As an error, we use the statistical error, combined in quadrature withthe error of Z P and r /a . We do not observe significant cut-off effects. The final valuein the continuum limit is 0.689(16). To account for the fitting range error, we performthe full analysis for tailored fitting ranges, excluding the first or last value of M R for eachensemble. This corresponds to a shift in M R to approx. 80 or 100 MeV, respectively. Whilethe extracted value of the condensate in the chiral limit should not depend on the fittingrange, in practice the results for different fitting ranges differ, which is due to using only4-5 values in the fits to extract the slope of ν R ( M R , µ R ) . The fits from tailored fittingranges yield 0.678(18) and 0.718(20), respectively. To be conservative, as our systematicerror from the fitting range we choose the larger difference of the two with respect to thecentral value 0.689(16). This finally gives: r Σ / N f =2 = 0 . , where the first error is the combined statistical error, the error of Z P and of r /a , while thesecond error originates from the choice of the fitting range. N f = 2 + 1 + 1 In this subsection, we present results for the N f = 2 + 1 + 1 case. By comparing to theresults of the 2-flavour case, we can investigate the role of the dynamical strange and charmquarks.We proceed as in the previous section. For each value of β , we have 3-5 sea quark masses,according to Tab. 2. We compute the mode number at 4 values of the renormalized spectralthreshold M R , from around 50 to 110 MeV, and extract the mass-dependent condensatefrom the slope of the ν R ( M R , µ R ) vs. M R dependence for each ensemble. We have alsocomputed the mode number for a fifth value of M R ≈ MeV. However, given the verygood statistical precision of the spectral projectors method of evaluating the mode number,we observe significant deviations from the linear dependence of ν R ( M R , µ R ) on M R whenthis fifth value of M R is included. Because of this, we decided not to include the results at M R ≈ MeV.In Tab. 6, we show all our results for r Σ / in the 2+1+1-flavour case. We also includethe results of a linear extrapolation to the chiral limit for each value of β , shown in themain plot of Fig. 4. As before, we plot only statistical errors, since all extrapolations areperformed at fixed β and the errors from Z P and r /a are the same for all quark masses(we use chirally extrapolated values of Z P and r /a ) – given in Tab. 6. Contrary to the N f = 2 case, the slope of the dependence of the mass-dependent condensate on the lightquark mass slightly decreases for increasing β (the change of slope is statistically significantwhen going from β = 1 . to β = 2 . ). This has the effect of lowering the systematic errorrelated to the choice of the fitting range for decreasing lattice spacing , which is for all β Moreover, we always have at least 3 sea quark masses for each β in the N f = 2 + 1 + 1 case, comparedto only 2 masses at β = 4 . ( N f = 2 )). – 13 –nsemble aµ r Σ / A . . . . . . . . . . . . Table 6 : Results for r Σ / for all considered N f = 2 + 1 + 1 ensembles. The given errorsare: statistical, resulting from Z P , resulting from r /a , respectively. We also show resultsin the chiral limit, where we also give the systematic error from the choice of the fittinginterval (4th error). See text for more details.the dominating source of error (although for β = 2 . other errors become comparable).The chirally extrapolated values at three lattice spacings are then used to perform anextrapolation to the continuum limit, which is again compatible with O ( a ) cut-off effects.To estimate the fitting range uncertainty, we again perform 3 separate continuum limitextrapolations, using different fitting ranges and different values of M R , corresponding toapprox. 80, 90 and 100 MeV. The values of r Σ / in the continuum limit are, respectively,0.668(24), 0.680(20) and 0.659(27). As our central value we take the result from the fullfitting range: r Σ / N f =2+1+1 = 0 . , with the larger of the differences with respect to values from tailored fitting intervals as thefitting range systematic error. This can be compared to the 2-flavour result which amountsto r Σ / N f =2 = 0 . and both results are compatible within errors. In this paper, we presented our results on the chiral condensate in QCD with N f = 2 and N f = 2 + 1 + 1 flavours of dynamical Wilson twisted mass quarks at maximal twist.Our final results are: r Σ / N f =2 = 0 . ,r Σ / N f =2+1+1 = 0 . , – 14 – r Σ r µ R β =1.90 β =1.95 β =2.10 r Σ / (a/r ) N f =2+1+1 cont.limit = 0.680(20) Figure 4 : Main plot: chiral extrapolations of the chiral condensate r Σ for N f = 2 + 1 + 1 ensembles and β = 1 . , 1.95, 2.1. The lines are extrapolations to the chiral limit, linearin the quark mass. The errors are statistical only. Inset: continuum extrapolation of thechirally extrapolated chiral condensate r Σ / vs. ( a/r ) . The errors include: statisticalerrors, errors from Z P and errors from r /a .which indicates that at the current level of precision, we cannot discriminate the influence ofthe dynamical strange and charm quarks on the value of the light quark chiral condensate.The main source of the error of our results is the systematic error related to the choiceof the fitting range in the dependence of the renormalized mode number on the renormalizedspectral threshold. The second most important source of the error is either the statisticalerror or the error related to the uncertainty in the values of r /a (which are inputs of ouranalysis). The error from Z P (also an input of our analysis), used to renormalize the quarkmasses and the spectral threshold parameter, is usually the smallest. However, we want toemphasize that in all cases the fitting range error is the largest one and in most cases it islarger by a factor of 2-4 than any other error. The rather small statistical errors that weobtain indicate that increasing statistics would not make our total error significantly smaller.This implies that a way to improve the total error would be to increase the number of valuesof the spectral threshold M R at which one computes the mode number. This would allowto identify more precisely the linear region of the mode number vs. M R dependence (seediscussion in the Appendix) – on the one hand sufficiently far away from the renormalized– 15 –esult method N f fermions r Σ / this work spectral proj. 2 twisted mass 0.689(16)(29)this work spectral proj. 2+1+1 twisted mass 0.680(20)(21)RBC-UKQCD [3] chiral fits 2+1 domain wall 0.632(15)(12)MILC [4] chiral fits 2+1 staggered 0.654(14)(18)MILC [5] chiral fits 2+1 staggered 0.653(18)(11)S. Borsanyi et al. [9] chiral fits 2+1 staggered 0.662(5)(20)ETMC [6] chiral fits 2 twisted mass 0.575(14)(52)ETMC [27] quark propagator 2 twisted mass 0.676(89)(14)HPQCD [28] quark propagator 2+1+1 staggered 0.673(5)(11) Table 7 : Comparison of our results for r Σ / with large-volume continuum limit resultsfound in literature. In the given references, the values of Σ are given in MeV. To convertto the dimensionless product r Σ / , we combine them with results for r . The first errorgiven is always from the computation of the value in MeV (when there are several errorsgiven, we combine them in quadrature) and the second one from conversion using physicalvalue of r . For RBC-UKQCD, Σ / = 256(6) MeV, r = 0 . fm [3]. For MILC [4], Σ / = 278(6) MeV, r = 0 . fm [4], r /r =1.46(1)(2) [60]. For MILC [5], Σ / =281 . . MeV, r = 0 . fm [5], r /r =1.46(1)(2) [60]. For S. Borsanyi et al. [9], Σ / = 272 . . . MeV, r = 0 . fm [61]. For ETMC [6], Σ / = 269 . . MeV, r = 0 . fm. However, newer analyses indicate a higher value of r ≈ . fm[51]. To take this into account, we added the spread of the new and old value as a systematicerror and used r = 0 . fm to calculate r Σ / from Σ / in MeV. For ETMC [27], Σ / = 299(26)(29) MeV, r = 0 . fm [6]. For HPQCD [28], Σ / = 283(2) MeV, r = 0 . fm [62], r /r =1.46(1)(2) [60].quark mass, on the other hand low enough such that there are no deviations from the linearbehaviour at the upper end of the fitting range (observed already at M R ≈ MeV in the N f = 2 + 1 + 1 case).In order to place our values of the chiral condensate in context of results of other col-laborations, we attempt in Tab. 7 a comparison. Given the large amount of approachesto compute the chiral condensate, as mentioned in the introduction, we make a selectionby only considering results that are given in the literature as continuum limit values fromlarge volume simulations. Note that the other available results that are listed in Tab. 7are obtained in a different way than the spectral projector method. They mostly use chiralperturbation theory fits to the quark mass dependence of light pseudoscalar meson observ-ables or determine the chiral condensate from the quark propagator. Not discussing theadvantages or disadvantages of the various fermion discretizations used, we see in Tab. 7 anoverall agreement for the dimensionless quantity r Σ / , which is reassuring and establishes,in our opinion, the spectral projector method as a valuable alternative to determine thechiral condensate. On the other hand, when looking at the chiral condensate in physicalunits, see the caption of Tab. 7, a spread of results is obtained. Thus, it seems that the scale– 16 –etting from the different lattice calculations introduces a systematic effect and it would bedesirable to clarify this uncertainty in future more precise calculations. Acknowledgments
We thank the European Twisted Mass Collaboration for generatinggauge field configurations ensembles used for this work. We are grateful to A. Shindlerfor discussions and in particular suggestions concerning O ( a ) improvement of the modenumber. We thank G. Herdoiza for discussions at various stages of this work and for hiscareful reading of the manuscript. We acknowledge useful discussions with P. Damgaard,V. Drach, M. Lüscher, K. Ottnad, G.C. Rossi, S. Sharpe, C. Urbach, F. Zimmermann.K.C. was supported by Foundation for Polish Science fellowship “Kolumb”. This workwas supported in part by the DFG Sonderforschungsbereich/Transregio SFB/TR9. K.J.was supported in part by the Cyprus Research Promotion Foundation under contract Π PO Σ E Λ KY Σ H/EM Π EIPO Σ /0311/16. The computer time for this project was madeavailable to us by the Jülich Supercomputing Center, LRZ in Munich, the PC cluster inZeuthen, Poznan Supercomputing and Networking Center (PCSS). We thank these com-puter centers and their staff for all technical advice and help. A Testing the method
To test our implementation of the method in the tmLQCD package [63], we compared thespectral projectors results for the mode number with ones from explicit computation of150 lowest eigenvalues of D † D for each gauge field configuration, using ensemble B85.24.The results of this comparison are shown in the right plot of Fig. 5, where the 4 pointscorrespond to the stochastically evaluated mode number using spectral projectors, whilethe continuous line (which has an error roughly of the order of the width of the line) is theresult from explicitly computing the eigenvalues. We observe very good agreement betweenthe two methods.Moreover, we used the results of explicit computation of eigenvalues to estimate theregion of renormalized spectral threshold of M R where we observe linear dependence be-tween the renormalized mode number and the threshold value (left plot of Fig. 5). Theonset of non-linear behaviour corresponds to approx. 130-150 MeV. On the other end of thespectrum, one clearly observes effects of M R close to the renormalized quark mass up toaround 10-20 MeV above the latter . This allows us to identify the linear region to extendbetween around 60 and 120 MeV, to allow for some safety margin. Therefore, we decidedto choose our values of M R for the computation of the chiral condensate roughly in thisinterval. We remark that values in this range were used in Ref. [31].Some parameters employed in the method of spectral projectors need to be tuned toobtain a compromise between the accuracy of results and computational cost.First of all, as mentioned in Ref. [31], the precision of the inverter can be chosen to berelatively sloppy without reducing the accuracy. In order to identify the optimal precision,which does not affect the correctness of the result, but still decreases the computational One expects that the behaviour close to the renormalized quark mass is modified in an important wayby lattice artefacts [24]. – 17 – pect. proj β = 1 . aµ = 0 . M [MeV] m o d e nu m b e r ν β = 1 . aµ = 0 . M [MeV] m o d e nu m b e r ν Figure 5 : The mode number as a function of M R for the ensemble B85.24. The solid linein both plots is the result of explicit computation of 150 eigenvalues for each gauge fieldconfiguration, while data points represent results of spectral projectors calculation of themode number for 4 values of M R . The right plot is the zoom of the left part. CG prec = 10 − CG prec = 10 − CG prec = 10 − CG prec = 10 − CG prec = 10 − CG prec = 10 − N f = 2 , β = 3 . , × m o d e nu m b e r ν Figure 6 : The mode number vs. the number of averaged stochastic sources using differentvalues of the inverter relative precision for a × lattice at β = 3 . , aµ = 0 . (ensemble b40.16).time, we computed the mode number for several values of the relative precision of theinverter. The results (for ensemble b40.16) are shown in Fig. 6, which shows that evenprecision of − gives reasonable result. However, we decided to be conservative and wechose a value of − for the relative precision of the inverter.We also checked the dependence of the mode number on the number of stochasticsources used for each gauge field configuration, shown again in Fig. 6. We observe thatall results are compatible within error, which matches the suggestion of Ref. [31] that onestochastic source should be enough. However, we observe that adding a second source– 18 –nsemble number step τ int (boot) τ int (UW) τ int (boot) τ int (UW)of confs HMC traj. smallest M largest M A30.32 116 20 0.7 0.6(2) 1.7 1.9(8)A40.20 197 16 0.8 0.9(3) 0.8 0.9(3)A40.24 185 20 0.6 0.7(2) 1.5 2.2(8)A40.32 201 8 0.8 0.7(2) 0.8 0.7(2)A50.32 201 20 0.5 0.7(2) 0.5 0.8(2)A60.24 161 8 0.7 0.5(1) 1.3 1.6(6)A80.24 200 8 0.6 0.7(2) 0.7 0.9(2)B25.32 200 20 1.0 0.7(2) 1.0 1.2(4)B35.32 199 20 0.4 0.5(1) 1.0 1.7(6)B55.32 201 20 0.4 0.4(1) 0.6 0.7(2)B75.32 199 8 0.7 0.6(1) 0.5 0.5(1)B85.24 196 20 0.4 0.4(1) 0.5 0.5(1)D15.48 119 20 0.5 0.5(1) 0.6 0.8(2)D20.48 193 20 0.6 0.5(1) 0.7 0.7(2)D30.48 203 20 0.4 0.4(1) 0.5 0.4(1)
Table 8 : Autocorrelations in the mode number, N f = 2 + 1 + 1 ensembles. We give thenumber of gauge field configurations used for each ensemble, the step in units of HMCtrajectories and the calculated values of τ int using two methods: bootstrap with blocking(boot) and the method proposed by U. Wolff [64] (UW).might help to considerably reduce the statistical error, which may be important for shorterMonte Carlo runs, when the available number of independent gauge field configurations israther small. Adding further sources does not change the error considerably, because ofcorrelations between results obtained from the same gauge field configuration. B Autocorrelations
In this appendix, we show the results of our autocorrelation analysis of the mode number(for the smallest and largest value of M that we used for chiral condensate extraction).We applied two methods: bootstrap with blocking (with block size of 10 measurements)and the method proposed by U. Wolff in Ref. [64]. Our results are shown in Tabs. 8( N f = 2 + 1 + 1 ) and 9 ( N f = 2 ). In general, both methods yield compatible results for theintegrated autocorrelation time τ int (in the case of the method of Ref. [64], we quote alsothe error of τ int and the two values of τ int that we obtain agree within this error).Our conclusions about the dependence of τ int on simulation parameters are the follow-ing: • at the smallest value of M , no autocorrelations are observed ( τ int compatible with0.5), – 19 –nsemble number step τ int (boot) τ int (UW) τ int (boot) τ int (UW)of confs HMC traj. smallest M largest M b .
134 8 0.3 0.4(1) 0.9 0.9(3)b .
544 10 0.5 0.5(1) 1.6 1.8(4)b .
265 20 0.6 0.5(1) 1.2 1.3(3)b .
465 20 0.5 0.7(2) 0.8 1.3(3)b .
232 16 0.4 0.4(1) 0.5 0.5(1)b .
272 20 0.4 0.4(1) 0.6 0.8(2)b .
187 20 0.3 0.4(1) 0.5 0.5(1)c .
183 8 0.7 0.5(1) 1.2 1.6(5)c .
123 20 0.5 0.5(1) 0.4 0.5(1)c .
201 10 0.4 0.5(1) 1.1 0.8(2)d .
77 10 0.7 0.5(1) 0.5 0.6(2)d .
199 20 0.6 0.5(1) 0.7 1.0(3)
Table 9 : Autocorrelations in the mode number, N f = 2 ensembles. We give the numberof gauge field configurations used for each ensemble, the step in units of HMC trajectoriesand the calculated values of τ int using two methods: bootstrap with blocking (boot) andthe method proposed by U. Wolff [64] (UW). • at the largest value of M , in some cases the autocorrelations become visible, with τ int between 1 and 2, • there is a tendency towards larger autocorrelations for smaller quark masses (e.g. B25compared to B85) and for smaller volumes (e.g. b40 ensembles at 4 volumes), • we don’t observe a tendency towards increased τ int for decreasing lattice spacing. References [1]
PACS-CS
Collaboration, S. Aoki et al., , Phys.Rev.
D79 (2009) 034503, [ arXiv:0807.1661 ].[2]
RBC-UKQCD
Collaboration, C. Allton et al.,
Physical Results from 2+1 Flavor DomainWall QCD and SU(2) Chiral Perturbation Theory , Phys.Rev.
D78 (2008) 114509,[ arXiv:0804.0473 ].[3]
RBC-UKQCD
Collaboration, Y. Aoki et al.,
Continuum Limit Physics from 2+1 FlavorDomain Wall QCD , Phys.Rev.
D83 (2011) 074508, [ arXiv:1011.0892 ].[4] A. Bazavov, D. Toussaint, C. Bernard, J. Laiho, C. DeTar, et al.,
Nonperturbative QCDsimulations with 2+1 flavors of improved staggered quarks , Rev.Mod.Phys. (2010)1349–1417, [ arXiv:0903.3598 ].[5] A. Bazavov, C. Bernard, C. DeTar, X. Du, W. Freeman, et al., Staggered chiral perturbationtheory in the two-flavor case and SU(2) analysis of the MILC data , PoS
LATTICE2010 (2010) 083, [ arXiv:1011.1792 ]. – 20 – ETMC
Collaboration, R. Baron et al.,
Light Meson Physics from Maximally Twisted MassLattice QCD , JHEP (2010) 097, [ arXiv:0911.5061 ].[7]
JLQCD, TWQCD
Collaboration, H. Fukaya et al.,
Determination of the chiral condensatefrom QCD Dirac spectrum on the lattice , Phys.Rev.
D83 (2011) 074501, [ arXiv:1012.4052 ].[8] V. Bernard, S. Descotes-Genon, and G. Toucas,
Topological susceptibility on the lattice andthe three-flavour quark condensate , JHEP (2012) 051, [ arXiv:1203.0508 ].[9] S. Borsanyi, S. Durr, Z. Fodor, S. Krieg, A. Schafer, et al.,
SU(2) chiral perturbation theorylow-energy constants from 2+1 flavor staggered lattice simulations , arXiv:1205.0788 .[10] S. Durr, Topological susceptibility in full QCD: Lattice results versus the prediction from theQCD partition function with granularity , Nucl.Phys.
B611 (2001) 281–310,[ hep-lat/0103011 ].[11]
TWQCD
Collaboration, T.-W. Chiu, T.-H. Hsieh, and P.-K. Tseng,
Topologicalsusceptibility in 2+1 flavors lattice QCD with domain-wall fermions , Phys.Lett.
B671 (2009)135–138, [ arXiv:0810.3406 ].[12] T. W. Chiu, T. H. Hsieh, and Y. Y. Mao,
Topological Susceptibility in Two Flavors LatticeQCD with the Optimal Domain-Wall Fermion , Phys.Lett.
B702 (2011) 131–134,[ arXiv:1105.4414 ].[13] F. Bernardoni, P. Hernandez, N. Garron, S. Necco, and C. Pena,
Probing the chiral regime of N f = 2 QCD with mixed actions , Phys.Rev.
D83 (2011) 054503, [ arXiv:1008.1870 ].[14]
ETMC
Collaboration, R. Frezzotti, V. Lubicz, and S. Simula,
Electromagnetic form factorof the pion from twisted-mass lattice QCD at N(f ) = 2 , Phys.Rev.
D79 (2009) 074506,[ arXiv:0812.4042 ].[15] P. Hernandez, K. Jansen, and L. Lellouch,
Finite size scaling of the quark condensate inquenched lattice QCD , Phys.Lett.
B469 (1999) 198–204, [ hep-lat/9907022 ].[16] P. Damgaard, U. M. Heller, R. Niclasen, and K. Rummukainen,
Eigenvalue distributions ofthe QCD Dirac operator , Phys.Lett.
B495 (2000) 263–270, [ hep-lat/0007041 ].[17] T. A. DeGrand and S. Schaefer,
Chiral properties of two-flavor QCD in small volume and atlarge lattice spacing , Phys.Rev.
D72 (2005) 054503, [ hep-lat/0506021 ].[18] C. Lang, P. Majumdar, and W. Ortner,
The Condensate for two dynamical chirally improvedquarks in QCD , Phys.Lett.
B649 (2007) 225–229, [ hep-lat/0611010 ].[19]
JLQCD
Collaboration, H. Fukaya et al.,
Two-flavor lattice QCD simulation in theepsilon-regime with exact chiral symmetry , Phys.Rev.Lett. (2007) 172001,[ hep-lat/0702003 ].[20] TWQCD
Collaboration, H. Fukaya et al.,
Two-flavor lattice QCD in the epsilon-regime andchiral Random Matrix Theory , Phys.Rev.
D76 (2007) 054503, [ arXiv:0705.3322 ].[21] A. Hasenfratz, R. Hoffmann, and S. Schaefer,
Low energy chiral constants fromepsilon-regime simulations with improved Wilson fermions , Phys.Rev.
D78 (2008) 054511,[ arXiv:0806.4586 ].[22] K. Jansen and A. Shindler,
The Epsilon regime of chiral perturbation theory with Wilson-typefermions , PoS
LAT2009 (2009) 070, [ arXiv:0911.1931 ].[23] K. Splittorff and J. Verbaarschot,
The Microscopic Twisted Mass Dirac Spectrum , Phys.Rev.
D85 (2012) 105008, [ arXiv:1201.1361 ]. – 21 –
24] S. Necco and A. Shindler,
Spectral density of the Hermitean Wilson Dirac operator: a NLOcomputation in chiral perturbation theory , JHEP (2011) 031, [ arXiv:1101.1778 ].[25] S. Necco and A. Shindler,
On the spectral density of the Wilson operator , PoS
LATTICE2011 (2011) 250, [ arXiv:1108.1950 ].[26] S. Necco and A. Shindler,
Corrections to the Banks-Casher relation with Wilson quarks , arXiv:1302.5595 .[27] F. Burger, V. Lubicz, M. Muller-Preussker, S. Simula, and C. Urbach, Quark mass and chiralcondensate from the Wilson twisted mass lattice quark propagator , arXiv:1210.0838 .[28] C. McNeile, A. Bazavov, C. Davies, R. Dowdall, K. Hornbostel, et al., Direct determinationof the strange and light quark condensates from full lattice QCD , arXiv:1211.6577 .[29] G. Colangelo, S. Durr, A. Juttner, L. Lellouch, H. Leutwyler, et al., Review of lattice resultsconcerning low energy particle physics , Eur.Phys.J.
C71 (2011) 1695, [ arXiv:1011.4408 ].[30] T. Banks and A. Casher,
Chiral Symmetry Breaking in Confining Theories , Nucl.Phys.
B169 (1980) 103.[31] L. Giusti and M. Luscher,
Chiral symmetry breaking and the Banks-Casher relation in latticeQCD with Wilson quarks , JHEP (2009) 013, [ arXiv:0812.3638 ].[32] K. Cichy, V. Drach, E. Garcia-Ramos, and K. Jansen,
Topological susceptibility and chiralcondensate with N f = 2 + 1 + 1 dynamical flavors of maximally twisted mass fermions , PoS
LATTICE2011 (2011) 102, [ arXiv:1111.3322 ].[33]
ETMC
Collaboration, P. Boucaud et al.,
Dynamical twisted mass fermions with lightquarks , Phys.Lett.
B650 (2007) 304–311, [ hep-lat/0701012 ].[34]
ETMC
Collaboration, P. Boucaud et al.,
Dynamical Twisted Mass Fermions with LightQuarks: Simulation and Analysis Details , Comput.Phys.Commun. (2008) 695–715,[ arXiv:0803.0224 ].[35]
ETMC
Collaboration, R. Baron, P. Boucaud, J. Carbonell, A. Deuzeman, V. Drach, et al.,
Light hadrons from lattice QCD with light (u,d), strange and charm dynamical quarks , JHEP (2010) 111, [ arXiv:1004.5284 ].[36]
ETMC
Collaboration, R. Baron et al.,
Computing K and D meson masses with N f =2+1+1 twisted mass lattice QCD , Comput.Phys.Commun. (2011) 299–316,[ arXiv:1005.2042 ].[37]
ETMC
Collaboration, R. Baron et al.,
Light hadrons from Nf=2+1+1 dynamical twistedmass fermions , PoS
LATTICE2010 (2010) 123, [ arXiv:1101.0518 ].[38] P. Weisz,
Continuum Limit Improved Lattice Action for Pure Yang-Mills Theory. 1. , Nucl.Phys.
B212 (1983) 1.[39] Y. Iwasaki,
Renormalization Group Analysis of Lattice Theories and Improved LatticeAction: Two-Dimensional Nonlinear O(N) Sigma Model , Nucl.Phys.
B258 (1985) 141–156.[40] Y. Iwasaki, K. Kanaya, T. Kaneko, and T. Yoshie,
Scaling in SU(3) pure gauge theory with arenormalization group improved action , Phys.Rev.
D56 (1997) 151–160, [ hep-lat/9610023 ].[41]
Alpha
Collaboration, R. Frezzotti, P. A. Grassi, S. Sint, and P. Weisz,
Lattice QCD with achirally twisted mass term , JHEP (2001) 058, [ hep-lat/0101001 ].[42] R. Frezzotti and G. Rossi,
Chirally improving Wilson fermions. 1. O(a) improvement , JHEP (2004) 007, [ hep-lat/0306014 ]. – 22 –
43] R. Frezzotti and G. Rossi,
Chirally improving Wilson fermions. II. Four-quark operators , JHEP (2004) 070, [ hep-lat/0407002 ].[44] A. Shindler,
Twisted mass lattice QCD , Phys.Rept. (2008) 37–110, [ arXiv:0707.4093 ].[45] R. Frezzotti and G. Rossi,
Twisted mass lattice QCD with mass nondegenerate quarks , Nucl.Phys.Proc.Suppl. (2004) 193–202, [ hep-lat/0311008 ].[46] T. Chiarappa, F. Farchioni, K. Jansen, I. Montvay, E. Scholz, et al.,
Numerical simulation ofQCD with u, d, s and c quarks in the twisted-mass Wilson formulation , Eur.Phys.J.
C50 (2007) 373–383, [ hep-lat/0606011 ].[47] F. Farchioni, C. Urbach, R. Frezzotti, K. Jansen, I. Montvay, et al.,
Exploring the phasestructure of lattice QCD with twisted mass quarks , Nucl.Phys.Proc.Suppl. (2005)240–245, [ hep-lat/0409098 ].[48] F. Farchioni, K. Jansen, I. Montvay, E. Scholz, L. Scorzato, et al.,
The Phase structure oflattice QCD with Wilson quarks and renormalization group improved gluons , Eur.Phys.J.
C42 (2005) 73–87, [ hep-lat/0410031 ].[49] R. Frezzotti, G. Martinelli, M. Papinutto, and G. Rossi,
Reducing cutoff effects in maximallytwisted lattice QCD close to the chiral limit , JHEP (2006) 038, [ hep-lat/0503034 ].[50]
XLF
Collaboration, K. Jansen, M. Papinutto, A. Shindler, C. Urbach, and I. Wetzorke,
Quenched scaling of Wilson twisted mass fermions , JHEP (2005) 071,[ hep-lat/0507010 ].[51]
ETMC
Collaboration, B. Blossier et al.,
Average up/down, strange and charm quark masseswith Nf=2 twisted mass lattice QCD , Phys.Rev.
D82 (2010) 114513, [ arXiv:1010.3659 ].[52]
ETMC
Collaboration, K. Ottnad et al., η and η ′ mesons from Nf=2+1+1 twisted masslattice QCD , JHEP (2012) 048, [ arXiv:1206.6719 ].[53]
ETMC
Collaboration, M. Constantinou et al.,
Non-perturbative renormalization of quarkbilinear operators with Nf = 2 (tmQCD) Wilson fermions and the tree-level improved gaugeaction , JHEP (2010) 068, [ arXiv:1004.1115 ].[54] C. Alexandrou, M. Constantinou, T. Korzec, H. Panagopoulos, and F. Stylianou,
Renormalization constants of local operators for Wilson type improved fermions , Phys.Rev.
D86 (2012) 014505, [ arXiv:1201.5025 ].[55] ETMC preliminary result for N f = 4 renormalization constants, private communication fromD. Palao.[56] ETMC
Collaboration, P. Dimopoulos et al.,
Renormalization constants for Wilson fermionlattice QCD with four dynamical flavours , PoS
LATTICE2010 (2010) 235,[ arXiv:1101.1877 ].[57]
ETMC
Collaboration, B. Blossier et al.,
Renormalisation constants of quark bilinears inlattice QCD with four dynamical Wilson quarks , PoS
LATTICE2011 (2011) 233,[ arXiv:1112.1540 ].[58]
ETMC
Collaboration, C. Alexandrou et al.,
Light baryon masses with dynamical twistedmass fermions , Phys.Rev.
D78 (2008) 014509, [ arXiv:0803.3190 ].[59] K. Cichy, E. Garcia-Ramos, K. Jansen, and A. Shindler in preparation.[60] C. Bernard, C. E. DeTar, L. Levkova, S. Gottlieb, U. Heller, et al.,
Status of the MILC lightpseudoscalar meson project , PoS
LAT2007 (2007) 090, [ arXiv:0710.1118 ]. – 23 –
61] Y. Aoki, S. Borsanyi, S. Durr, Z. Fodor, S. D. Katz, et al.,
The QCD transition temperature:results with physical masses in the continuum limit II. , JHEP (2009) 088,[ arXiv:0903.4155 ].[62]
HPQCD
Collaboration, R. Dowdall et al.,
The Upsilon spectrum and the determination ofthe lattice spacing from lattice QCD including charm quarks in the sea , Phys.Rev.
D85 (2012) 054509, [ arXiv:1110.6887 ].[63] K. Jansen and C. Urbach, tmLQCD: A Program suite to simulate Wilson Twisted massLattice QCD , Comput.Phys.Commun. (2009) 2717–2738, [ arXiv:0905.3331 ].[64]
ALPHA
Collaboration, U. Wolff,
Monte Carlo errors with less errors , Comput.Phys.Commun. (2004) 143–153, [ hep-lat/0306017 ].].