Topological susceptibility from the twisted mass Dirac operator spectrum

We present results of our computation of the topological susceptibility with $N_f=2$ and $N_f=2+1+1$ flavours of maximally twisted mass fermions, using the method of spectral projectors. We perform a detailed study of the quark mass dependence and discretization effects. We make an attempt to confront our data with chiral perturbation theory and extract the chiral condensate from the quark mass dependence of the topological susceptibility. We compare the value with the results of our direct computation from the slope of the mode number. We emphasize the role of autocorrelations and the necessity of long Monte Carlo runs to obtain results with good precision. We also show our results for the spectral projector computation of the ratio of renormalization constants $Z_P/Z_S$.


Introduction
The topological susceptibility in gauge theories, e.g. in QCD, expresses the fluctuations of the topological charge. As such, it describes non-trivial topological properties of the underlying gauge field configurations. Such properties have far-reaching phenomenological consequences, in particular topological effects are to a large extent responsible for the mass of the flavour-singlet pseudoscalar η ′ meson, making it distinct from the octet of pions, kaons and η. The relation between the topological susceptibility and the η ′ mass is expressed in the Witten-Veneziano formula [1,2].
There exist many definitions of the topological charge on the lattice 1 and there has been a debate in the literature about the validity of different approaches. One of the main problems is the appearance of non-integrable short distance singularities in some definitions, which require regularization.
To avoid such theoretical problems, a possible solution is to use the definition of the topological charge as the index of the overlap Dirac operator [4], which is by construction integer-valued. However, this is very demanding in terms of computing time and hence impractical when large lattice sizes are used. Using Ginsparg-Wilson fermions, it is also possible to derive an expression for the topological susceptibility which does not have any power divergences [5,6]. This has been further generalized by Lüscher, leading to a definition employing the so-called density chain correlation functions [7]. The latter can be evaluated efficiently using the method of spectral projectors [8]. This definition of the topological susceptibility was subject to numerical analysis in the quenched case [9] and it is the aim of the present paper to analyze the results of its usage in the case with N f = 2 and N f = 2 + 1 + 1 active flavours of twisted mass fermions.
The outline of the paper is as follows. In section 2, we describe the theoretical principles of the adopted approach. Section 3 presents our lattice setup. In section 4, we show our results for the renormalization constants ratio Z P /Z S and in section 5 for the topological susceptibility. We conclude in section 6. In an appendix, we show our tests concerning the number of stochastic sources.

