Chiral condensate from the twisted mass Dirac operator spectrum

We present the results of our computation of the chiral condensate 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 spectral projector method proposed by Giusti and Luscher. We use 3 lattice spacings and several quark 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 results for $N_f=2$ and $N_f=2+1+1$ dynamical flavours.

The chiral condensate is related to the spectral density of the Dirac operator via the Banks-Casher relation [30]: 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 to evaluate the chiral condensate on the lattice.
However, 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 to compute the chiral condensate. The method has also other applications, e.g. it allows to compute the topological susceptibility or renormalization constants. Briefly, the method consists in stochastically evaluating the mode number, i.e. the number of eigenmodes of the Dirac operator below some spectral threshold value and using the dependence of this number of eigenmodes on the threshold value to calculate the observable of interest. In the following, we will refer to this method as spectral projectors. One of its essential advantages for computing the mode number is the fact that it is very effective in terms of computational cost -the required computational effort grows linearly with the lattice volume instead of quadratically, as is the case for a direct computation of eigenmodes and counting their number 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 our computations for the N f = 2 + 1 + 1 case were presented in Ref. [32]. The paper is organized as follows. In the second section we provide a short description of the spectral projector method. In section 3, we describe our lattice setup. Section 4 presents our results for 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 in literature. An appendix presents our tests of the method.

Spectral projectors and chiral condensate
Many interesting properties of the chiral regime of QCD can be understood from the behaviour of quantities related to the low-lying spectrum of the Dirac operator. One of such spectral quantities, essential in the determination of the chiral condensate, is the mode number, i.e. the number of eigenvectors of the massive Hermitian Dirac operator D † D with eigenvalue magnitude smaller than some threshold value M 2 . We will denote this mode number by ν(M, µ), where µ is the quark mass.
Here we provide a short description of the spectral projector method for the computation 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 the particular form of the lattice Dirac operator.
If È M is the orthogonal projector to the subspace of fermion fields spanned by the lowest lying eigenmodes of the massive Hermitian Dirac operator D † D with eigenvalues below some threshold value M 2 , the mode number ν(M, µ) can be represented stochastically by: where η 1 , . . ., η N are pseudo-fermion fields added to the theory.
The orthogonal projector È M can be approximated by a rational function of D † D: where M * is a mass parameter related to the spectral threshold value M 1 . The function: is an approximation to the step function θ(−x) in the range −1 ≤ x ≤ 1, where P (y) is in our case the Chebyshev polynomial (of some adjustable degree n) that minimizes the deviation: for some ǫ > 0. Computing the approximation to the spectral projector È M requires solving the following equation an appropriate number of times: for a given source field η. Solving this equation is the main computational cost in the calculation 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.: where the subscript R denotes renormalized quantities. Note that the spectral threshold parameter M renormalizes in the same way as the light quark mass (i.e. M R = Z −1 P M for Wilson twisted mass fermions).
Finally, we give here the relation between the mode number and the mass-dependent renormalized chiral condensate: [31] which is defined to match the chiral condensate to leading order of chiral perturbation theory.

Lattice setup
In this section, we will specify the lattice Dirac operator that is used for our work, i.e. the Wilson twisted mass Dirac operator. For our computations of the chiral condensate, we have 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][36][37] dynamical quarks. The gauge action is: with β = 6/g 2 0 , g 0 the bare coupling and P 1×1 , P 1×2 are the plaquette and rectangular Wilson loops, respectively. For the N f = 2 case, we use the tree-level Symanzik improved action [38], i.e. we set b 1 = − 1 12 , with the normalization condition b 0 = 1 − 8b 1 . In the case of N f = 2 + 1 + 1, we use the Iwasaki action [39,40], i.e. b 1 = −0.331.
The Wilson twisted mass fermion action for the light, up and down quarks for both the N f = 2 and N f = 2 + 1 + 1 cases, is given in the so-called twisted basis by: [41][42][43][44] where m 0,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 −1 P µ. The matrix τ 3 acts in flavour space and χ l = (χ u , χ d ) is a two-component vector in flavour space, related to the one in the physical basis by a chiral rotation. The standard massless Wilson-Dirac operator D W reads: where ∇ µ and ∇ * µ are the forward and backward covariant derivatives. The twisted mass action for the heavy doublet is given by: [43,45] where m 0,h denotes the bare untwisted heavy quark mass, µ σ the bare twisted mass with the twist along the τ 1 direction and µ δ the mass splitting along the τ 3 direction, introduced to make the strange and charm quark masses non-degenerate. The mass parameters µ σ and µ δ are related to the physical renormalized strange m s R and charm m c R quark masses by m s,c R = Z −1 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 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 ≈ 0.05 to a ≈ 0.085 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. The values of the renormalization constant Z P for different ensembles 2 [53][54][55], used to convert bare light quark masses µ and bare spectral threshold parameters M to their renormalized values in the MS scheme (at the scale of 2 GeV), are given in Tab. 3, where we also show the values of r 0 /a (in the chiral limit), used to express our results for the condensate as a dimensionless product r 0 Σ 1/3 . Our physical lattice extents L for extracting physical results   (N f = 2) and β = 1.9, aµ = 0.004 (N f = 2 + 1 + 1).

