Universality of the topological susceptibility in the SU(3) gauge theory

The definition and computation of the topological susceptibility in non-abelian gauge theories is complicated by the presence of non-integrable short-distance singularities. Recently, alternative representations of the susceptibility were discovered, which are singularity-free and do not require renormalization. Such an expression is here studied quantitatively, using the lattice formulation of the SU(3) gauge theory and numerical simulations. The results confirm the expected scaling of the susceptibility with respect to the lattice spacing and they also agree, within errors, with computations of the susceptibility based on the use of a chiral lattice Dirac operator.


Introduction
In QCD and other non-abelian gauge theories, the discussion of the effects of the topological properties of the classical field space tends to be conceptually non-trivial, because the gauge field integrated over in the functional integral is, with probability 1, nowhere continuous. The topological susceptibility, for example, is only formally given by the two-point function of the topological density at zero momentum, unless a prescription is supplied of how exactly the non-integrable short-distance singularity of the two-point function is to be treated.
In lattice gauge theory, the problem was reexamined some time ago [1][2][3] starting from a formulation of lattice QCD which preserves chiral symmetry. An important result of this work was that the topological susceptibility can be written as a ratio of expectation values of other observables which remains well-defined in the continuum limit. A particular choice of regularization is then not required, i.e. the new formula provides a universal definition of the susceptibility. Moreover, this definition is such that the anomalous chiral Ward identities are fully respected.
The aim in the present paper is to complement these theoretical developments by demonstrating the suitability of the universal definition for the computation of the topological susceptibility in lattice gauge theory. In this study, the pure SU(3) gauge theory is considered and a recently proposed version [4] of the universal formula is used (see sect. 2). As far as the feasibility of the calculation is concerned, the results are however expected to be directly relevant for QCD too.

Singularity-free expressions for the topological susceptibility
The formula for the susceptibility obtained in [4] is not very complicated, but some preparation is required to be able to write it down. From the beginning, the theory is considered on a finite hypercubic lattice with spacing a, volume V and periodic boundary conditions. While some particular choices have to be made along the way, these details are expected to be irrelevant in the continuum limit in view of the fact that the expression is renormalized and free of short-distance singularities.

Spectral-projector formula
The construction starts by adding a multiplet of valence quarks with bare mass m 0 to the theory. On the lattice, the added fields are taken to be of the Wilson type [5] and the associated massive Dirac operator D m is assumed to include the Pauli term required for O(a) improvement [6,7] (the relevant improvement and renormalization constants are collected in appendix A).
The hermitian operator D m † D m has a complete set of orthonormal eigenmodes with non-negative eigenvalues α. On average there are only few eigenvalues below some threshold α th proportional to the square of the valence-quark mass (see fig. 1). Above the threshold, the spectrum has an approximately constant density with a slight downward trend in the range considered in the figure. Such a trend is absent in two-flavour QCD [4], but is qualitatively in line with the behaviour of the spectral density at next-to-leading order of quenched chiral perturbation theory [8].
The topological susceptibility is now given by [4]  where P M denotes the orthogonal projector to the subspace spanned by the eigenmodes of D m † D m with eigenvalues α < M 2 . It is taken for granted in this formula that M 2 is above the effective threshold α th of the spectrum and that the renormalized valence-quark mass m R as well as the renormalized value M R of M are held fixed when the lattice spacing is taken to zero (cf. appendix A).

Alternative expressions
Equation (2.1) derives from a study of the renormalization and symmetry properties of the n-point correlation functions of the scalar and pseudo-scalar densities of the valence quarks. There exist various representations of the topological susceptibility of a similar kind, all having the same continuum limit. In particular, for any function R M of D m † D m which is equal to unity in the vicinity of the spectral threshold and rapidly decaying above M 2 . The shape of the function can otherwise be chosen arbitrarily and only affects the size of the O(a 2 ) corrections.
In the numerical work reported later, R M is set to the rational approximation of the projector P M previously used in [4] for the computation of the mode number in For the reader's convenience, the function is given explicitly in appendix B.