Theoretical principles
The method that we follow in this paper was introduced in Refs. [8,9] and we refer to these papers for a comprehensive description. Here, we summarize only the main points to render the paper self-contained.
Let us define an orthogonal projector È M to the subspace of fermion fields spanned by the lowest lying eigenmodes of the operator D † D with eigenvalues below some threshold value M 2 . In practice, if the projector È M is approximated by a rational function of D † D, denoted by Ê M (see Refs. [8,9] for the details of this approximation), the following equation for the topological susceptibility χ holds: The calculation of the topological susceptibility from this expression requires an evaluation of three gauge field ensemble averages. However, if the value of the scheme-and scaleindependent ratio Z P /Z S is available from another computation, the above expression can be rewritten as: where the numerator can be expressed using two stochastic observables defined in Ref. [9]: where N is the number of randomly generated pseudofermion fields η i added to the theory 2 and (2.4) 2 We use the Z(2) random noise, i.e. (ηi)r = (±1 ± i)/ √ 2, where r spans the set of source degrees of freedom (space-time, colour, spin) and all signs ± are chosen randomly. (2.5) The term B /N is a correction to the result given by C 2 needed if the number of stochastic sources N is finite and if one computes: using the same stochastic sources for C {η k } and C {η l } . In chiral symmetry preserving formulations of Lattice QCD (e.g. using overlap fermions), the observable C is just the index Q of the Dirac operator, i.e. the difference in the number of zero modes with positive and negative chirality, since (η, γ 5 η) = ±1 if η is a zero mode and 0 otherwise. Moreover, in such theories Z P = Z S and in the limit N → ∞ Eq. (2.3) becomes just the well-known formula χ = Q 2 /V . The distribution of Q is expected to be of the Gaussian type (with Q = 0) and the topological susceptibility is then alternatively given by the width of this distribution. In theories where chiral symmetry is explicitly broken at finite lattice spacing, e.g. for Wilson fermions, the observable C is in general non-integer and counts the number of zero modes only approximately (up to cut-off effects). However, as we will show, C is still compatible with a Gaussian-shaped distribution and the renormalized C ren ≡ Z S Z P C can be thought of as a proxy for the topological charge. As it is well known, the topological charge is an observable which is particularly susceptible to autocorrelations in Monte Carlo (MC) time [10]. Hence, to obtain reliable estimates of the topological susceptibility, it is essential that MC histories are long enough, such that all topological sectors are correctly probed.
Since the observable C is strongly related to the topological charge, its autocorrelation time and the quality of its distribution provides a criterion of MC history being "long enough". In particular, we demand the distribution of C to be compatible with a Gaussian and C should be compatible with zero.
We have mentioned above that the full renormalized topological susceptibility can be obtained from expression (2.1). This means that the ratio of renormalization constants Z P /Z S can be calculated with spectral projectors, as first noticed in Ref. [8]. The formula reads: where B is given by Eq. (2.5) and A is: i.e. it is the mode number ν(M ) -the number of eigenmodes of the operator D † D with eigenvalues below the threshold value M 2 .
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 twisted basis by: [20][21][22][23] where τ 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. m 0 and µ l are the bare untwisted and twisted light quark masses, respectively. The renormalized light quark mass is µ R = Z −1 P µ l . 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: [22,24] where µ σ is the bare twisted mass with the twist along the τ 1 direction and µ δ the mass splitting along the τ 3 direction that makes the strange and charm quark masses non-degenerate.
The physical renormalized strange m s R and charm m c R quark masses are related to the bare parameters µ σ and µ δ via m s R = Z −1 P (µ σ − (Z P /Z S )µ δ ) and m c R = Z −1 P (µ σ + (Z P /Z S )µ δ ). The heavy quark doublet in the twisted basis χ h = (χ c , χ s ) is related to the one in the physical basis by a chiral rotation.
The details of the gauge field ensembles considered for this work are presented in Tab. 1 for N f = 2 and Tab. 2 for N f = 2 + 1 + 1. They include lattice spacings from a ≈ 0.045 fm 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 3 [33,34,36] Parameters of the N f = 2 gauge field ensembles [11][12][13]. We show the inverse bare coupling β, lattice size (L/a) 3 × (T /a), bare twisted light quark mass aµ l , renormalized quark mass µ R in MeV, critical value of the hopping parameter at which the PCAC mass vanishes and physical extent of the lattice L in fm and the product m π L.
Ensemble β lattice  Table 2. Parameters of the N f = 2 + 1 + 1 gauge field ensembles [14][15][16]. We show the inverse bare coupling β, lattice size (L/a) 3 × (T /a), bare twisted light quark mass µ l , renormalized quark mass µ l,R in MeV, critical value of the hopping parameter at which the PCAC mass vanishes, physical extent of the lattice L in fm and the product m π L.
quark masses µ l 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. There we also show the values of r 0 /a (in the chiral limit), used to express our results for the topological susceptibility as a dimensionless product r 4 0 χ. Our physical lattice extents L for extracting physical results range from 2 fm to 3 fm (in the temporal direction, we always have T = 2L). To check for the size of finite volume effects, we included different lattice sizes for β = 3.9, aµ l = 0.004   [30][31][32], the schemeand scale-independent renormalization constants ratio Z P /Z S and the renormalization constant Z P in the MS scheme at the scale of 2 GeV [33][34][35][36][37], for different values of β and N f = 2 and N f = 2 + 1 + 1 flavours.

