Topology via Spectral Projectors with Staggered Fermions

The spectral projectors method is a way to obtain a theoretically well posed definition of the topological susceptibility on the lattice. Up to now this method has been defined and applied only to Wilson fermions. The goal of this work is to extend the method to staggered fermions, giving a definition for the staggered topological susceptibility and testing it in the pure $SU(3)$ gauge theory. Besides, we also generalize the method to higher-order cumulants of the topological charge distribution.


INTRODUCTION
The properties related to the existence of field configurations with non-trivial topology, and the associated non-trivial dependence on the topological parameter θ, represent some of the most significant non-perturbative aspects of QCD and QCD-like theories. Monte Carlo simulations on a lattice are the most natural first principle tool to investigate such properties, however the fact that on a discrete space-time homotopy classes are not well defined makes the issue non-trivial. In principle, many definitions of topological charge can be assigned in the discretized theory, all consistent with each other in the continuum limit, however discretization errors can be different depending on the choice.
The topological charge in continuum Yang-Mills theories is defined in terms of gluon fields as follows: where q(x) is the topological charge density, and is integer valued when proper boundary conditions are taken (e.g., periodic boundary conditions for a finite box, or vanishing action density at infinity). The index theorem [1] then relates Q to fermion field properties, in particular where n ± are, respectively, the number of left-handed and right-handed zero-modes of the Dirac operator / D. The possible lattice discretizations can be divided essentially into two different classes: gluonic and fermionic.
Gluonic definitions are based on a straightforward rewriting of Eq. (1) in terms of lattice gauge links. Despite having the correct naïve continuum limit, they are non-integer valued and subject to renormalizations induced by ultraviolet (UV) fluctuations. In particular, correlation functions of Q L must be renormalized both additively [2] and multiplicatively [3] in order to match the corresponding continuum quantities apart from finite O(a) corrections (where a is the lattice spacing): this is * Electronic address: claudio.bonanno@pi.infn.it † Electronic address: giuseppe.clemente@pi.infn.it ‡ Electronic address: massimo.delia@unipi.it § Electronic address: sanfilippo@roma3.infn.it the case for the topological susceptibility, χ ≡ Q 2 /V , as well as for the higher order cumulants of Q [4] which enter the coefficients of the Taylor expansion of the free energy density in θ. Alternatively, one can make use of smoothing methods which dampen UV fluctuations of gauge fields while leaving the global topological background unchanged, thus leading to an approximately integer valued topological charge: various similar methods have been proposed, such as cooling [5] or the gradient flow [6,7], all leading to equivalent results [8,9].
Fermionic definitions, being based on a counting of zero modes, are in principle better founded, and would be realized in practice by a simple evaluation of the trace of the γ 5 operator on a basis of eigenvectors of the Dirac operator. However, also in this case one has to face problems related to the difficulty in implementing fermions with the correct chiral properties on a lattice. The best approximation is provided by discretizations of the Dirac operator satisfying the Ginsparg-Wilson relation [10], an example being the overlap operator [11], which satisfies an approximate chiral symmetry [12], leads to a counting of exact zero modes and has been widely used as a tool to extract the topological content of gauge configurations [13,14].
Alternatively, one can use a more standard fermion discretization, either Wilson or staggered based. However, in this case zero modes are non-exact, or do not have a well defined chirality, or both. The counting is not well defined and, if one tries to define Q L as the trace of the discretized γ 5 operator, one needs again to take into account proper renormalizations [15][16][17][18]. This, however, is not an obstruction, and the method based on spectral projectors relies exactly on this strategy. Indeed, in Ref. [19,20] spectral projectors have been used to obtain a theoretically well-posed definition of the continuum topological susceptibility, which is also easily adaptable for numerical simulations on the lattice. In particular, the method has been derived for Wilson fermions and successfully tested both in pure Yang-Mills [21,22] and in full QCD [23].
All the methods exposed above, either fermionic or gluonic, are theoretically well founded or designed so as to match the correct definition of homotopy classes when these become well defined: as a matter of fact, all methods provide consistent results when the continuum limit is taken. The question about which of the methods one should adopt can then be answered based on two considerations: numerical convenience, i.e. the computational effort required by a given definition of Q L , and the magnitude of residual corrections to the continuum limit.
The second issue can be particularly relevant for numerical simulations involving light dynamical quarks. Actually, in this case most of the lattice artifacts stem from the discretization of the fermion determinant: the presence of zero modes should suppress configurations with non-zero topological charge, eventually leading to the absence of θ-dependence in the limit of massless quarks; however, a lattice discretization with non-exact chiral properties typically leads to a less efficient suppression, because of large would-be zero modes, thus leading to somewhat larger values of the topological susceptibility at finite lattice spacing. This problem can make the approach to the continuum limit particularly difficult, both at zero and finite temperature [24][25][26][27][28], making it necessary to perform simulations at lattice spacings much smaller than those usually adopted in quenched simulations. A possible heuristic solution adopted in the recent literature has been to reweight gauge configurations by hand, according to the lowest eigenvalues of the dynamical fermion operator [26] Even if the above problem is related to the discretized path integral measure, rather than to the discretized observable, it is not inconceivable that the choice of a proper fermion discretization for the topological charge could ameliorate the convergence to the continuum, especially if the discretization matches the one adopted in the measure. This possibility is actually supported by a recent study [23], investigating full QCD with twisted mass Wilson fermions, in which strongly reduced lattice artifacts are observed for the zero-temperature topological susceptibility if a definition based on twisted mass spectral projectors is adopted, instead of other standard gluonic definitions.
We can now come to the main point of our study: we would like to extend the definition of topological quantities based on spectral projectors to the case of staggered fermions, the main motivation being to adopt it in ongoing lattice investigations of θ-dependence in full QCD with staggered fermions [28]. Most of our discussion will focus on how to properly define and renormalize the definition of the topological susceptibility based on spectral projectors in the case of staggered fermions; we will also present some numerical results which however, given the goals of this paper, will be limited to measurements taken on quenched ensembles; the case of full QCD ensembles will be treated separately in an upcoming work. In addition, we will also show how the spectral projectors method can be exploited to define and evaluate cumulants of the topological charge higher than just the topological susceptibility.
The paper is organized as follows. In Section 2, after a brief review of the method for Wilson fermions, we extend the spectral projectors method to the case of staggered fermions, deriving spectral expressions for the topological susceptibility and for all higher-order cumulants; moreover, in the same section, we also describe the numerical strategies we adopted to test spectral projectors on the lattice. In Section 3 we present numerical results for the pure SU (3) gauge theory and finally, in Section 4, we draw our conclusions and discuss future perspectives. In this Section, after a brief review of the main ideas underlying the method of spectral projectors for Wilson fermions, we show how they can be extended to the staggered case, obtaining a similar expression for the topological susceptibility. We also discuss about the extension of the method to higher order cumulants and about a practical way to fix the cut-off scale adopted in the method.
A. Topological susceptibility via spectral projectors: the Wilson case As for other definitions of topological charge based on the index theorem, the starting point is to write it in terms of the trace of the γ 5 operator, Q 0 = Tr{γ 5 }. When the trace is taken over eigenvectors of the lattice Wilson fermion operator, this definition is subject to a multiplicative renormalization, because chiral symmetry is explicitly broken by Wilson fermions. In particular, making use of non-singlet chiral Ward identities (see Ref. [15] for more details), one can show that the renormalized charge can be expressed as [19,20] where Z (ns) S and Z (ns) P are, respectively, the renormalization constants of the non-singlet scalar and pseudo-scalar fermionic densities The correct renormalization factor of the charge can be easily obtained from its bare expression once the renormalization constants of the densities are chosen according to the non-singlet Ward identities written for Wilson fermions (for further details about this topic, we refer to Ref. [15,16]). Besides, note that the ratio Z (ns) S /Z (ns) P is different from 1 at finite lattice spacing because the Wilson operator D W explicitly breaks chiral symmetry [16].
The renormalization constants Z can be obtained by a variety of non-perturbative means. In the spectral projectors approach, one writes directly their ratio in terms of the so-called "spectral sums": where since spectral sums can be expressed in terms of density chains, whose renormalization properties are known [20] Note that spectral sums (7) and (8) lead to pseudo-scalar density chains (9) and (10) because of the γ 5 -hermiticity of the Wilson operator: γ 5 D † W γ 5 = D W . Eq. (6) holds for high enough values of k and l, but also if the inverse powers of D † W D W are traded for a generic, fast-decreasing function f (D † W D W ) (see Ref. [20] for more details). Choosing the Heaviside function f (x) = θ(M 2 − x), one can easily evaluate the traces in the spectral basis and obtain where P M is the orthogonal projector on eigenspaces of the Wilson operator with eigenvalues |λ| ≤ M Analogous considerations can be applied to the fermionic definition of the bare charge Q 0 , which can be rewritten in terms of a spectral expression as well: in the continuum, it would suffice to project just on the kernel of the operator, however this is of course not true at finite lattice spacing and for a fermion operator with non-exact zero modes. Finally, one can write the renormalized definition of the lattice topological susceptibility via spectral projectors as follows: This definition presents only O(a) or O(a 2 ) corrections (depending on the explicit discretization) if the cut-off mass M is properly tuned (i.e. its renormalized value is kept fixed) as the continuum limit is taken [19][20][21]. An important point, which is worth stressing, is that, because of the fast decreasing behavior in the ultraviolet of the functions appearing in the spectral sums (and in particular of the projector P M ), the above expressions are free of short distance singularities, so that no further renormalizations are needed apart from the multiplicative one. In particular, no additive renormalization appears for the topological susceptibility, contrary to what happens for the standard gluonic definition because of contact terms. In this sense, the spectral projectors method has some analogies with filtering methods [29][30][31], in which a projection onto the eigenspace of the lowest eigenvectors of the Dirac operator is used as a smoothing technique.
B. Topological susceptibility via spectral projectors: the staggered case A bare version of the index theorem can be written for the staggered Dirac operator D st by just taking into account that, in the continuum, it describes 2 d/2 degenerate flavors of dynamical fermions, where d is the space-time dimension, so that the number of zero modes corresponds to 2 d/2 times the topological charge. Therefore, one can start from the bare definition Q 0st = (−2) −d/2 Tr{Γ 5 }, where Γ 5 is the staggered version of γ 5 (more precisely Γ 5 → γ 5 ⊗ Id in the continuum limit, see, e.g., Ref. [17] for an explicit expression). Chiral symmetry is partially broken also for staggered fermions, so that Q 0st renormalizes multiplicatively as in the Wilson case, however the renormalization constants are different, since the breaking pattern is not the same. Indeed, the staggered lattice action is invariant under a remnant of the chiral symmetry, where Γ 55 = γ 5 ⊗ γ 5 . We refer the reader to Refs. [17,18] for a detailed discussion of the anomalous Ward identities for staggered fermions. The final result for the renormalized staggered charge, which has been obtained writing the Witten-Veneziano equation starting from the renormalized singlet axial Ward identity, is the following where the constants Z (s) S and Z (s) P appear in the inverse order with respect to Wilson fermions, and refer respectively to the scalar and pseudo-scalar flavor-singlet (compared to non-singlet in the Wilson case) bare densities, S 0 and P 0 : meaning that the corresponding renormalized quantities read 1 In the staggered case the ratio Z and 1 Note that our notation differs from the one employed in [17,18]. The bare singlet scalar and pseudo-scalar densities are written as S 1 and P 1 , cf. Eqs. (2.5) and (2.6) in [17], and their renormalization constants are written in terms of the quark mass one m R = Z −1 m m as Z  [17] and also in [18]. through the ratio We note that also in this case the ratio of spectral sums is inverted with respect to Eq. (6) for Wilson fermions. This is due to the fact that we are dealing with singlet, rather than non-singlet, densities.
Following the same line of reasoning of Ref. [20], one can trade also in this case the inverse powers of D † st D st for a fast-decreasing function, in particular for P M , where now P M is the orthogonal projector on eigenspaces of the Dirac operator with purely imaginary eigenvalues −iλ such that λ 2 ≤ M 2 . That leads to the following expression for the ratio and finally, using also in this case a spectral projector definition for the bare topological charge we obtain the following expression for the topological susceptibility of staggered fermions in d dimensions: which coincides with Eq. (14) apart from the factor 2 −d , related to the taste degeneration, as previously explained. Also in this case, the same considerations apply regarding the fast decreasing behavior of the projector in the ultraviolet, leading to the absence of short distance singularities and additive renormalizations. As for Wilson fermions, the cut-off mass M appears as a free parameter of the definition. If the zero modes were exact, one could take M arbitrarily small. Exact zero modes are obtained for the overlap version of the staggered operator [32,33], however, for the standard staggered operator D st they are shifted by lattice artifacts because of the explicit breaking of the chiral symmetry [18], so that one must extend the sum up to a certain cut-off eigenvalue M , keeping its renormalized value fixed as the continuum limit is approached (see Subsection 2 D for more details).

