Topological susceptibility from twisted mass fermions using spectral projectors and the gradient flow

We compare lattice QCD determinations of topological susceptibility using a gluonic definition from the gradient flow and a fermionic definition from the spectral projector method. We use ensembles with dynamical light, strange and charm flavors of maximally twisted mass fermions. For both definitions of the susceptibility we employ ensembles at three values of the lattice spacing and several quark masses at each spacing. The data are fitted to chiral perturbation theory predictions with a discretization term to determine the continuum chiral condensate in the massless limit and estimate the overall discretization errors. We find that both approaches lead to compatible results in the continuum limit, but the gluonic ones are much more affected by cut-off effects. This finally yields a much smaller total error in the spectral projector results. We show that there exists, in principle, a value of the spectral cutoff which would completely eliminate discretization effects in the topological susceptibility.

(Dated: July 9, 2018) We compare lattice QCD determinations of topological susceptibility using a gluonic definition from the gradient flow and a fermionic definition from the spectral projector method.
We use ensembles with dynamical light, strange and charm flavors of maximally twisted mass fermions. For both definitions of the susceptibility we employ ensembles at three values of the lattice spacing and several quark masses at each spacing. The data are fitted to chiral perturbation theory predictions with a discretization term to determine the continuum chiral condensate in the massless limit and estimate the overall discretization errors. We find that both approaches lead to compatible results in the continuum limit, but the gluonic ones are much more affected by cut-off effects. This finally yields a much smaller total error in the spectral projector results. We show that there exists, in principle, a value of the spectral cutoff which would completely eliminate discretization effects in the topological susceptibility.

I. INTRODUCTION
Topology is essential in understanding the physics of QCD. Due to the fact that the lattice inherently has different topology from the continuum, an accurate and precise measurement of topological charge is a difficult task in lattice QCD, with a long history of failed attempts and practical and theoretical issues to overcome [1][2][3][4][5][6]. Yet, a good definition of topological charge is essential to modern lattice calculations, e.g. ones involving CP-violation. In particular, the neutron electric dipole moment (nEDM) depends on a topological charge dependent parameter called α 1 , defined in Eq. (31) of [7]. A significant source of error in the nEDM calculation and the calculation of the mass of the η meson, is due to uncertainty in the topological charge [8].
There are several definitions of topological charge that belong to two broad classes, fermionic and gluonic definitions. In the continuum, they are all equivalent, but on the lattice they can have very different discretization effects. The gluonic definition, utilizing various smearing techniques to damp UV fluctuations, such as the gradient flow (GF) [9][10][11] that we use in this paper, is computationally relatively cheap to calculate and is currently widely used, for example, in the aforementioned nEDM calculation [8]. In this paper we focus on a fermionic definition using the method of spectral projectors [2,12], but we also compare the results to ones from a GF-smeared gluonic definition.
Given the topological charge, Q, on every gauge field configuration, the topological susceptibility, χ, is given by the ensemble average Q 2 , divided by the lattice volume V . Hence, it is essentially a quantity which provides a measure of topological charge fluctuations. Using leading order chiral perturbation theory ( χ PT) and parametrizing the leading discretization term, the topological susceptibility can be written as a simple equation depending on the quark mass and lattice spacing.
In this paper, the coefficient of the mass term, which is a low energy constant (LEC) known as the chiral condensate, Σ, is also extracted, from both the spectral projector definition and the GF-smeared gluonic definition. Our aim is, in particular, a quantitative investigation of cut-off effects present in both formulations, motivated by our observations that the gluonic topological charge is subject to large discretization effects, at least in the setup of our lattice study (see below). This paper is organized as follows: Section II presents the lattice setup. In Section III, we briefly review the spectral projector method and the gluonic definition of the topological charge (with GF as the smoother), as well as the χ PT formula for topological susceptibility. Section IV contains results from the two methods, together with fits used to determine the continuum value of the chiral condensate. Lastly, we conclude in Section V reviewing our results and discussing future directions.

II. LATTICE SETUP
The following calculations were performed using twisted mass configurations with N f = 2+1+1 dynamical flavors provided by the European Twisted Mass Collaboration (ETMC) [13][14][15]. The Iwasaki action [16] is used in the gauge sector. The action in the fermionic sector for the light flavors (in the twisted basis) reads: [17][18][19][20] where τ 3 is the third Pauli matrix, acting in flavor space, χ l = (χ u , χ d ) is a vector in flavor space and D W is the standard Wilson-Dirac operator. The bare mass parameters, m 0 and µ l , are the untwisted and twisted light quark masses, respectively. The renormalized light quark mass, which we denote by µ l,R , is related to the bare twisted mass µ l,R = Z −1 P µ l , where Z P is the renormalization function of the pseudoscalar density.
The action for the heavy flavors is: [18,20] where µ σ is the bare twisted mass with the twist along the τ 1 direction in charm-strange flavor space and µ δ gives the splitting along the τ 3 direction, i.e. generates unequal strange and charm quark masses. χ h = (χ c , χ s ) is a vector in charm-strange flavor space, written in the twisted basis.
Three different lattice spacings were used, ranging from ∼ 0.062 − 0.089 fm [25]. For each spac-  Tables I and II, while in Table III