Numerical studies
The expression on the right of eq. (2.2) is a ratio of well-defined expectation values that can in principle be computed through numerical simulations. In practice, the traces Tr{. . .} can normally not be evaluated exactly, but as explained in subsect. 3.3, they can be estimated stochastically with a moderate computational effort and without compromising the correctness of the final results.

Simulation parameters
The studies reported in this paper are based on simulations of the lattice theory at three values of the inverse bare gauge coupling β = 6/g 2 0 (see table 1). A well-known deficit of all currently available simulation algorithms for non-abelian gauge theories (including the link-update algorithms used here) is the fact that the integrated autocorrelation times of quantities related to the topological charge are rapidly growing when the lattice spacing decreases [9,10]. In order to guarantee the statistical independence of the N cnfg gauge-field configurations used for the "measurement" of the topological susceptibility, the distance in simulation time of subsequent configurations was required to be at least 10 times larger than the relevant autocorrelation times.
Physical units are defined through the Sommer reference scale r 0 = 0.5 fm [11]. In the range of the gauge coupling covered here, the conversion factor r 0 /a from lattice to physical units was accurately determined by Guagnelli et al. [12]. The spacings of the three lattices thus decrease from roughly 0.1 to 0.05 fm by factors of 1/ √ 2, while the lattice sizes in physical units are approximately constant.

Spectral projector parameters
As already mentioned, the operator R M is taken to be a rational approximation to the projector P M . It thus depends on the valence-quark mass, the mass M and the parameters n and ǫ that determine the accuracy of the approximation (cf. appendix B). A reasonable choice of the latter, previously made ref. [4], is n = 32 and ǫ = 0.01. In the range of eigenvalues of D m † D m below 0.85 × M 2 , the approximation error is then smaller than 2.2 × 10 −4 , which is by far small enough to guarantee the absence of significant systematic effects in eq. (2.2). Moreover, the contribution of the high modes is safely suppressed.
The valence-quark mass and the mass parameter M were adjusted such that their renormalized values in the MS scheme at normalization scale µ = 2 GeV are about 25 and 100 MeV, respectively. Using the information collected in appendix A, the corresponding values of the bare mass parameters, κ = (8 + 2am 0 ) −1 and aM , can be worked out and are listed in table 1. On the lattices considered, there are then 57 − 70 eigenmodes of D m † D m with eigenvalues below M 2 and an average density of roughly 1 such mode per fm 4 .
As already emphasized, the calculated values of the topological susceptibility are not expected to strongly depend on all these details and should in any case always extrapolate to the same value in the continuum limit. The lattice effects are unlikely to be small, however, if aM is not much smaller than 1 or if the expectation values on the right of eq. (2.2) would be dominated by the modes up to and slightly above the spectral threshold, where the effects are kinematically enhanced. Both of these unfavourable situations are avoided by the above choice of the mass parameters.

Random-field representation
In lattice QCD, random field representations were introduced many years ago [13] and are now widely used. The application of the method in the present context requires a set η 1 , . . . , η N of N pseudo-fermion fields to be added to the theory with action where the bracket (η, ζ) denotes the obvious scalar product of such fields. For every gauge field configuration, these fields are generated randomly so that one obtains a representative ensemble of fields for the complete theory. In the rest of this section, expectation values are always taken with respect to both the gauge field and pseudofermion fields. The stochastic observables may now be introduced and a moment of thought then shows that the expectation values on the right of eq. (2.2) are given by The topological susceptibility can thus be computed by calculating the expectation values of A, B and C 2 . For a given gauge-field configuration, the evaluation of these observables requires the fields R M η k , R 2 M η k and R M γ 5 R M η k to be computed, i.e. the total numerical effort per configuration is roughly equivalent the one required for 3N applications of the operator R M to a given pseudo-fermion field.
From this point of view, small values of N are favoured, but a good choice of N must also take into account the fact that the variance of the stochastic observables decreases with N . Some experimenting then shows that setting N = 6 is a reasonable compromise at the specified values of the mass parameters. Since R M is a rational function of D m † D m of degree [2n + 1, 2n + 1], the measurement of the stochastic observables requires the (twisted-mass) Dirac equation to be solved for altogether 2340 source fields. The computational load thus tends to be heavy, but the problem is well suited for the application of highly efficient solver techniques such as local  (26) deflation [14] (see ref. [15] for a recent review of the subject). In particular, when these are used, the effort scales linearly with the lattice size and is nearly independent of the values of the mass parameters.