Results
In this section, we show our results of the calculation of the chiral condensate. First, we illustrate the procedure of extraction of the chiral condensate and discuss the influence of the various errors that enter the computation. Then, we analyze finite volume effects in our simulations. Finally, we move on to our chiral and continuum extrapolations.

Procedure and errors
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, we computed the dependence of the mode number on the renormalized spectral threshold parameter M R for 5 values of M R , from around 2.5 times the renormalized quark mass to around 120 MeV. Shortly above the latter value one starts to see deviations from the linear regime of ν R (M R , µ R ) vs. M R (see Appendix A). Fig. 1 shows the dependence of the mode number on the renormalized spectral threshold 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 according to Eq. (2.7). The error of this slope includes two sources: the error of the slope of the 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 take its value to be the middle point of the chosen fitting interval, see below for details of the fitting intervals considered. Finally, Eq. (2.7) yields, after taking the cubic root 3 : where the first error is the one of the slope ∂ν(M, µ)/∂M and the other one comes from Z P = 0.437 (7) and is dominated by systematic effects -hence, we take it as a systematic error of our computation.  Table 4: Values of r 0 Σ 1/3 extracted from different fitting ranges. Every fit includes at least 3 values of M R . The fit labeled "1 -5" is the full fit. We estimate the error from the choice 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Σ 1/3 can be further converted to a dimensionless product r 0 Σ 1/3 (which will be the final result of this paper, after taking the chiral and continuum limits) or to a physical value in MeV. For the former, we use the value in the chiral limit r 0 /a = 5.35 (4). Since the error of this value is again mostly systematic, we quote it as another systematic error of r 0 Σ 1/3 : where the two errors are as above and the third one comes from r 0 /a. For a conversion to MeV, one needs to choose a value of the lattice spacing in physical units. There are several such estimates for ETMC 2-flavour ensembles, giving for β = 3.9 values including 0.079 fm [6], 0.085 fm [51] and 0.089 fm [58]. Taking this spread into account, the relative error on the lattice spacing is around 7%, which leads to a similar relative uncertainty in the value of the chiral condensate Σ 1/3 in physical units, which amounts to about 20 MeV. This is roughly an order of magnitude more than other errors entering our computation. Even being less conservative and using for the error the value quoted in Ref.
[51] -0.085(2) stat (1) syst fmthe error that it yields is still of the order of 10 MeV. Therefore, we decided to give our final results as the dimensionless product r 0 Σ 1/3 and we chose not to quote any value in MeV 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. Finally, our estimate of r 0 Σ 1/3 for ensemble B40.32, including all sources of error, is: We note that the total error is dominated by systematic errors. This means that increasing statistics would not essentially change our total error. It should be considered an important advantage of the method of spectral projectors that rather moderate statistics (in our case around 230 independent gauge field configurations for this ensemble) leads to a practically negligible statistical error. Let us also mention that the quoted statistical error takes autocorrelations fully into account. We performed an analysis of autocorrelations using two methods and found that in general the autocorrelations are small, even at our smallest lattice spacings. For the details of our autocorrelation analysis, we refer to Appendix B.