III. THEORY
The method of spectral projectors was introduced and described in detail in Refs. [2,12]. In these papers, the projector was constructed stochastically, thus limiting the computational cost with respect to an explicit computation of eigenmodes from O(V 2 ) to O(V ), at the expense of introducing stochastic noise. For the current investigation, we opted for an explicit computation of eigenmodes, in order not to have this stochastic contamination. In practice, we used the ARPACK (ARnoldi PACKage) eigensolver library [28] to calculate the lowest 400 eigenmodes of the normal massless Wilson-Dirac operator, D † W D W . We performed the calculations on P100 Nvidia GPU architecture, using the ETMC's code library [29], which uses the QUDA linear solver library as its base [30].
Here, we summarize the relevant details for calculating the topological susceptibility from the Hermitian Wilson-Dirac operator (D † W D W ) spectrum. The bare topological charge can be defined and the bare topological susceptibility as The renormalized quantities are In the above formulae, u i are eigenvectors of the normal operator D † W D W and M 2 is the renormalized spectral threshold, i.e. all the renormalized eigenvalues taken into account have magnitude smaller than M 2 . The 400 lowest eigenmodes allow M up to 160 MeV for all ensembles described in Table I. For chiral fermions, such as overlap fermions, the topological charge receives contributions only from exact zero modes of D † W D W . The zero modes correspond to λ = 0, such that R λ=0 = ±1 and R λ>0 = 0. This would indicate that a spectral threshold slightly larger than zero would be enough to obtain a topological charge exactly equivalent to the one from the index of the overlap Dirac operator [31][32][33]. For non-chiral fermions, such as Wilson twisted mass fermions which we use in this paper, the would-be zero and non-zero modes are shifted by O(a 2 ) effects. Moreover, different values of the spectral threshold may correspond to different cut-off effects. An optimal value of M , i.e. one minimizing discretization effects, can be chosen by analyzing several values of the lattice spacing. This allows one to find such a value that the slope of the topological susceptibility vs. the lattice spacing is minimal, as we show explicitly below. The value that we find is appropriate for the here considered setup and it may be different when using e.g. a different fermion discretization.
If one works with a single ensemble, it is, of course, not possible to find such an optimal value. However, one can still check the dependence of the topological susceptibility on M and if a region of small variation of the susceptibility is found, it likely corresponds to small cut-off effects, as we also observe in our data. Note that the topological susceptibility from the spectral projector method is automatically O(a)-improved when using twisted mass fermions at maximal twist [34].
The other definition of the topological charge that we employ in this paper is the gluonic (field theoretic) one. The calculation of the topological charge with this definition proceeded along the lines of Ref. [5], which we refer to for more details. The field theoretic definition reads: The definition of the field strength tensor used is tree-level Symanzik improved, i.e. O(a 2 ) improved by including both clover terms and rectangular Wilson loops, as described in Ref. [5]. UV fluctuations are filtered using the GF procedure, as described in Ref. [9]. The smoothing action in the GF is the tree-level Symanzik improved gauge action.
Next, we derive an expression for the topological susceptibility in chiral perturbation theory, following the arguments used in Ref. [35] but adapted for twisted mass. First, we consider the partition function of QCD with non-zero Θ QCD : The topological susceptibility can be written as the second derivative of the energy density, evaluated at Θ = 0, The lowest order, two flavor χ PT Lagrangian can be written as: where U (x) is an element of SU (2) that can be expressed as wheren is a unit vector in su (2)  The partition function can now be written in terms of L (2) χ PT : For volumes that are large relative to the quark mass, V µ l Σ 1, the group integral is dominated by U (x) and by extension φ(x), which is x independent and minimizes −L (2) χ PT . For maximal twist and non-zero Θ, this minimum occurs at φ = ±π/2,n = (0, 0, 1). The partition function simplifies to: which yields: which is correct to leading order in the quark mass. In the lattice calculation, there will be corrections of order a 2 , as maximal twist removes any O(a) effects.