Simulation results
The simulation data discussed in the following paragraphs are summarized in table 2. In all cases, the statistical errors were estimated using the jackknife method and were combined in quadrature with the quoted scale errors (where appropriate).
(a) Mode number. The average number ν of eigenmodes of D m † D m with eigenvalues below M 2 is an extensive quantity and is therefore normalized by the lattice volume in table 2. At the specified bare masses, the renormalized masses m R and M R are practically equal to 25 and 100 MeV, respectively, on all three lattices considered. In view of the renormalization properties of the mode number [4], the calculated values of ν/V are thus expected to be the same up to O(a 2 ) effects.
Within errors, the figures listed in the second column of table 2 in fact coincide with one another. Note that the quoted errors do not take into account the fact that the mass renormalization factors and thus the renormalized values of the masses are only known up to an error of about 2% (see appendix A). Once this error is included in the analysis, one can still conclude, however, that the simulation results confirm the expected scaling of the mode number to the continuum limit at a level of precision of 3% or so.
(b) Topological susceptibility. The values of the susceptibility calculated along the lines of the present paper are listed in the third column of table 2. Again one observes no statistically significant dependence on the lattice spacing. Finite-volume effects are, incidentally, known to be negligible with respect to the statistical errors on all three lattices [16].
Fits of the data by a constant or a linear function in a 2 yield consistent results in the continuum limit. Since the slope in a 2 turns out to vanish within errors, the (more accurate) number (c) Charge sectors & the Wilson flow. An understanding of how exactly the topological charge sectors emerge in the continuum limit has recently been achieved using the Wilson flow [18]. The definition of the topological susceptibility suggested by the sector division is geometrically appealing and computationally far less demanding than the spectral-projector formula (2.2). Presumably the two definitions agree in the continuum limit, but there is currently no solid theoretical argument that would show this to be the case. The values of the susceptibility computed using the Wilson flow are listed in the fourth column of table 2 (see ref. [18] for the details of the calculation). While they appear to be systematically lower than the ones obtained using the spectral-projector formula, the differences are statistically insignificant on each lattice. Moreover, there could be lattice effects of size up to the level of the statistical errors.
Since the same ensemble of representative gauge-field configurations was used in the two cases, the quoted errors are correlated to some extent (not completely so in view of the use of random fields). The ratio listed in the last column of table 2 is therefore obtained with slightly better precision than the susceptibilities. Fits of the ratio by a constant and linear function in a 2 are both possible, the values in the continuum limit being 1.048(14) and 1.036(31), respectively. The spectral-projector and the Wilson-flow definition of the susceptibility thus coincide to a precision of a few percent. While there is some tension in the data, there is no clear evidence for the definitions to be different at this level of accuracy.