Results -Z P /Z S
We first present our results for the renormalization constants ratio Z P /Z S , which is a scaleand scheme-independent quantity. Nevertheless, in order to avoid problems with e.g. cutoff effects or dependence on the threshold parameter We perform our N f = 2 analysis using small volume ensembles (b40.16, c30.20, d20.24 and e17.32) at a fixed pion mass of around 300 MeV (in infinite volume). In this way, we can keep the computational cost rather low and at the same time investigate a wide range of values of M R to control the systematic effects of varying M R . For all these ensembles, we can compare the values of Z P /Z S with an alternative computation -in the framework of the RI-MOM renormalization scheme (β = 3.9, 4.05, 4.2 [34]) or the X-space renormalization scheme (β = 4.35 [35]). The dependence of Z P /Z S on M R is shown in Fig. 1. For small values of M R , we observe a significant dependence of Z P /Z S on the threshold parameter M R . For larger values of M R , the dependence of Z P /Z S on M R flattens and we observe a tendency to approach a plateau. This signals that we obtain the above discussed lattice window, where we can extract the scale-and scheme-independent value. However, the lattice data shows that the strong dependence of Z P /Z S on M R below about 1 GeV is related to the value of the lattice spacing a involved. In addition, also the size of variation of Z P /Z S for 1 GeV < M R < 2 GeV is getting smaller for decreasing values of the lattice spacing. For β = 3.9, the change in Z P /Z S when going from around 1 to 2 GeV is approx. 6%, while for finer lattice spacings this change decreases to approx. 3%, 2% and below 1% (for β = 4.05, β = 4.2 and β = 4.35, respectively). A reassuring observation is, however, that the plateau value is consistent with values from the RI-MOM or X-space renormalization schemes, shown in Fig. 1    remark also that at large values of M R , finite volume effects are small. We have explicitly checked that the gauge ensemble averages of the observables A and B are always compatible with each other for ensembles b40.16 and b40.24, provided that M R 1 GeV. Moreover, the ratio of these two quantities, which gives Z P /Z S , is compatible between the two ensembles at all values of M R that we investigated -see Fig. 2. Therefore, the conclusions from our small volume results are valid in general.
To summarize, the method of spectral projectors allows in principle a computation of the ratio of Z P /Z S . However, to obtain the universal scale-and scheme-independent value, the calculation of the observables A and B (with N = 1 stochastic source) has to be performed at a rather large number of threshold parameters M R to be able to explore the significant M R -dependence we observe. To account for this M R -dependence, we followed the strategy to take the central value of Z P /Z S at some fixed (in physical units) value of M R , e.g. 1.5 GeV (sufficiently far away from the low-energy scales and sufficiently below the inverse lattice spacing for typical parameters of contemporary simulations), and assign a systematic error related to the difference of Z P /Z S across a range of scales (e.g. M R between 1 and 2 GeV). If we follow this strategy, we obtain the results of Tab. 4 and the horizontal bands in Fig. 1. The first given error of the spectral projectors result is statistical and the second one comes from the residual M R -dependence of Z P /Z S . Note that only the RI-MOM results, given for comparison in Tab. 4, were chirally extrapolated. However, non-zero quark mass corrections to the chiral limit value were found to be small in our setup, both in the RI-MOM scheme and in the X-space scheme [33][34][35]. In particular, the values of Z P /Z S in the chiral limit and at the pion mass of 300 MeV are always compatible in both of these schemes. Given our experience that Z P /Z S has a very mild quark mass dependence [33][34][35], we consider the values obtained here at a fixed but small pion mass to be an appropriate estimate. The overall agreement between the spectral projector method and other renormalization schemes is certainly reassuring. However, it would be very good to understand the M R -dependence of Z P /Z S better and to disentangle effects that lead to it. For example, a lattice perturbative calculation within the framework used would be very helpful to learn about the role of cut-off effects.
4.2 N f = 2 + 1 + 1 We repeated the computation of the M R -dependence of Z P /Z S also for one chosen ensemble with N f = 2 + 1 + 1 flavours (B55.32). The chiral limit value from RI-MOM is 0.697(7) [36]. The residual M R -dependence originating from spectral projectors is rather large in this case and the prescription from the previous subsection leads to the value 0.637(1) (21). The systematic error is comparable to the one for β = 3.9 with N f = 2, which corresponds β Z P /Z S (spec.proj.) Z P /Z S (RI-MOM) Z P /Z S (X-space) 3 Table 4. The values of the scale-and scheme-independent ratio Z P /Z S for N f = 2 ensembles, extracted from spectral projectors (the first error given is statistical and the second one systematic from varying the threshold value M R ), as compared to RI-MOM [34] and X-space results [35]. All RI-MOM results and the X-space results at β = 3.9 and β = 4.05 were chirally extrapolated.  to a similar lattice spacing. Although there is some tension between the spectral projector result and RI-MOM, the observed difference between the two results is still plausible, given a finite and rather large value of the lattice spacing. Note, however, that in the following we do not rely on the values of Z P /Z S from spectral projectors -we rather use the RI-MOM values to evaluate the topological susceptibility.