IV. EXTRAPOLATION OF THE CHIRAL CONDENSATE
A. Three-parameter fit One approach for determining the chiral condensate is a global 3-parameter fit to data at different values of the quark masses and the lattice spacing. Specifically, we fit to the expression: where r 0 is the Sommer parameter, used to make all quantities dimensionless, Σ is the renormalized chiral condensate, and c and α are the coefficients of the O(a 2 ) and O(µ l a 2 ) discretization terms.
The square brackets indicate the actual fitting parameters. Other expressions were tested, including a two parameter fit, where α was set to zero. Ultimately, the three parameter fit in Eq.  Table IV.
The first observation from Fig. 3 is that a consistent value of Σ 1/3 is obtained for any value of M in the spectral projector method, with a small tendency for a smaller value for M below around 50 MeV. The value from the gluonic GF definition is consistent with all the spectral projector values.
The statistical error from the spectral projector definition is a factor of 2-4 smaller than the one from the gluonic definition. It should be noted that the statistical error noticeably increases for M 170 MeV due to fewer ensembles being included in the fits.
We also observe that the spectral projector definition produces a topological susceptibility with   Table I. Resulting fit parameters are shown in Table IV. Finally, for the largest mass, µ l,R = 26.0 MeV, M is optimal at ∼ 170 MeV (compatible with zero for 80 − 230 MeV).
It appears that the spectral projector method has larger mass dependent discretization effects compared to the gluonic definition, as α/r 0 has larger magnitude, and smaller mass independent discretization effects, as c is smaller. It should be noted, however, that the statistical error in α/r 0 is large enough that the values for both spectral projectors and the gluonic definition are compatible.  Table I. Resulting fit parameters are shown in Table IV.  Table I and II. Fit parameters are shown in Table IV.

B. χ(a → 0) LO Fit
In order to obtain an estimate of a systematic error resulting from the choice of the ansatz for the fit, we have used a second method for extracting the continuum values of the chiral condensate.
While each ensemble used has a different quark mass, each lattice spacing has a small, medium, and large quark mass ensemble. First, we choose one representative small, medium, and large quark mass and find the comparable values of topological susceptibility for each lattice spacing by performing a small linear interpolation to the same mass. For example, for the small quark mass, we choose µ l,R = 12.6 MeV, which is the mass of the A30.32 ensemble, seen in Table I The continuum limit of the topological susceptibility is taken for these three masses and obtained by performing a fit linear in (a/r 0 ) 2 , using the ansatz: where r 4 0 χ(a = 0) is the continuum value of the susceptibility. The square brackets again denote the actual fit parameters. The extracted values of the susceptibility can be found in Table V. The results of this procedure are collected in Figs. 5a, 5b, and 5c, for the small, medium, and large quark masses, respectively. A comparison of the discretization effects, contained in the parameter K, is shown for a range of values of M in Fig. 7.
To calculate the chiral condensate, each extracted continuum value of the topological susceptibility is plotted against the quark mass and fitted to the continuum lowest order χ PT formula Eq. (12), or in dimensionless form: This fit is shown in Fig. 6 and the corresponding values of the renormalized chiral condensate are presented in Table VI.
From Figs. 5a, 5b, and 5c, it is apparent that for a given quark mass all values of M and the gluonic definition converge to the same value in the continuum limit, within errors. This is as expected, as the definitions should only differ at finite lattice spacing.
We observe that the spectral projector definition produces a topological susceptibility with much smaller discretization effects compared to the gluonic definition, regardless of the choice   The fit values of r 0 Σ 1/3 seen in Table VI are larger but within errors of the values found using the global fit (cf. Table IV). The statistical errors are smaller, likely due to being more constrained, only fitting to one parameter. However, the χ 2 per degree of freedom is larger for the spectral projectors results, for which reason we prefer the values from the global fit.

V. CONCLUSIONS
In this paper, we have compared the calculation of the topological susceptibility and the chiral condensate using the spectral projector method to an O(a 2 )-improved gluonic definition using the gradient flow. The topological susceptibility computed using spectral projectors was found to have much smaller discretization effects in the considered setup. It should be noted, however, that in the case of the gradient flow an improved definition could be derived which reduces cut-off effects significantly [36]. It would therefore be interesting, whether such an improved definition can also be found for the here considered topological susceptibility. We have also determined that Where we have used r 0 = 0.474 fm [25] to convert to physical units. Comparing to our previously calculated values, this result is slightly larger. One value, r 0 Σ 1/3 = 0.651(61) [34] is from a fit of the leading order χ PT to the quark mass dependence of the topological susceptibility evaluated with stochastic spectral projectors on a similar set of ensembles. Another, r 0 Σ 1/3 = 0.689 (16)(29) [37] resulted from direct extraction of the condensate from the mode number evaluated also with stochastic spectral projectors on a similar set of ensembles. While slightly larger, this result is still compatible within errors.
By choosing an appropriate value of M , we have shown the method of spectral projectors eliminates discretization artifacts in topological susceptibility up to our computed percision. This is highly relevant for calculating topological charge dependent quantities, when there is only one lattice spacing available and a continuum extrapolation cannot be performed. This is at present the case with the physical point ensembles recently generated by the ETMC [38][39][40]. We believe that this approach will significantly reduce discretization errors due to the topological charge, while preserving reasonable statistical errors, allowing more accurate computation of experimentally relevant physical parameters.

VI. ACKNOWLEDGEMENTS
We thank the European Twisted Mass Collaboration for generating ensembles of gauge field configurations that we have used for this work and for a very enjoyable collaboration. We acknowledge helpful discussions with Dean Howarth on technical details. K.C. was supported in part