Finite volume effects
One of the main sources of systematic effects in Lattice QCD simulations are finite volume effects (FVE). In Ref. [31], theoretical arguments were provided that FVE should be small for the chiral condensate computed from the mode number -with exponentially small difference between the finite volume and infinite volume results and F is the pion decay constant in the chiral limit. Since in practice the mass-dependent chiral condensate is extracted at Λ ≫ µ, the mass M Λ is much higher than the pion mass, which typically governs FVE. Hence, one expects that for the 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 ≥ 2 fm should be small, O( 1%) for M R ≈ O(60 − 120) MeV and renormalized quark masses of O(10 − 20) MeV (with larger FVE at smaller M R ). Indeed, in practice, it was shown in Ref. [31] that for L ≥ 2 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 the mode number and the chiral condensate for: In Fig. 2, we show the volume dependence of the mode number density ν/V for 4-5 different values of the renormalized spectral threshold M R . The mode number density can be computed very precisely. It hence provides a strong test of finite size effects.
The 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 mode number density for L/a = 24 is compatible with the one for L/a = 32 for 3 intermediate values of M R , while it differs by 2-3σ for the lowest and highest M R , thus changing the slope of the mode number vs. M R dependence and the extracted chiral condensate (see the inset of Fig. 2 (left)). This change of slope is statistically significant, but it is still a relatively small, 1-1.5% effect. Taking into account the uncertainty from the choice of the fitting range, the final results for the chiral condensate are compatible for all cases -including the ones for small volumes, indicating that even if the mode number density goes systematically down, the slope of the whole ν(M R , µ R ) dependence is less affected. We have also tried a description of FVE in the framework of the formula derived in Ref. [24]. We conclude that it provides a reasonable agreement with actual lattice data for L/a 24, while FVE for smaller volumes are somewhat underestimated (by a factor of O(2) at L/a = 16, compared to the actually observed FVE).
In the right plot, we show analogous data for N f = 2 + 1 + 1. Similarly, we observe significant 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 reaches a linear lattice extent of around 2 fm (L/a = 24 at both β = 3.9 (N f = 2) and β = 1.9 (N f = 2 + 1 + 1)). Therefore, we used in our analysis all available ETMC ensembles with a linear lattice extent of at least 2 fm 5 .