Results -topological susceptibility
In this section, we discuss our results for the topological susceptibility. We first show the details of our analysis for two of the 2+1+1-flavour ensembles -B55.32 and B75.32. Then, we investigate finite volume effects and finally we present results for the cases with N f = 2 and N f = 2+1+1 flavours of twisted mass fermions and perform chiral perturbation theory fits to the quark mass dependence of the topological susceptibility.

Examples -ensemble B55.32 and B75.32
We start with the ensemble B55.32, see Tab. 6, for which we performed measurements on 538 independent gauge field configurations separated by 40 MC trajectories, using N = 6 stochastic sources for each configuration. For a discussion about the optimal number of stochastic sources per configuration, we refer to Appendix A. The MC history of the observable C (whose fluctuations determine the topological susceptibility) is shown in the left panel of Fig. 3. We observe that different topological sectors are sampled and the magnitude of fluctuations seems to be rather uniform for different regions of MC time. As we have stated above, the sampling is correct if the histogram of C is close to Gaussian and if the ensemble average C = 0. The histogram of the observable C for the ensemble B55.32 is shown in Fig. 4 (left). It is almost ideally symmetric and it is almost perfectly Gaussian. We have therefore fitted the following Gaussian ansatz: where N is a normalization constant and σ is related to the topological susceptibility: χ = (Z S /Z P ) 2 (σ 2 − B /N ), i.e. σ 2 = C 2 . The 3 fitting parameters are then: N , C and σ. There is very good agreement between C extracted from the histogram and computed directly by averaging -the former yields 0.02 (20) and the latter -0.06 (16), which implies that both the negative and positive topological charge sectors are sampled equally often. The bare topological susceptibility extracted from the direct computation and using σ 2 from the fit of the histogram is 3.56(51) · 10 −6 (histogram) and 3.65(33) · 10 −6 (direct). This agreement implies that indeed the observable C is Gaussian distributed and can be interpreted to play the role of the topological charge. We also note that the constructed histograms and the extracted values of C and a 4 χ depend very little on the chosen bin size. Using bin sizes of 0.5, 1, 2 and 3, our results for a 4 χ change only by a few percent and are fully compatible within errors. Hence, we decided to use such bin size that the number of bins with non-zero number of gauge field configurations is around 10. We emphasize that the good properties of the histogram ( C ≈ 0 and Gaussian shape) hold only if the MC history is long enough. We think that both properties can provide a good benchmark whether the topological charge sectors are sampled in a correct way. The statistics that we have for ensemble B55.32 is significantly higher than for other ensembles. Let us show the details for a more typical ensemble B75.32 with around 100 independent measurements. The Monte Carlo history (right panel of Fig. 3) indicates a correct sampling of topological sectors, however it is not long enough to build a fully symmetric histogram (Fig. 4 (right)). For example, the number of configurations for which −3 ≤ C < −1 and 1 ≤ C < 3 is, respectively, 47(7) and 31 (5), where the error comes from bootstrap with blocking analysis and takes into account autocorrelations. Hence, in the generated ensemble, the samples with slightly negative topological charge are somewhat overrepresented with respect to the ones with slightly positive topological charge, although statistically they are still compatible. As a consequence, the peak of the Gaussian fit is for C below zero. Nevertheless, within the computed errors we observe that the shape of the histogram is close to Gaussian and the topological susceptibility and C are within large errors compatible between the fit and the direct computation and read for ensemble B75.32: C = 0.04(35) (direct) and -0.20(37) (histogram), bare topological susceptibility: a 4 χ = 4.13(48) · 10 −6 (direct), a 4 χ = 4.80(1.10) · 10 −6 (histogram). However, we would like to give a warning that the rather low statistics we have for the ensemble B75.32 may lead to an underestimation of the error, i.e. the error of the error might be large. To reach full confidence for the obtained results, statistics of at least the size we have for the ensemble B55.32 would be necessary.