C. Higher-order terms of the θ-expansion via spectral projectors
The θ-dependence of the vacuum energy (free energy) density can be parametrized around θ = 0 as follows [34]: where the b 2n coefficients, which parametrize the corrections to the quadratic behavior of f (θ), are defined as: where Q k c denotes the k th -order cumulant of the probability distribution P (Q).
These quantities can be computed in terms of spectral projectors as well, exploiting in particular the fact that, because of the absence of short distance singularities, only multiplicative renormalizations have to be taken into account. In particular, following the same line of thoughts of the topological susceptibility, it is easy to prove the following general expression The above expression has been written explicitly for the case of staggered fermions but, of course, the final expression holds for Wilson fermions as well, after omitting the factor 2 −dn , again related to taste degeneracy.

D. Numerical implementation and remarks on the choice of the cut-off mass M
In the case of Wilson fermions, different strategies have been adopted in the literature for the evaluation of the traces appearing in Eq. (14), either by means of noisy estimators [21,22], or by an explicit computation, configuration by configuration, of all relevant eigenvectors of the Dirac operator entering the traces [23]. In our numerical implementation we have followed the second strategy, i.e. we evaluated the traces appearing in Eqs. (24) and (27) expressed the projector P M through the eigenvectors of D st : limiting the sum over eigenvalues up the threshold λ max = aM . The relevant quantities entering the expressions for the topological susceptibility and the higher order cumulants are then where ν(M ) is the number of eigenvalues with |λ| ≤ aM .
As already explained in Subsection 2 B, the choice of the cut-off mass M is irrelevant in the continuum limit, since index theorem states that only zero-modes contribute to topology. However, corrections to the continuum limit do depend on it, and it is possible to show that lattice artifacts are O(a 2 ) if the renormalized value of the cut-off mass, M R , is kept fixed as the lattice spacing is varied [17,20]): where χ is the continuum value of the topological susceptibility. Therefore, one has to tune M as a function of a in order to keep M R fixed. Most of the following discussion applies specifically to the case of staggered fermions, for which the mass renormalizes as follows [17] where Z (s) S is the renormalization constant of the singlet scalar density mentioned above. This quantity is not accessible separately in terms of spectral sums, however in many lattice studies it is already known by other means. For instance, in numerical simulations of full QCD performed on a line of constant physics, one already tunes the bare quark masses as a function of the lattice spacing so as to keep the physics, hence the renormalized quark masses, unchanged: in these cases it will suffice to keep the ratio of M to any of the bare quark masses unchanged as the continuum limit is approached.
However, in the present study we consider as a numerical test-bed the case of the pure gauge theory at zero temperature, for which the strategy above cannot be applied. In this case, trying to avoid a direct computation of Z (s) S for each lattice spacing, we have devised the following strategy. The number of eigenmodes which are found below a given threshold scales proportionally to the total lattice volume, i.e. the density of eigenmodes with |λ| < M , ν(M )/V , is a constant as the thermodynamical limit V → ∞ is approached. Moreover, if the renormalized threshold M R is kept fixed, the density of eigenmodes is expected to be independent of the lattice spacing. This is supported by leading order chiral perturbation theory and the Banks-Casher relation: in the large-volume and chiral limit, and for small enough M , one has [20] where Σ is minus the chiral condensate in the thermodynamic and chiral limit, and we have used the fact that Z (s) S is the renormalization constant for both the singlet scalar density and the inverse mass, so that ΣM is a renormalization group invariant quantity. Therefore, our prescription in the following will be to keep the bare quantity ν(M ) /V fixed (with the spacetime volume expressed in physical units) in order to maintain M R constant as the lattice spacing is changed.