Chiral and Continuum Limit -N f = 2
We now show our results for the 2-flavour case. For each value of β, we have 2-4 sea quark masses, according to Tab. 1. For each ensemble, we perform computations of the mode number 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 condensate from the slope ν R (M R , µ R ) as function of M R for each ensemble.
The results for r 0 Σ 1/3 for all considered ensembles are gathered in Tab. 5. These results are then used to extrapolate to the chiral limit for each value of β. The chiral corrections to the mass-dependent condensate were calculated at the next-to-leading order of chiral perturbation theory in Ref. [31]. The obtained formula suggests that the mass-dependent condensate 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 a linear extrapolation of the mass-dependent condensate to the chiral limit is fully justified and we follow this conclusion in our chiral extrapolations. As a check, we tried fits of the NLO formula (inserting values of the low energy constants used in Ref. [ Table 5: Results for r 0 Σ 1/3 for all considered N f = 2 ensembles. The given errors are: statistical, resulting from Z P , resulting from r 0 /a, respectively. We also show results in the chiral 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 our errors.
Our extrapolations for all three values of β are shown in the main plot of Fig. 3 (we plot r 3 0 Σ vs. r 0 µ R to allow comparisons between different values of β). The plotted errors are only statistical, since in extrapolations at fixed β, the relative errors from Z P and r 0 /a are the same for all quark mass values (we use chirally extrapolated values of Z P and r 0 /a) -we give them in Tab. 5. To estimate the systematic error originating from the choice of the fitting range in ν R (M R , µ R ) vs. M R fits, we repeated all chiral extrapolations for two tailored fitting ranges -excluding the first value of M R (to account for effects of coming too close to the renormalized quark mass) or the last value thereof (to account for possible deviations 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 latter increases when approaching the continuum limit. The reasons for this behaviour include the fact that the number of quark masses that we use decreases for smaller lattice spacings and at β = 4.2 the slope of the quark mass dependence of the chiral condensate is apparently larger than at coarser lattice spacings 6 , making the final chiral limit value more susceptible to changes in the fitting interval.
Finally, we can use the chirally extrapolated values of the condensate to perform an extrapolation to the continuum limit. We start by discussing the O(a)-improvement of the chiral condensate. even under the R 5 parity transformation: ψ → iγ 5 τ 1 ψ,ψ → iψγ 5 τ 1 [42]. Let us consider the spectral sums [31]: σ k (µ, m q ) = Tr (D † m D m + µ 2 ) −k , where k ≥ 3 for reasons explained 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 the chiral condensate. Representing the spectral sums as a density chain correlation function (for k = 3): it is straightforward to show that the object on the right-hand side is even under R 5 transformation, since the number of densities is even. However, one also needs to consider contact 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 0 ) terms in the mode number -hence they vanish at maximal twist. In this way, the contact terms do not spoil automatic O(a)-improvement of the chiral condensate.
Hence, our continuum limit extrapolation is performed linearly in a 2 , using results at three lattice spacings, with fixed fitting range of the mode number vs. spectral threshold dependence, corresponding to M R ≈ 90 MeV (entering the square root in Eq. (2.7)) for all values of β. As an error, we use the statistical error, combined in quadrature with the error of Z P and r 0 /a. We do not observe significant cut-off effects. The final value in the continuum limit is 0.689 (16). To account for the fitting range error, we perform the full analysis for tailored fitting ranges, excluding the first or last value of M R for each ensemble. This corresponds to a shift in M R to approx. 80 or 100 MeV, respectively. While the extracted value of the condensate in the chiral limit should not depend on the fitting range, in practice the results for different fitting ranges differ, which is due to using only 4-5 values in the fits to extract the slope of ν R (M R , µ R ). The fits from tailored fitting ranges yield 0.678 (18) and 0.718 (20), respectively. To be conservative, as our systematic error from the fitting range we choose the larger difference of the two with respect to the central value 0.689 (16). This finally gives: where the first error is the combined statistical error, the error of Z P and of r 0 /a, while the second error originates from the choice of the fitting range.

Chiral and Continuum
Limit -N f = 2 + 1 + 1 In this subsection, we present results for the N f = 2 + 1 + 1 case. By comparing to the results of the 2-flavour case, we can investigate the role of the dynamical strange and charm quarks.
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 spectral threshold M R , from around 50 to 110 MeV, and extract the mass-dependent condensate from the slope of the ν R (M R , µ R ) vs. M R dependence for each ensemble. We have also computed the mode number for a fifth value of M R ≈ 130 MeV. However, given the very good 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 when this fifth value of M R is included. Because of this, we decided not to include the results at M R ≈ 130 MeV.
In Tab. 6, we show all our results for r 0 Σ 1/3 in the 2+1+1-flavour case. We also include the results of a linear extrapolation to the chiral limit for each value of β, shown in the main plot of Fig. 4. As before, we plot only statistical errors, since all extrapolations are performed at fixed β and the errors from Z P and r 0 /a are the same for all quark masses (we use chirally extrapolated values of Z P and r 0 /a) -given in Tab. 6. Contrary to the N f = 2 case, the slope of the dependence of the mass-dependent condensate on the light quark mass slightly decreases for increasing β (the change of slope is statistically significant when going from β = 1.9 to β = 2.1). This has the effect of lowering the systematic error related to the choice of the fitting range for decreasing lattice spacing 7 , which is for all β  : Results for r 0 Σ 1/3 for all considered N f = 2 + 1 + 1 ensembles. The given errors are: statistical, resulting from Z P , resulting from r 0 /a, respectively. We also show results in the chiral limit, where we also give the systematic error from the choice of the fitting interval (4th error). See text for more details.
the dominating source of error (although for β = 2.1 other errors become comparable). The chirally extrapolated values at three lattice spacings are then used to perform an extrapolation to the continuum limit, which is again compatible with O(a 2 ) cut-off effects. To estimate the fitting range uncertainty, we again perform 3 separate continuum limit extrapolations, using different fitting ranges and different values of M R , corresponding to approx. 80, 90 and 100 MeV. The values of r 0 Σ 1/3 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 full fitting range: (20)(21), with the larger of the differences with respect to values from tailored fitting intervals as the fitting range systematic error. This can be compared to the 2-flavour result which amounts (16) (29) and both results are compatible within errors.

Conclusions
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: which indicates that at the current level of precision, we cannot discriminate the influence of the 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 choice of the fitting range in the dependence of the renormalized mode number on the renormalized spectral threshold. The second most important source of the error is either the statistical error or the error related to the uncertainty in the values of r 0 /a (which are inputs of our analysis). The error from Z P (also an input of our analysis), used to renormalize the quark masses and the spectral threshold parameter, is usually the smallest. However, we want to emphasize that in all cases the fitting range error is the largest one and in most cases it is larger by a factor of 2-4 than any other error. The rather small statistical errors that we obtain 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 values of the spectral threshold M R at which one computes the mode number. This would allow to identify more precisely the linear region of the mode number vs. M R dependence (see discussion in the Appendix) -on the one hand sufficiently far away from the renormalized  [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 0 Σ 1/3 with large-volume continuum limit results found in literature. In the given references, the values of Σ are given in MeV. To convert to the dimensionless product r 0 Σ 1/3 , we combine them with results for r 0 . The first error given is always from the computation of the value in MeV (when there are several errors given, we combine them in quadrature) and the second one from conversion using physical value of r 0 . For RBC-UKQCD, Σ 1/3 = 256(6) MeV, r 0 = 0.487(9) fm [3]. For MILC [4], Σ 1/3 = 278(6) MeV, r 1 = 0.318 (7) fm [4], r 0 /r 1 =1.46(1)(2) [60]. For MILC [5], Σ 1/3 = 281.5(7.9) MeV, r 1 = 0.3133(23) fm [5], r 0 /r 1 =1.46(1)(2) [60]. For S. Borsanyi et al. [9], [61]. For ETMC [6], Σ 1/3 = 269.9(6.5) MeV, r 0 = 0.420 (14) fm. However, newer analyses indicate a higher value of r 0 ≈ 0.45 fm [51]. To take this into account, we added the spread of the new and old value as a systematic error and used r 0 = 0.420 (38) fm to calculate r 0 Σ 1/3 from Σ 1/3 in MeV. For ETMC [27], Σ 1/3 = 299(26)(29) MeV, r 0 = 0.446(9) fm [6]. For HPQCD [28], Σ 1/3 = 283(2) MeV, r 1 = 0.3209(26) fm [62], r 0 /r 1 =1.46(1)(2) [60].
quark mass, on the other hand low enough such that there are no deviations from the linear behaviour at the upper end of the fitting range (observed already at M R ≈ 130 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 collaborations, we attempt in Tab. 7 a comparison. Given the large amount of approaches to compute the chiral condensate, as mentioned in the introduction, we make a selection by only considering results that are given in the literature as continuum limit values from large volume simulations. Note that the other available results that are listed in Tab. 7 are obtained in a different way than the spectral projector method. They mostly use chiral perturbation theory fits to the quark mass dependence of light pseudoscalar meson observables or determine the chiral condensate from the quark propagator. Not discussing the advantages or disadvantages of the various fermion discretizations used, we see in Tab. 7 an overall agreement for the dimensionless quantity r 0 Σ 1/3 , which is reassuring and establishes, in our opinion, the spectral projector method as a valuable alternative to determine the chiral condensate. On the other hand, when looking at the chiral condensate in physical units, see the caption of Tab. 7, a spread of results is obtained. Thus, it seems that the scale setting from the different lattice calculations introduces a systematic effect and it would be desirable to clarify this uncertainty in future more precise calculations.
Acknowledgments We thank the European Twisted Mass Collaboration for generating gauge field configurations ensembles used for this work. We are grateful to A. Shindler for discussions and in particular suggestions concerning O(a) improvement of the mode number. We thank G. Herdoiza for discussions at various stages of this work and for his careful 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 work was 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 made available to us by the Jülich Supercomputing Center, LRZ in Munich, the PC cluster in Zeuthen, Poznan Supercomputing and Networking Center (PCSS). We thank these computer 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 the spectral projectors results for the mode number with ones from explicit computation of 150 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 points correspond to the stochastically evaluated mode number using spectral projectors, while the continuous line (which has an error roughly of the order of the width of the line) is the result from explicitly computing the eigenvalues. We observe very good agreement between the two methods.
Moreover, we used the results of explicit computation of eigenvalues to estimate the region of renormalized spectral threshold of M R where we observe linear dependence between the renormalized mode number and the threshold value (left plot of Fig. 5). The onset of non-linear behaviour corresponds to approx. 130-150 MeV. On the other end of the spectrum, one clearly observes effects of M R close to the renormalized quark mass up to around 10-20 MeV above the latter 8 . This allows us to identify the linear region to extend between around 60 and 120 MeV, to allow for some safety margin. Therefore, we decided to choose our values of M R for the computation of the chiral condensate roughly in this interval. 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 to obtain 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 be relatively 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   time, we computed the mode number for several values of the relative precision of the inverter. The results (for ensemble b40.16) are shown in Fig. 6, which shows that even precision of 10 −2 gives reasonable result. However, we decided to be conservative and we chose a value of 10 −6 for the relative precision of the inverter.
We also checked the dependence of the mode number on the number of stochastic sources used for each gauge field configuration, shown again in Fig. 6. We observe that all results are compatible within error, which matches the suggestion of Ref. [31] that one stochastic source should be enough. However, we observe that adding a second source  might help to considerably reduce the statistical error, which may be important for shorter Monte Carlo runs, when the available number of independent gauge field configurations is rather small. Adding further sources does not change the error considerably, because of correlations 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 the integrated autocorrelation time τ int (in the case of the method of Ref. [64], we quote also the 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 following: • at the smallest value of M , no autocorrelations are observed (τ int compatible with 0.5),  Table 9: Autocorrelations in the mode number, N f = 2 ensembles. We give the number of gauge field configurations used for each ensemble, the step in units of HMC trajectories and the calculated values of τ int using two methods: bootstrap with blocking (boot) and the 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. B25 compared 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.