Finite volume effects
Before we show results for all our ensembles, we shortly discuss finite volume effects (FVE) in our simulations. We show the bare topological susceptibility for three N f = 2 + 1 + 1 ensembles at β = 1.90, aµ l = 0.004 and four N f = 2 ensembles at β = 3.90, aµ l = 0.004 in Fig. 5. All ensembles give compatible results (with some tension between A40.20 and A40.32), but given the precision we have for these ensembles, i.e. statistical errors of the order of 10-20%, we can not conclude about the size of FVE from numerical data. However, general arguments involving the size of FVE imply that they should be exponentially small if m π L 4 (see e.g. Ref. [39]). Since this condition is satisfied for almost all of our N f = 2 + 1 + 1 ensembles (see the last column of Tab. 2), we are confident that FVE are much smaller than our statistical errors. For N f = 2, we analyze the quark mass dependence of the topological susceptibility only at one lattice spacing (β = 3.9) and the product m π L > 4 for all of them. However, even in a small volume (L ≈ 1.3 fm, with m π L ≈ 2.4), FVE are not larger than the statistical errors (cf. β = 3.9, L/a = 16 and L/a = 32 in Fig. 5).

N f = 2 results
Tab. 5 provides our results for the observables A , B , C and the topological susceptibility in the case of N f = 2 flavours. In Fig. 6, we show our results at a single lattice spacing corresponding to β = 3.9 and a physical volume such that the condition m π L > 4 is satisfied. In order to test whether the obtained values of the topological susceptibility could, in principle, be used to obtain a value for the chiral condensate, we apply the leading order (LO) Chiral Perturbation Theory (χPT) expression for N f flavours of light quarks:  where Σ is the chiral condensate. In particular, we impose that the topological susceptibility vanishes at zero quark mass. Working with the assumption that LOχPT can be applied, the slope of this fit gives the following result for the renormalized condensate (MS scheme at 2 GeV): where the error is mostly statistical, but takes into account also the uncertainties of Z P /Z S , Z P and r 0 /a. The error decomposition is as follows: r 0 Σ 1/3 = 0.650(21)(6)(3) (5), where the first error is statistical, the second comes from the uncertainty of Z P /Z S (entering via Eq. (2.3)), the third one from the uncertainty of Z P (entering the renormalized quark mass µ l ) and the fourth one from the final conversion of aΣ 1/3 to r 0 Σ 1/3 . The used values of Z P /Z S , Z P and r 0 /a, together with their uncertainties, are shown in Tab. 3. The final error quoted is the sum of the individual errors, combined in quadrature. We recall here respective values from our direct determination from the mode number of the Dirac operator: 0.696(20) (at β = 3.9 in the chiral limit) or 0.689(33) (in the continuum limit and in the chiral limit) [40]. The fact that we observe an agreement indicates a posteriori the validity of our assumption about the applicability of LOχPT, at least within the large errors of our present results 4 . Figure 6. Renormalized quark mass dependence of the renormalized topological susceptibility (normalized with r 4 0 ) for N f = 2 ensembles at β = 3.9. The fit is to a LO χPT expression, 16 Table 5. Our results for N f = 2 flavours. We give the ensemble label, the number of stochastic sources N , the number of configurations used (cnfs), the step between measurements (in units of molecular dynamics trajectories) and the values of A , B , C and the topological susceptibility, together with integrated autocorrelation times τ int (in units of measurements). The error for r 4 0 χ is, respectively, statistical, resulting from the uncertainty of Z P /Z S (from the RI-MOM method) and resulting from the uncertainty of r 0 /a. In all other cases the error is statistical only.