NUMERICAL TESTS IN THE PURE GAUGE THEORY
In order to test the definition of topological quantities via staggered spectral projectors, we have considered the pure SU (3) Yang-Mills theory. Configurations have been generated using the standard Wilson plaquette action for N c = 3: and a standard local algorithm consisting of a 4:1 mixture of over-relaxation and over-heatbath. For simulations at zero temperature, we have considered 4 different lattice spacings, corresponding to β = {5.9, 6.0, 6.125, 6.25}, and symmetric lattices L × L × L × L with L in the range 1.2 -1.8 fm (see Table I), which for the pure gauge theory is large enough to ensure the absence of significant finite size effects, in particular for topological quantities. In the following we will express physical quantities, as well as the lattice spacing, in terms of the Sommer parameter r 0 ≃ 0.5 fm.
In all cases, we have collected 300 well decorrelated configurations, on which topological quantities have been measured both by staggered spectral projectors and, in order to make a comparison, with a standard gluonic definition of the topological charge. We note that statistics are not large, because the main purpose of our numerical simulations is to test the staggered definition of spectral projectors and not to make a precision study about topology in the pure gauge theory.
Concerning the gluonic definition, in this work we adopted the clover discretization: As for the smoothing method, we decided to apply the standard cooling procedure, performing 80 cooling sweeps for each configuration (the topological susceptibility was stable already after 30 sweeps). The cooled topological charge has been further rounded to the closest integer following the procedure described in Refs. [35,36] and then used to compute the gluonic definition of topological susceptibility via A. Spectral determination of χ at T = 0 To start with, in Figure 1 we show the topological susceptibility χ SP obtained via spectral projectors for β = 6.25 and different values of the bare cut-off mass M , comparing it with the gluonic determination on the same sample of configurations. We observe an approximate plateau in a wide range of M , where χ SP is in good agreement with the gluonic definition. Such a plateau is not required a priori, however it is reasonable to expect it: the cut-off mass M filters away fluctuations at the UV scale, in particular M −1 can be viewed as the analogous of the smoothing radius for smoothing techniques, so that the appearance of the plateau is the signal of a well defined separation between the UV scale and the physical scale of topological excitations.
In order to extrapolate χ SP towards the continuum, we considered determinations at fixed values of the renormalized mass M R . To keep M R fixed, we measured the dependence of ν(M ) /V on M so that, once fixed a particular value of the mode density, we could find, for each β, the value of M corresponding to the same renormalized mass M R . Fig. 2 illustrates this procedure for two different values of M R employed for the continuum extrapolation. In Table I we report, for each β, the corresponding value of the lattice spacing, the volume in lattice units and the measures of the topological susceptibility  obtained with spectral projectors and with the gluonic definition.
As shown in Table II, the continuum value of χ obtained with spectral projectors is independent of the choice of M R and well compatible with the gluonic determination within the errors. For the sake of completeness, we also report determinations of χ obtained by other fermionic methods, in particular using the overlap operator and Wilson spectral projectors. They all agree, within errors, with our staggered spectral determination. Fig. 3 shows the extrapolation towards the continuum both for χ SP and χ gluo . Lattice artifacts have a slight dependence on the cut-off mass M R , which is however well contained within errors and comparable in magnitude to that affecting the gluonic definition. This is similar to what happens  Table I: Determinations of χSP and χ gluo for each β. The lattice spacings in units of the Sommer parameter r0 were taken from Ref. [37]. The values of the renormalized cut-off masses M1 and M2 correspond, respectively, to r 4 0 ν /V = 1 · 10 −3 and 3 · 10 −3 . Assuming Eq. (34) and the values of r0 and ΣR measured respectively in [38] and [39], their values in physical units are M1 ≃ 33 MeV and M2 ≃ 98 MeV.
in the case of Wilson spectral projectors [21,22].    Table I. B. Spectral determination of b2 at high T The next-to-leading coefficient in the θ-expansion of the free energy energy is (see Eqs. (25) and (26)) The gluonic discretization simply yields: Instead, from Eq. (27) we get the staggered spectral expression: The measure of b 2 at zero temperature requires in general quite large statistics, because it is necessary to detect deviations from gaussianity of the topological charge distribution P (Q), which are small [4,35,36,[40][41][42] and become less and less visible as the lattice volume is increased. For this reason, we decided to test the numerical determination of b 2 via spectral projectors in the hightemperature, deconfined phase of the SU (3) pure gauge theory, since in that regime its value is larger than in the T = 0 case, approaching the prediction b 2 = −1/12 from the Dilute Instanton Gas Approximation (DIGA), while at the same time the width of the distribution (proportional to the topological susceptibility) is smaller [43]. We have considered, in particular, a determination at β = 6.305 on a 30 3 × 10 lattice, corresponding to a temperature T ≃ 338 MeV ≃ 1.145 T c , for which a determination of b 2 by the gluonic method has been already reported in Ref. [43]. For a finite temperature implementation of the spectral projectors method, an ambiguity could emerge as to whether the first ratio in Eq. (40), corresponding to the multiplicative renormalization (Z (s) P /Z (s) S ) 2 , should be computed in the finite T simulation or instead at zero T . In principle, renormalization constants should be independent of infrared (IR) conditions such as the temperature scale. In order to check that, we have computed the ratio both from the finite temperature simulation and from a dedicated simulation on a symmetric 30 4 lattice at the same value of β: results are shown and compared in Fig. 4, where it clearly appears that, apart from the lowest values of M , for which the sensibility to IR conditions is large, the two determinations are in reasonable agreement with each other.
Finally, in Fig. 5, we show results obtained for b SP 2 , and using the zero temperature renormalization constants, as a function of the bare cut-off mass M . In this case we do not fix a particular value of M , since we have data at a single value of β, hence we do not aim at performing the continuum extrapolation; however we notice that results are in good agreement, over a wide range M , with the gluonic determination of b 2 performed on the same configuration sample, as well as with the determination of Ref. [43] with the same β and lattice size (−12 b 2 = 1.10 (7)).

DISCUSSION AND CONCLUSIONS
In this work we have defined the extension of the spectral projectors method to the case of staggered fermions. Despite the different patterns of chiral symmetry breaking, the final formula for the topological susceptibility, Eq. (24), turns out to be practically identical to the one for Wilson fermions, when the proper staggered discretization of the γ 5 operator and the fourfold degeneracy of staggered fermions are taken into account. Moreover, the method has been extended to all higher-order cumulants of the topological charge distribution, which enter the Taylor expansion in θ of the free energy density. The method has then been tested in the pure SU (3) gauge theory, both at zero temperature and, for the fourth order cumulant, at finite T , with results in agreement with previous results in the literature obtained by other fermionic or gluonic definitions of the topological charge.
Corrections to the continuum limit turn out to be of the same order of magnitude as those observed for the gluonic