Conclusions
The numerical studies reported in this paper confirm the universality of the spectralprojector formula (2.2) for the topological susceptibility. In particular, no statistically significant lattice-spacing effects were observed and the calculated values agree † In ref. [16], a different convention for the conversion from lattice to physical units was used and the value for the susceptibility quoted there is therefore slightly different from the one printed here. with the result obtained by Del Debbio, Giusti and Pica [16], where a chiral lattice Dirac operator was used.
The numerical effort required for the stochastic evaluation of the spectral-projector formula increases proportionally to the number V /a 4 of lattice points, but is a flat function of all other parameters. On large lattices, computations of the susceptibility along these lines thus tend to be more feasible than those based on a chiral lattice Dirac operator (which scale roughly like V 2 ). Even less computer time is required if the susceptibility is defined via the Wilson flow, but a formal proof for this definition to be in the same universality class as the spectral-projector formula is still missing.
With respect to the pure gauge theory, the application of the spectral-projector formula in QCD is not expected to run into additional difficulties. Accurate calculations of the topological susceptibility however require representative ensembles of, say, a few hundred statistically independent gauge-field configurations to be generated. This part of the calculation usually consumes most of the computer time and may rapidly become prohibitively expensive at small lattice spacings [10]. At present, computations of the susceptibility on lattices similar to the ones considered here are therefore not easily extended to QCD with light sea quarks.
We wish to thank Leonardo Giusti for helpful discussions on various issues related to this work. All numerical calculations were performed on a dedicated PC cluster at CERN. We are grateful to the CERN management for providing the required funds and to the CERN IT Department for technical support. F. P. acknowledges financial support by an EIF Marie Curie fellowship of the European Community's Seventh Framework Programme under contract number PIEF-GA-2009-235904.

A.1 Dirac operator and renormalization constants
The lattice theory considered in this paper is set up as usual, using the Wilson gauge action and the standard O(a)-improved Wilson-Dirac operator. The notation and normalizations are as in ref. [7]. In particular, c sw and c A denote the coefficients of the Pauli term in the Dirac operator and the O(a) term required for the improvement of the axial current. Here these coefficients were set to the values given by the nonperturbatively determined interpolation formula quoted in ref. [19] (see table 3).
The values of the renormalization constant Z A of the axial current listed in table 3 were obtained by evaluating the interpolation formula given in ref. [20]. In the case  (13) of the renormalization constant Z P of the pseudo-scalar quark density, the quoted values are the ones required to pass from the lattice normalization of the density to the one in the MS scheme of dimensional regularization at normalization scale µ = 2 GeV. The constant was calculated in two steps, first passing from the lattice to the renormalization-group-invariant normalization [21] and then from there to the MS scheme [22].

A.2 Quark masses
The renormalized quark mass in the MS scheme is given by where m is the bare current-quark mass, m q = m 0 − m c the subtracted bare mass and m c the critical bare mass. Here and below, b X (where X = A, P, . . .) denotes an improvement coefficient required to cancel lattice effects proportional to am q . At fixed gauge coupling, the current quark mass is related to the subtracted bare mass through Using the Schrödinger functional, the coefficients Z and b m −b A +b P were determined non-perturbatively by Guagnelli et al. [23] (see table 4). Also shown in table 4 is the current quark mass at one value of the hopping parameter κ = (8 + 2am 0 ) −1 . Together with eq. (A.2), these data allow the current quark mass to be estimated at larger values of κ, where a direct computation on the lattices considered in this paper tends to be compromised by the presence of accidental near-zero modes of the Dirac operator.

A.3 Renormalization of the mode number
The renormalization and improvement properties of the mode number were discussed in detail in ref. [4]. In particular, it was shown there that is only known to one-loop order of perturbation theory [24,4]. for some specified (small) value of ǫ. This choice ensures that h(x) provides a uniform approximation to the step function in the range |x| ≥ √ ǫ, with maximal absolute deviation equal to 1 2 δ. Moreover, inspection shows that h(x) decreases monotonically in the transition region − √ ǫ ≤ x ≤ √ ǫ.
For a given degree n and transition range ǫ, the coefficients of the minmax polynomial P (y) can be computed numerically using standard techniques. An efficient procedure was described in ref. [25], for example. The mass M * ∝ M is then determined through As explained in appendix B of ref. [4], this convention is intended to minimize the deviation Tr{P M − R 4 M } . In the present context, other choices of M * would however do just as well, since eq. (2.2) is expected to hold for any M .
Small approximation errors δ are achieved with moderately high degrees n if ǫ is not very small. For n = 32 and ǫ = 0.01, for example, one obtains