N f = 2 + 1 + 1 results
In this subsection, we discuss our data for the case with 2+1+1 active flavours. Our results for the observables A , B , C and the topological susceptibility are collected in Tab. 6 and Fig. 7 shows the results for the topological susceptibility. We required that the autocorrelations for the topological charge are kept under control, i.e. can be measured with reasonable accuracy using the method proposed in Ref. [41] (UW method). This method allows for an estimate of the integrated autocorrelation time τ int and of its error. We also made an independent error analysis using the method of bootstrap with blocking. In all cases, we found results compatible with the UW method, given in Tab. 6. In particular, we  Our results for N f = 2 + 1 + 1 flavours. We give the ensemble label, the number of stochastic sources N , the number of configurations used (cnfs), the step between measurements (in units of molecular dynamics trajectories) and the values of A , B , C and the topological susceptibility, together with integrated autocorrelation times τ int (in units of measurements). The error for r 4 0 χ is, respectively, statistical, resulting from the uncertainty of Z P /Z S (from the RI-MOM method) and resulting from the uncertainty of r 0 /a. In all other cases the error is statistical only.  found that the autocorrelation time for the observable C is τ int 1. Typically, we have O(200) configurations per ensemble, although for our ensembles at the finest lattice spacing, we only have around 100 configurations. Thus, the histograms that we can build have large statistical errors and within these large errors the deviation from a zero-centered Gaussian is insignificant. Few exceptions to this rule occur -e.g. for ensemble A80.24 C is more than 3σ away from zero. The typical error of the computed topological susceptibility is of the order of 15% and we manage to go below 10% only for ensemble B55.32. In this way, we conclude that the precision one can reach for the topological susceptibility is only modest. However, we want to emphasize here that this is a consequence of too short lengths of typical Monte Carlo simulations in Lattice QCD and do not originate from the spectral projector method itself. Especially with finer lattice spacings, autocorrelations are such that to obtain truly independent gauge field configurations one has to perform measurements skipping several trajectories. Our experience shows that to obtain a 10% precision in the computation of the topological susceptibility, we need around 300-400 truly independent configurations, which implies Monte Carlo runs of 10000-20000 trajectories, which is somewhat longer than is typically needed for most other applications.
In order to describe the quark mass dependence of the topological susceptibility, we follow the same strategy as discussed above for N f = 2 flavours. Using only the LOχPT formula, we decided to apply a mass cut on our data, excluding points for which the pion mass is larger than 400 MeV, i.e. keeping points for which r 0 µ R < 0.07. The fits of the LO formula to our data are shown in Fig. 8(a,b,c).
As in the N f = 2 case, we can extract the chiral condensate from the dependence of χ on the quark mass. We have performed an analysis separately for each lattice spacing, taking the individual uncertainties of Z P /Z S , Z P and r 0 /a into account (in the way described in the previous section) to propagate them to the values of r 0 Σ 1/3 at finite lattice spacings (shown in Fig. 8(d)). Such obtained values are then extrapolated to the continuum limit, yielding the value: We mention here that it is possible to prove that the topological susceptibility computed using the spectral projector method and twisted mass fermions at maximal twist is O(a)improved. This is not guaranteed a priori by standard arguments for the automatic O(a)improvement at maximal twist [21], since the topological susceptibility is defined via density chains that include integrals (in the continuum) or sums (on the lattice) over all space time points, which leads to contact terms with short distance singularities. Such contact terms can, in principle, spoil automatic O(a)-improvement. The proof that this is not the case is sketched in Refs. [42,43], while for the details of this proof we refer to an upcoming publication [44]. In general, the quality of our LOχPT fits is reasonable (see the values of χ 2 /d.o.f. in the caption of Fig. 8), with the exception of β = 1.95, for which χ 2 /d.o.f. ≈ 5. This may signal the presence of effects beyond the ones captured in our LOχPT fitting ansatz. However, with our current precision we are not able to address this issue. The fact that some data points are off the fit line may well be a statistical fluctuation at this level of precision. As a check of the robustness of our result, we performed also another LOχPT fit including all our data, i.e. also pion masses between 400 and 500 MeV. This leads to a value for the chiral condensate in the continuum limit: r 0 Σ 1/3 = 0.619(58). The result from this additional fit is slightly lower, although still compatible with the one from fits applying a mass cut. The values of χ 2 /d.o.f. for the LOχPT fits without pion mass cuts are: 1.70 (β = 1.9), 4.52 (β = 1.95), 0.49 (β = 2.1), i.e. they are comparable to the ones for fits without pion mass cuts (see the caption of Fig. 8).
The error that we give is dominated by statistical uncertainties, but the contribution from the systematic errors related to r 0 /a and Z P /Z S is also included. However, it does not include the main source of systematic effects coming from χPT: the use of only the leading order expression. As we mentioned above, our precision is not enough to use an NLO fitting ansatz. Still, our result is in agreement with the direct determination from the mode number on the same set of gauge field ensemblesr 0 Σ 1/3 = 0.680(20)(21) [40], indicating that LOχPT describes the quark mass dependence of the topological susceptibility at least within the rather large errors of our results.
It is worth emphasizing that at β = 1.9 and β = 1.95 the data for r 4 0 χ do not show a clear tendency to assume a zero value when the quark mass is decreased. Only at β = 2.1 and hence closer to the continuum limit, the data seem to approach zero linearly in the quark mass. Thus, in order to cleanly identify this expected behaviour of the topological susceptibility, smaller quark masses and a significantly increased precision are required.

Conclusions
We have computed the topological susceptibility in dynamical Lattice QCD simulations using the method of spectral projectors. This method has two important advantages that we want to emphasize here: • it relies on a theoretically sound definition of the topological susceptibility from density chain correlators that is free of short distance singularities, • it is significantly less computer time expensive than the topological susceptibility computation from the index of the overlap Dirac operator.
One main result of our work is that the topological susceptibility is affected by substantial statistical fluctuations necessitating long Monte Carlo histories. With typical parameter values of Lattice QCD simulations nowadays, i.e. lattice spacings of 0.05 fm a 0.1 fm and lengths of Monte Carlo runs of O(5000) trajectories with autocorrelation times τ int = O(10) trajectories, it is very difficult to obtain errors smaller than 10-15% for a given ensemble. We emphasize that this is not a property of the method used here, but of the gauge field configurations themselves and as such can not be easily overcome, i.e. without running very long simulations. In addition, the topological properties of gauge fields -here characterized by the quantity C of Eq. (2.4), which is closely related to the topological charge -tend to be particularly susceptible to autocorrelation effects, which increase with decreasing lattice spacing. This is indeed observed with the present method and implies that very high statistics is needed (in particular at small lattice spacings) to overcome this problem, unless one works with open boundary conditions that naturally allow to move the problem to at present unachievably small lattice spacings [10]. Despite these difficulties, we were able to demonstrate that by imposing LO chiral perturbation theory as a description of our data for the topological susceptibility, values of the chiral condensate could be determined, which read: r 0 Σ 1/3 = 0.650 (22) (N f = 2, no continuum extrapolation) and r 0 Σ 1/3 = 0.651(61) (N f = 2 + 1 + 1). These results, although having large errors for the reasons discussed above, are fully compatible with the ones of our direct calculation using spectral projectors [40]. We estimate that a meaningful test of the NLO chiral perturbation theory prediction for the quark mass dependence of the topological susceptibility would require a factor 3-10 longer runs (than typical ones, as specified above), which would bring the errors down below 10%. Nevertheless, we have shown that such calculations are becoming feasible with present-day computing resources and the advantages of computations with a theoretically sound definition of the topological susceptibility using density chains, promote the here used method to one of the most promising ways to address topological properties of QCD in the future.