Comparison of topological charge definitions in Lattice QCD

In this paper, we show a comparison of different definitions of the topological charge on the lattice. We concentrate on one small-volume ensemble with 2 flavours of dynamical, maximally twisted mass fermions and use three more ensembles to analyze the approach to the continuum limit. We investigate several fermionic and gluonic definitions. The former include the index of the overlap Dirac operator, the spectral flow of the Wilson–Dirac operator and the spectral projectors. For the latter, we take into account different discretizations of the topological charge operator and various smoothing schemes to filter out ultraviolet fluctuations: the gradient flow, stout smearing, APE smearing, HYP smearing and cooling. We show that it is possible to perturbatively match different smoothing schemes and provide a well-defined smoothing scale. We relate the smoothing parameters for cooling, stout and APE smearing to the gradient flow time τ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau $$\end{document}. In the case of hypercubic smearing the matching is performed numerically. We investigate which conditions have to be met to obtain a valid definition of the topological charge and susceptibility and we argue that all valid definitions are highly correlated and allow good control over topology on the lattice.


Introduction
QCD gauge fields can have non-trivial topological properties, manifested in non-zero and integer values of the socalled topological charge. Such topological properties are believed to have important phenomenological implications, a e-mail: andreas.athinodorou@pi.infn.it b e-mail: kcichy@amu.edu.pl (corresponding author) e.g. the fluctuations of the topological charge are related to the mass of the flavour-singlet pseudoscalar η meson [1,2]. The topology of QCD gauge fields is a fully non-perturbative issue, hence lattice methods are well-suited to investigate it. Naively, lattice gauge fields are topologically trivial, since they can always be continuously deformed to the unit gauge field. However, it can be shown that for smooth enough gauge fields (sufficiently close to the continuum limit), the notion of topology of lattice gauge fields is still meaningful [3]. Historically, the first investigation aimed at studying the topological properties of a non-Abelian gauge theory was reported in Refs. [4,5] for the SU(2) gauge group case and then extended to SU(3) in Ref. [6].
Over the years, many definitions of the topological charge of a lattice gauge field were proposed [7,8]. It is clear that the definitions differ in terms of their computational cost and convenience, but also theoretical appeal. These definitions can be characterized either as fermionic or gluonic. During the last decade, a number of efforts have revealed important aspects of the topological susceptibility, which reflects the fluctuations of the topological charge. Universality of the topological susceptibility in fermionic definitions has been demonstrated [9], giving an insight in the basic properties a definition has to obey so that the topological susceptibility is free of short-distance singularities, which have plagued some of the earlier attempts at a proper and computationally affordable fermionic definition [10,11]. On the other hand, considering the gluonic definition, a theoretically clean understanding on how the topological sectors emerge in the continuum limit has been attained through the gradient flow [12][13][14][15][16]. Further projects have also divulged the numerical equivalence of the gradient flow with the smoothing tech-nique of cooling at finite lattice spacing. Previous investigations, on the other hand, have shown numerically that the field theoretic topological susceptibility extracted with several smoothing techniques such as cooling, APE and HYP smearing give the same continuum limit. Although no solid theoretical argument suggests so, it is believed that all different definitions of the topological charge agree in the continuum limit. Good agreement has also been found between the gluonic definiton and a fermionic one in finite-temperature studies [17]. For an overview of topology-related issues on the lattice, we refer to the review paper by Müller-Preussker [8].
This aim of the paper is two-fold. The first purpose is to investigate the perturbative equivalence between the smoothing schemes of the gradient flow, cooling as well as APE, stout and HYP smearing. In Refs. [18,19] it was demonstrated that gradient flow and cooling are equivalent if the gradient flow time τ and the number of cooling steps n c are appropriately matched. By expanding the gauge links perturbatively in the lattice spacing a, at subleading order, the two methods become equivalent if one sets τ = n c /(3 − 15b 1 ) where b 1 is the Symanzik coefficient multiplying the rectangular term of the smoothing action. It is, thus, interesting to extend the study of Ref. [18] and explore whether similar matching can also be derived for APE, stout and HYP smearings. Hence, we carry out such investigation and support it with the appropriate numerical results.
The second motivation of this paper is to attempt a systematic investigation of different topological charge definitions. We have computed them on selected ensembles generated by the European Twisted Mass Collaboration (ETMC) 1 . The included definitions are: 2 -index of the overlap Dirac operator on HYP-smeared and non-HYP-smeared configurations, -Wilson-Dirac operator spectral flow (SF), -spectral projector definition, -field theoretic (gluonic) definition, with gauge fields smoothed using: gradient flow (GF) with different smoothing actions and at different flow times. Namely, we smooth the gauge fields using the Wilson plaquette, Symanzik tree-level and Iwasaki actions at flow times t 0 , 2t 0 and 3t 0 ; t 0 is defined in Sect. 3.5.1. cooling (cool) with the three different smoothing actions and cooling steps matched to GF time for t 0 , 2t 0 and 3t 0 , 1 In 2019, the name of the collaboration has been changed to Extended Twisted Mass Collaboration. 2 For references on each of the following definitions, we refer to Sect. 3.4.
stout smearing with three different values of the stout parameter ρ st and smearing steps matched to GF time at flow times t 0 , 2t 0 and 3t 0 , -APE smearing with three different values of the parameter α APE and smearing steps matched to GF time for t 0 , 2t 0 and 3t 0 , -HYP smearing for a given set of parameters α HYP1 , α HYP2 , α HYP3 , and smearing steps numerically matched to GF time at flow times t 0 , 2t 0 and 3t 0 .
The outline of the paper is as follows: Sect. 2 describes our lattice setup as well as the relevant details regarding the production of the N f = 2 configurations. Section 3 introduces the topological charge definitions that we are using and includes the derivation of matching conditions between different smoothing schemes. In Sect. 4, we discuss and compare different definitions of the topological charge, we analyze the approach to the continuum limit and we show results for the topological susceptibility. Finally, in Sect. 5 we summarize and conclude.

Lattice setup
Motivated by the desire to cover as many definitions of the topological charge as possible, including the costly overlap definition, we performed our comparison of different topological charge definitions on small volume ensembles (to keep the computational cost affordable and at the same time incorporate all definitions) generated by ETMC with N f = 2 [20][21][22] dynamical twisted mass fermions. The action in the gauge sector is Re Tr 1 − P 1×1 x;μ,ν +b 1 4 μ,ν=1 μ =ν Re Tr 1 − P 1×2 with β = 6/g 2 0 , g 0 being the bare coupling and P 1×1 , P 1×2 the plaquette and rectangular Wilson loops respectively. The configurations were generated with the tree-level Symanzik improved action [23], i.e. b 1 = − 1 12 , b 0 = 1 − 8b 1 . Note that for smoothing the gauge fields using the gradient flow and cooling, in addition to the tree-level Symanzik improved action, we use the Wilson plaquette action which corresponds to b 1 = 0 and b 0 = 1 as well as the Iwasaki improved action with b 1 = −0.331 and b 0 = 3.648.
The fermionic action for the light quarks is the Wilson twisted mass action [24][25][26][27], given in the so-called twisted Table 1 Parameters of the employed ETMC N f = 2 gauge field configuration ensembles [20][21][22]. The columns contain: the inverse bare coupling β, the approximate values of the lattice spacing a [28][29][30], r 0 /a [28,30], the scheme-and scale-independent renormalization con-stants ratio Z P /Z S [31][32][33], the lattice size (L/a) 3 × (T /a), the bare twisted light quark mass aμ l , the critical value of the hopping parameter (where the PCAC mass vanishes), physical extent L of the lattice in fm and the product m π L xχ l (x) D W +m 0 +iμ l γ 5 τ 3 χ l (x) , (2) 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, ψ = exp(iωγ 5 τ 3 /2)χ , with ω being the twist angle (ω = π/2 at maximal twist). The bare untwisted and twisted quark masses are, respectively, m 0 and μ l , while the multiplicatively renormalized light quark mass is μ R = Z −1 P μ l . D W is given by: where ∇ μ and ∇ * μ represent the forward and backward covariant derivatives, respectively.
The details of the gauge field configuration ensembles that were used in this work are shown in Table 1. Most of our investigations are performed using the ensemble b40. 16. However, we also investigate how the correlation between different topological charge definitions changes when approaching the continuum limit, thus using also ensembles c30.20, d20.24 and e17.32. For all of these ensembles the pion mass is close to 340 MeV

Definitions of the topological charge
In this section, we introduce the definitions of the topological charge that we use below for numerical studies. We attempted to include the most commonly used, theoretically sound definitions of the topological charge. The relevant characteristics of each definition are summarized in Table 2. 3.1 Index of the overlap Dirac operator For many years, it was considered impossible that chiral symmetry can be realized on the lattice without violating certain essential properties, like locality, translational invariance and the absence of doublers. This feature of lattice Dirac operators was formulated in terms of a no-go theorem, the Nielsen-Ninomiya theorem [39]. Only after several years, it was realized that this theorem can be overcome by allowing a modified definition of chiral symmetry on the lattice. It was shown by Lüscher [40] that if the lattice Dirac operator satisfies the so-called Ginsparg-Wilson relation [41], there is a corresponding exact symmetry and this symmetry becomes just the standard chiral symmetry in the continuum limit. Thus, any Dirac operator that satisfies the Ginsparg-Wilson relation is chirally symmetric. One of such operators is the overlap Dirac operator, introduced by Neuberger [42,43].
The overlap operator, as a chirally symmetric Dirac operator, can have exact zero modes [44]. The famous Atiyah-Singer index theorem [45] relates in a simple way the number of these zero modes to the topological charge Q of a given gauge field configuration: where n ± denotes the number of zero modes with positive/negative chirality. This remarkable result, thus, links a property of gauge fields to a fermionic observable. By construction, it gives integer values of Q. Note, however, that the definition of the overlap operator is not unique -it depends on the details of the construction of the operator. In common notation, the massless overlap operator is with D W being the standard Wilson-Dirac operator, given by Eq. (3). The s parameter, appearing in the kernel operator A, can be tuned to optimize locality properties of the overlap operator D [46][47][48]. It effectively introduces a dependence of the index obtained on a given configuration on the used value of s. This dependence vanishes towards the continuum limit,   ) cFT HYP (GF Wplaq 3t 0 ) G but at practically used lattice spacings, Q evaluated from the zero modes of the overlap operator shows a dependence on the value of the parameter s. In a sense, this reflects the general property that topology is uniquely defined only for continuum gauge fields. In Sect. 4, we will comment more on the dependence of Q on s by explicitly comparing results obtained for different values of the latter. The overlap index definition of the topological charge is theoretically clean and very appealing, because it provides integer values of Q, while for most other definitions discussed in this paper, the Q values at non-zero lattice spacing are driven away from integers by cut-off effects, ultraviolet fluctuations and/or stochastic noise. However, it has a severe practical drawback -the cost of using the overlap operator is around one to two orders of magnitude larger than the one of e.g. variants of Wilson fermions [49].

Wilson-Dirac operator spectral flow
Closely related to the overlap index is the index derived from the spectral flow of the hermitian Wilson-Dirac operator [50,51]. Its definition is derived from the fact that the continuum hermitian Euclidean Dirac operator H (m 0 ) = γ 5 D(m 0 ) with bare mass m 0 = 0 has a gap, i.e. it has no eigenvalues in the region (− |m 0 |, |m 0 |). As a consequence, eigenvalues crossing zero in the spectral flow of H (m 0 ) can only occur at m 0 = 0, i.e. they correspond to the zero modes of D, and, hence, the net number of crossings is related to the topological charge of the background gauge field [45].
On the lattice, zero crossings in the spectral flow of the hermitian Wilson-Dirac operator H W (m 0 ) = γ 5 (D W + m 0 ) can occur for any value of m 0 in the region −8 ≤ m 0 ≤ 0 [52,53] and counting the net number of crossings in the region −(1 + s) ≤ m 0 ≤ 0 enables one to associate an index to the Wilson-Dirac operator as a function of s. The interpretation from the overlap formalism is essential to make this connection [51]. In fact, the correspondence between the index of the overlap operator and the index from the spectral flow is exact, so the parameter s in Eq. (5) is the same as the one used here, and all the good properties of the overlap index carry over to the index from the spectral flow.
To be more specific, we consider the hermitian Wilson-Dirac operator and its eigenvalues λ H W k (m 0 ). Their dependence on m 0 defines the spectral flow. Since the Wilson-Dirac operator D W is non-normal, the eigensystems of H W (m 0 ) and D W + m 0 are related in a non-trivial way, except for the modes of H W (m 0 ) which are zero for a particular value of m 0 , It follows from this equation that the real modes λ W k ∈ R of D W correspond to zero modes of H W (m 0 = −λ W k ), while the chirality of the modes is given through first order perturbation theory by the derivative of the spectral flow at m 0 = −λ W k [50], Finally, summing up the chiralities of the real modes λ W k ∈ R of D W up to 1 + s yields an index of the Wilson-Dirac operator, and hence the topological charge from the spectral flow, where the sum is over λ W k < 1 + s only, i.e. it excludes the real doubler modes.
Smoothing the gauge fields in the covariant derivatives of D W reduces the non-normality of D W , and hence improves the chirality of the real modes [47]. In addition, it also improves the separation of the physical modes from the doubler modes and in this way reduces the ambiguity of the charge definition due to the choice of the parameter s. Interestingly, this ambiguity can be quantified in the context of Wilson Random Matrix Theory [54][55][56][57][58].

Spectral projectors
Another fermionic definition of the topological charge was introduced by Giusti and Lüscher [9,59]. One introduces the projector P M to the subspace of eigenmodes of the squared Hermitian Dirac operator D † D with eigenvalues below M 2 . The projectors P M can be evaluated stochastically and for chirally symmetric fermions, the topological charge can be defined in terms of it as Q = Tr {γ 5 P M }. This definition is then equivalent to the index definition, apart from the fact that the counting of modes proceeds stochastically, instead of determining it from zero modes. For non-chirally symmetric fermions, the chirality of modes is no longer ±1, but it can be schematically written as ±1 + O(a). Thus, the above definition of Q still holds, but it gives in general non-integer values, contaminated by cut-off effects and by noise from the stochastic evaluation. In practice, the spectral projector computation of the topological charge proceeds in the following way for Wilson-type fermions [9]. One introduces in the theory a set η 1 , . . ., η N src of N src pseudofermion fields with the action S η = N src j=1 (η j , η j ), where the bracket denotes the scalar product. The fields are generated randomly, thus their gauge ensemble average is distributed according to the introduced action. Then, one defines the observable which plays the role of the topological charge. To compute the topological susceptibility from this definition, one needs a correction to account for a finite number of stochastic noise samples N src and the ratio of renormalization constants Z P /Z S . This is done using other observables Then, the topological susceptibility, χ , is defined as If the ratio Z P /Z S is known from another computation, one can replace A 2 / B 2 in the above equation with Z 2 S /Z 2 P . We can, therefore, define as a proxy of the topological charge the quantity with Q = lim N src →∞ Q eff . It can be shown that the spectral projector definition is manifestly ultraviolet finite [11,[59][60][61] and hence theoretically very appealing, especially for the computations of the topological susceptibility, as done in Refs. [9,60,[62][63][64][65]. However, if the aim is to e.g. separate topological sectors, i.e. to choose configurations from a given sector, then the stochastic noise present in the spectral projector evaluated observables strongly contaminates the results, if using a relatively small N src ≈ 6, while for large N src → ∞, the method becomes expensive. As we show in the results section, the stochastic ingredient also makes the correlation with respect to other definitions only moderate and much smaller than e.g. the correlation between the field theoretic definitions (evaluated with different kinds of smearing).
The obtained result also depends on the spectral threshold M chosen for the projector P M . As stated in Ref. [59], M can be chosen arbitrarily, but it is wise to avoid its large values in lattice units (that enhance cut-off effects) and also values close to the quark mass. We will check a few values of M and investigate the implications of choosing different values.
Recently, the spectral projector formalism was extended to staggered fermions [66]. In this case, the bare topological charge is also multiplicatively renormalizable, but the renormalization constants are different, which is due to a different pattern of chiral symmetry breaking.

Field theoretic definition
The topological charge of a gauge field can be naturally defined as the four-dimensional integral over space-time of the topological charge density. In the continuum, this reads where q(x) denotes the topological charge density defined as On the lattice, one has to choose a valid discretization q L (x) of q(x) in order to evaluate Eq. (15), which now takes the form of the sum In practice, any discretization which gives the right continuum limit can be used for the evaluation of Eq. (17), but depending on the discretization q L (x), lattice artifacts affecting the total topological charge Q can vary. This means that using such an operator, on smoothed configurations as we explain in Sect. 3.5, we do not expect to obtain an exact integer value for the total topological charge Q, but rather that the obtained value of Q would be approaching an integer as we tend to the continuum limit i.e. Q = integer ± O(a 2 ). In addition, we expect that the total topological charge for some definitions of q L (x) converges faster and closer to an integer than that obtained by others. One can build such operators from closed path-ordered products of links which lead to the field strength tensor F μν if we perturbatively expand them in a. Namely, by using a number of different Wilson loop shapes and sizes, we cancel, step by step, the leading lattice artifacts contributions. Examples of such operators are demonstrated in the next paragraph. The simplest lattice discretization of q L is based on the simple plaquette and can be noted as with where the square pictorializes the path ordered product of the links lying along plaquette sides in the directionsμ and ν. This definition of q L (x) has a low computational cost and leads to lattice artifacts of order O(a 2 ). Furthermore, it has been used in several determinations of the topological susceptibility and investigations of the instanton properties [67,68]. Without question, the most commonly used definition of q L is the symmetrized clover leaf noted as with Like in the plaquette definition, clover includes lattice artifacts of order O a 2 . This can be viewed easily by perturbatively expanding C plaq μν (x) and C clov μν (x) and obtaining 1 + a 4 F μν (x) + O(a 6 ). Nevertheless, one can also construct improved definitions of topological charge density operators by including additional Wilson loop shapes in the definition of q L (x) and then perturbatively canceling the terms which contribute to higher powers of a. Such a definition is the Symanzik tree-level improved expressed as with and q rect L (x) is the clover-like operator where instead of squares we make use of horizontally and vertically oriented rectangular Wilson loops of size 2 × 1. We remove the discretization errors at tree-level using the Symanzik tree-level coefficients b 1 and b 0 as these were previously used in Eq. (1). Thus, this definition of the topological charge density, by a semiclassical inspection, converges as O(a 4 ) in the continuum limit 3 . Hence, a way to obtain topological quantities with small lattice artifact contributions is by using improved topological density operators. 4 Ultraviolet fluctuations of the gauge fields entering in the definition of the topological charge density lead to contamination of the topological charge. Hence, we employ methods to suppress these UV fluctuations. Such techniques include the gradient flow, the extensively used cooling and several smearing schemes such as APE, HYP and stout. We examine all the above smoothers and investigate their analytic as well as numerical relations.
It is worth mentioning that one can recover the correct definition of topological quantities from unsmoothed configurations by subtracting additive and multiplicative renormalization constants at zero smoothing step. However, the extraction of these renormalization constants may be subject to more systematic effects. Of course, this method is proven useful for particular studies such as the investigation of the θ -dependence of quantities extracted on configurations produced with the topological density weighted by an imaginary theta term iθ being included in the QCD action [70]. Although this method of defining the topological charge is out of the scope of this work, the reader can find more information in Refs. [71] and [72].

Smoothing procedures
Smoothing a gauge link U μ (x) can be accomplished by its replacement by some other link that minimizes a local gauge action. To this purpose, it makes more sense to rewrite the lattice gauge action as where X μ (x) is the sum of all the path ordered products of link matrices, called the "staples", which interact with the link U μ (x). If we consider the Wilson gauge action, the main components in X μ (x) are the staples extending over 1 × 1 squares (in lattice units). We can, therefore, write X μ (x) as According to the above equation, for a given link U μ (x), the total number of plaquette staples interacting with it is 6. There are several ways to iteratively minimize the local action. These include procedures such as the gradient flow, cooling, APE and stout smearing, which make use of the original staples X μ (x) to minimize the local action; this provides the opportunity to perturbatively relate these smoothers and obtain a more concrete understanding on the numerical equivalence among them. Furthermore, other more sophisticated smearing procedures such as HYP smearing, which makes use of a more complicated construction of staples, also exist. However, the latter, although it leads to numerical equivalence with the other smoothing techniques, prohibits us from relating it perturbatively at tree-level order with other smoothing techniques. In the next subsections, we give a brief overview of these most commonly used smoothing techniques used for the calculation of the topological charge.

Gradient flow
Modern non-perturbative studies of QCD have been employing the gradient flow, which has been proven to be a perturbatively and numerically well-defined smoothing procedure. It has good, perturbatively proven renormalization properties and the fields which have been smoothed via gradient flow do not need to be renormalized. From a historic point of view, the gradient flow is related to the streamline idea of Refs. [74][75][76] and its lattice counterpart was previously introduced in the context of Morse theory [77]. The gradient flow is defined as the solution of the evolution equations [12][13][14][15][16] where τ is the dimensionless gradient flow time. In the above equation, the link derivative is defined as with and T a (a = 1, · · · , 8) the Hermitian generators of the The last equation provides all we need in order to smooth the gauge fields according to the Eqs. (27). Evolving the gauge fields via the gradient flow requires the numerical integration of Eqs. (27) manifested by an integration step . This is performed using the third order Runge-Kutta scheme, as explained in Ref. [16]. We set = 0.01 for the integration step, since this has been shown to be a safe option [18]. For the exponentiation of the Lie-algebra fields required for the integration, we apply the algorithm described in Ref. [78]. An important use of the gradient flow is the determination of a reference scale t 0 , defined as the gradient flow time t = a 2 τ in physical units for which where E(t) is the action density To evaluate numerically E(t), we use the clover discretization similarly to Eq. (20).
Having defined the reference scale t 0 , a question is raised: for what value of t shall we read an observable? It has been argued that cut-off effects in some observables can be reduced by chosing larger flow times as reference length scales [79,80]. We therefore chose to commit a numerical comparison of the topological charge obtained with gradient flow with that extracted using different smoothers for three different gradient flow times, namely, t 0 , 2t 0 and 3t 0 .

Cooling
The smoothing technique of cooling [81][82][83][84] was one of the first methods used to remove ultraviolet fluctuations from gauge fields. Cooling is applied to a link variable U μ (x) ∈ SU (3) by updating it, from an old value U old μ (x) to U new μ (x), according to the probability density The basic step of the cooling algorithm is to replace the given link U old μ (x) by an SU (3) group element, which minimizes locally the action, while all the other links remain untouched. This is done by choosing a matrix U new In the case of an SU (2) gauge theory, the maximization is achieved by For SU (3), the maximization can be implemented using the Cabibbo-Marinari algorithm [85] according to which one has to iterate the maximization over all the SU (2) subgroups embedded into SU (3).
We iterate this procedure so that all the links on all sites are updated. Such a sweep over the whole lattice is called one cooling step n c = 1. Traditionally, during such a sweep the link variables, which have already been updated, are subsequently used for the update of the links still retaining their old value. Nevertheless, one can also consider to use the updated links only after the whole lattice is covered, increasing the CPU time by a factor of two.

APE (Array Processor Experiment) smearing
An alternative way to smooth the gauge fields is to apply APE smearing [86] on the gauge configurations. According to this smoothing procedure, we create fat links by adding to the original links the neighbouring staples weighted by a relative strength α APE , which represents the smearing fraction and can be tuned according to its use. This operation breaks unitarity of the resulting "fat" links and shifts them away from det U = 1, thus, we should project back to SU (3). The above operation is noted as The APE smearing scheme can be iterated n APE times to produce smeared links. In addition to the simple APE smearing, variations which make use of "chair" and "diagonal" staples have been proposed from time to time. For the purposes of this investigation, we have considered just the simple APE smearing for different values of the α APE parameter: One can project back to SU (3) by maximizing the expression ReTr{Ũ being the unprojected smeared link. Nonetheless, one can project onto SU (3) using other iterative procedures which suggest that there is not a unique way to do so. This subtlety, however, becomes irrelevant in analytic smearing schemes such as the stout.

Stout smearing
A method which allows analytical derivation of smoothed configurations in SU (3) is the so-called stout smearing proposed in Ref. [78]. This smoothing scheme works in the following way. Once again, let X μ (x) denote the weighted sum of the perpendicular staples which begin at lattice site x and terminate at neighboring site x + aμ. Now, we give a weight ρ st to the staples according to The weight ρ st is a tunable real parameter. Then, the matrix is Hermitian and traceless, where Thus, we define an iterative, analytic link smearing algorithm in which the links U (n st ) μ (x) at stout smearing step n st are mapped into links U (n st +1) μ (x) at stout smearing step n st + 1 using This step can be iterated n st times to finally produce link variables in SU (3) which we call stout links. The structure of the stout smearing procedure resembles the exponentiation steps of the gradient flow and, as a matter of fact, the gradient flow using the Wilson gauge action can be considered as a continuous generalization of stout smearing [87], using gauge paths more complicated than just staples. Hence, for a small enough lattice spacing, one would expect the two smoothing schemes to provide extremely similar results. The level of this similarity is investigated in this work.

HYP (Hypercubic) smearing
We turn now to the discussion of the HYP (Hypercubic) smearing which has been introduced in Ref. [88]. The smeared links of the HYP smoothing procedure are constructed in three steps. These steps are described in the next bullet points.
3. The final step of the HYP smearing consists of applying an APE smearing routine in which the staples are constructed by decorated links which have undergone HYP smearing levels 1 and 2: with The link U (n HYP ) μ (x) is the original link inμ direction which has been smoothed n HYP times.Ũ μ;ν (x) are fat links along the directionμ resulting from the second step of the HYP smearing procedure with the staples extending over directionν being not smeared. The parameter α 3 is tunable and real. 2. The second step of HYP smearing creates the decorated linksŨ μ;ν (x) by applying a modified APE smearing procedure according tõ with Once more, the link U (n HYP ) μ (x) is the link inμ direction which has been smoothed n HYP times. Now, U μ;ρ,ν (x) are fat links along directionμ resulting from the first step of the HYP smearing procedure with staples extending in the directionsρ,ν being non smeared. The parameter α 2 is again tunable and real.
1. The decorated links U ρ;ν,μ (x) are built from the links which have been smeared n HYP using the modified APE smearing step: with In the above expression, only the two staples orthogonal to directionsμ,ν,ρ are included in the smearing procedure.
For the purpose of our work, we choose the values [88]

Perturbative relation between smoothing techniques
The gradient flow, cooling as well as APE, stout and HYP smearing schemes can be used to remove the ultraviolet fluctuations and should lead to the same topological properties, provided that we are close enough to the continuum limit.
Assuming that a is small enough so that we are in the perturbative regime, we can carry out a comparison between the different smoothing procedures in order to obtain an analytic relation among the associated smoothing scales involved, following Refs. [18,19]. Since the gradient flow has the advantage of being the only smoothing scheme with good, perturbatively proven, renormalizability properties, we first relate all other smoothing schemes with the gradient flow and, subsequently, with each other. Gradient flow. The perturbative investigation of the relation between the gradient flow using the Wilson action and cooling has been studied in Ref. [18]. The authors demonstrated that the two smoothing schemes alter the links of a gauge configuration by the same amount if one rescales the flow time and the number of cooling steps according to In the following few lines, we sketch the extraction of the derivative evolution of a gauge link in order to compare it to other smoothing schemes.
In the perturbative regime, a link variable which has been smoothed via the Wilson flow 5 for a finite flow time τ can be expanded as with u a μ (x, τ ) ∈ R assumed to be infinitesimal. Using Eq. (26), the plaquette staples are written as where w a μ (x, τ ) is an infinitesimal quantity. The leading coefficient with the value 6 appearing in the above equation is just the number of plaquettes interacting with the link on which the Wilson flow evolution is applied. We can, therefore, write (52) Hence, Eq. (30) becomes Using the above expression, the evolution of the gradient flow by an infinitesimally small flow time can be approximated as Since the gradient flow is the only smoothing scheme with a concrete theoretical foundation and good renormalizability properties, we will attempt to relate it to other smoothing schemes. Hence, through the resulting matching formulae of any smoothing technique with the gradient flow, we will be able to relate all the different smoothing schemes with each other.
Cooling In the cooling procedure, the link U μ (x, n c ) is substituted with the projection of X μ (x, n c ) over the gauge group. Namely, for the case of the SU (2) gauge theory, this projection is manifested by Eq. (35) where we substitute X μ (x, n c ) by Eq. (51). In the perturbative approximation, this leads to 5 For convenience, here and in the following we refer to the gradient flow using the Wilson gauge action in short as the 'Wilson flow'. This is not to be confused with the spectral flow of the Hermitian Wilson Dirac operator which in the past sometimes has also been called 'Wilson flow'.
The above update corresponds to the substitution Comparing Eqs. (54) and (56), one sees that the flow would evolve the same as cooling if one chooses a step of = 1/6. In addition, during a whole cooling step, the link variables which have already been updated are subsequently used for the update of the remaining links that await update; this corresponds to a speed-up of a factor of two. Therefore, the predicted perturbative relation between the flow time τ and the number of cooling steps n c so that both smoothers have the same effect on the gauge field is given by Eq. (49). This has been shown analytically and demonstrated numerically in Ref. [18]. Moreover, the authors in Ref. [19] have generalized this equivalence for the case of the gradient flow and cooling employing smoothing actions which in addition to the square term multiplied by a factor of b 0 , also included a rectangular term multiplied by In Sect. 4.1, we investigate and confirm that the equivalence for the Wilson smoothing action is manifested by Eq. (49).
APE smearing. We now move to the case of the APE smearing and relate it perturbatively to the Wilson flow and consequently to cooling and stout smearing. Once more, we express the gauge link U μ (x, n APE ) in terms of elements of its Lie algebra Subsequently, we apply this expansion to Eq. (36) and obtain the evolution equation Comparing Eqs. (54) and (59), we observe that the Wilson flow would evolve the gauge links the same as APE smearing if one chooses a flow step of = α APE /6. Hence, the perturbative relation between the Wilson flow time τ and the number of APE smearing steps n APE so that both smoothers have the same effect on the gauge field is The above perturbative matching relation is investigated numerically in Sect. 4.1.
Stout smearing. Let us now turn to the stout smearing smoothing procedure and check whether we could demonstrate analytic equivalence with the Wilson flow and consequently with cooling. According to equation Eq. (38) C μ (x, n st ) can be written as Hence, Ξ μ (x, n st ) is written as The above leads to If we now apply the exponantiation and multiplication according to Eq. (41), we obtain that and in terms of u a Comparing Eqs. (54) and (65), we observe that the Wilson flow would evolve the same as stout smearing if one chooses a step of = ρ st . Therefore, the predicted perturbative relation between the flow time τ and the number of stout smearing steps n st so that both smoothers have the same effect on the gauge field is The above perturbative correspondence is also studied numerically in Section 4.1. The result that we found is in accordance with previous investigations on the matching between stout and APE smearing [89].
HYP smearing Let us now consider the HYP smearing procedure. We have already mentioned that the construction of the decorated staples in the case of HYP smearing prohibits the extraction of a tree-level perturbative relation with the other four smearing procedures. Indubitably, HYP smearing is a valid numerical scheme for smoothing gauge links and removing the ultraviolet fluctuations. Hence, the average action density decreases as we iterate the smoothing procedure. We can, therefore, obtain a numerical equivalence between HYP and another smoother. We attempt to relate the Wilson flow with HYP smearing by calculating the function τ (n HYP ) and interpolating it with an ansatz. The function τ (n HYP ), once more, is defined as the Wilson flow time τ for which the average action density changes by the same amount as when n HYP smearing steps are performed. The perturbative nature of the equivalence and the fact that we need at least three full cooling sweeps to relate the ordinary staple of Eq. (26) with all the components of the HYP staple in Eq. (43), as well as the dependence on three HYP parameters, suggests that the ansatz could be a polynomial such as τ (n HYP ) = An HYP + Bn 2 HYP + Cn 3 HYP .
Since the above equation involves the numerical determination of the associated coefficients, we proceed with this in Sect. 4.1.
Fixing the smoothing scale As the continuum limit is approached, one has to tune the smoothing scale, i.e. gradient flow time, cooling and smearing steps, so that the physics under investigation does not change. When applying a smoothing procedure, the ultraviolet properties of the theory are modified up to some length scale λ S because the ultraviolet fluctuations at smaller length scales are suppressed. In order to have a well defined smoothing procedure towards the continuum limit, one has to make sure that changing the ultraviolet part of the theory leaves the continuum results unaltered, thus, the underlying physics should not depend on λ S . This can be successfully applied by fixing the length scale λ S , which depends on the smoothing parameters. For APE, HYP and stout smearing schemes as well as for cooling, this scale was chosen in the past using different kind of arguments. Nevertheless, for the case of gradient flow, this length scale is quantified. Namely, it has been demonstrated that we can renormalize composite operators at fixed physical length scale related to the flow time by The above equation enables us to translate the length scale λ S to the number of smearing and cooling steps. For cooling, as was shown in Refs. [18,19], we can define the length scale as a function of n c according to the formula For APE smearing, we can write λ S as a function of n APE as In the case of stout smearing, λ S takes the form Finally, for HYP smearing, we can use the numerically extracted τ HYP (n HYP ) in Eq. (67) to define λ S as a function of n HYP according to λ S = a 8τ HYP (n HYP ). Using the above four equations, we can extract a matching relation between two different smoothing schemes with the corresponding matching coefficients given in Table 3.
To explain how we can use the length scale in order to obtain a continuum observable, let us consider the calculation of the continuum limit of the topological susceptibility in quenched QCD. One calculates the topological charge using a given smoothing scheme at a fixed value (in physical units) The value of λ S should be chosen such that it is not too small so that ultraviolet contamination is adequately suppressed, as well as not too large so that the underlying topological structure of the gauge fields is preserved. In most cases, λ S corresponds to a plateau for the topological susceptibility reflecting the scale invariance of the observable. We, therefore, extract the topological susceptibility at fixed λ S for a sequence of lattice spacings and then extrapolate it to the continuum limit.
Of course, the above scenario cannot hold if we consider unquenched QCD, where the physical reference scale depends on the pion mass. In such case, one needs to keep fixed a reference scale such as t 0 using Eq. (31). We can generalize the above procedure for any different smoothing scheme with an effective smoothing flow time defined as a 2 τ (n) where τ (n) corresponds to the matching condition between gradient flow and a given smoothing scheme with smooting scale n.
In our investigation, we evaluated t 0 for the ensemble b40.16 and found that it is equal to t 0 = a 2 τ 0 = a 2 × 2.5; in other words the dimensionless flow time τ 0 = 2.5. For each individual smoother, we can extract a smoothing scale which matches this flow time.

Numerical equivalence between different smoothers
We start the presentation of our results by investigating the perturbative matching between the different smoothing schemes which can be used to remove the ultraviolet diver-gences. We do so by exploring the relation between the average action extracted via both smoothers, looking at the correlation coefficient as well as comparing the topological charge and topological susceptibilities obtained using the two different smoothing schemes. In this section we investigate how the average action density reduces as a function of the two smoothing scales and in the next sections of this paper we investigate the correlation coefficient, the topological charge as well as the topological susceptibility.
Cooling vs. Wilson flow First, we consider the numerical results which have been shown in Refs. [18,19]. We confirm the equivalence realized by Eq. (49) by investigating how the average action density reduces as a function of the two smoothing scales τ and n c . In the left panel of Fig. 1 we show the function τ (n c ) defined as the Wilson flow time τ for which the average action density changes by the same amount as when n c cooling steps are performed. After a few cooling steps, the data appear to lie on the line τ = n c /3. Furthermore, in the right panel of Fig. 1, we present the average action density as a function of τ or the perturbatively determined values of cooling step n c /3. We observe that after approximately 20 cooling steps, the two sets of data coincide. This confirms that the relation τ = n c /3 leads to equivalent results for the average action density between the Wilson flow and cooling for small values of τ and n c .

APE smearing vs. Wilson flow
We move now to the investigation of the numerical equivalence between the APE smearing and the Wilson flow. To test the formula of Eq. (60), we smoothed the gauge configurations via APE smearing for three different values of α APE , namely α APE = 0.4, 0.5 and 0.6. Subsequently, we calculated the function τ (α APE , n APE ) defined as the Wilson flow time τ for which the average action density reduces by the same amount as when n APE smearing steps, for a fixed value of α APE , are performed. In the left panel of Fig. 2, we demonstrate τ (α APE , n APE ) for the three different values of α APE . The data points for the three values of α APE appear to agree with the lines τ = (α APE /6) × n APE providing strong evidence that Eq. (60) provides the right rescaling for which Wilson flow and APE smearing become numerically equivalent. Furthermore, in the right panel of Fig. 2, we provide the average action density as a function of τ and (α APE /6) × n APE demonstrating that the four sets of data perfectly agree with each other.

Stout smearing vs. Wilson flow
We now move to the numerical correspondence between the two smoothing schemes of stout smearing and the Wilson flow. In order to test this HYP smearing vs. Wilson flow As we have already mentioned in Section 3.5.6, the peculiar construction of the HYP smearing staples prohibits the extraction of a linear pertur-bative rescaling between the Wilson flow time τ and the number of HYP smearing steps n HYP . Thus, instead, we attempted a numerical fit using a parametrization of τ (n HYP ) in n HYP according to Eq. (67). In the left panel of Fig. 4, we provide the function τ (n HYP ). Obviously, the sketched behaviour deviates from a linear response (green line) such as those observed for cooling, stout and APE smearing. This suggests that it is impossible to extract a tree-level perturbative expression which relates this smoother with the others. We fit the data using Eq. (67) and extract the coefficients A = 0.25447(32), B = −0.001312(90) as well as C = 1.217(91) × 10 −5 . Of course, these numbers depend on the parameters α HYP1, 2, 3 . In the right panel of Fig. 4, we show the average action density for the HYP smearing as a function of the rescaling equation of Eq. (67) as well as the average action density for the Wilson flow. Clearly, the two lines coincide, demonstrating the realization of a numerical equivalence through Eq. (67).

Field theoretic topological charges on a single configuration
The behaviour of the topological charge Q for single configurations as a function of the gradient flow time τ and for matched smoothing scales for cooling, APE, stout as well as Strikingly, the topological charge obtained with APE smearing as well as stout smearing appears to exhibit significant agreement with that extracted via the Wilson flow. This, of course, occurs if n APE and n st are rescaled according to Eq. (60) and Eq. (66), respectively. Namely, the approximate plateaus observed in Fig. 5 for the three smoothing schemes appear to coincide. An interesting phenomenon is the fine structure occuring when a small instanton or anti-instanton (dislocations) start to drop off the lattice (in case one considers the semiclasical instanton picture). For instance, in the uppermost panel of Fig. 5, and between τ = 6 − 8, the approximate plateau shifts from Q 2 to Q 3. During this transition, we observe that the topological charge Q between the two smoothing schemes and the Wilson flow diverge. Nonetheless, this disagreement appears to vanish as we choose smaller values of ρ st and α APE . In addition, we observe that Q obtained via stout smearing is closer to Wilson flow than APE. The above suggest that the effect of APE and stout smearing on the gauge fields resemble, to a high extent, the Wilson flow. In fact, one does not expect two smoothing procedures to provide equal topological charges since different smoothers carry different lattice artifacts and do not need to agree at non-zero values of the lattice spacing. Indubitably, the topological charges will become closer as the lattice spacing decreases. Thus, as one approaches the continuum limit, any two different procedures converge. Nevertheless, for APE and more strikingly for stout smearing and at finite lattice spacing, the topological charge is, in good approximation, equal to the value extracted via the Wilson flow. We note that the topological charge itself is not the main quantity of interest -the physically relevant observable is the topological susceptibility, which measures the fluctuations of the topological charge.
Turning now to cooling, we demonstrate that the topological charge for a given configuration yields not necessarily the same values as the gradient flow even if we rescale n c according to Eq. (49). Of course, this is not a new observation. Similar comparison which reveals such a possible difference has been published in [19]. Once again, we emphasize that the topological charge for both smoothing procedures will become equivalent if one decreases enough the lattice spacing. As we already know from Ref. [19], both smoothers yield approximately the same topological susceptibility although Finally, we discuss the comparison of the results on Q obtained with HYP smearing and the Wilson flow. Similarly with cooling, if we rescale n HYP with the numericallyextracted formula of Eq. (67), the topological charge Q exhibits approximate equivalence with Q resulted by the Wilson flow for some short range of low values of τ ; however, for this range of τ the value of Q we get is highly dominated by the UV noise. The structure of HYP smearing includes staples which extend beyond the nearest neighbouring links to the original link. This may lead to a supposition that the topological charge obtained via HYP would differ enough from that obtained by the other four smoothers. Interestingly, this does not occur in a noteworthy manner. Of course, once more, the fluctuations of the topological charge are just a measure of the topological susceptibility. Hence, it would be interesting to investigate the response of HYP in this physical quantity; we discuss this topic in Sect. 4.6.

Monte Carlo histories and distribution histograms
In this subsection, we show some typical features of the topological charge evaluated with one representative fermionic  Fig. 6. All relevant topological sectors seem to be scanned correctly and there are no excessive autocorrelations. For the latter, we used the bootstrap procedure with blocking and we find that for different definitions, measurements with a step of 5 configurations (10 Monte Carlo trajectories, i.e. after each saved trajectory, one was unsaved in the generation) yield an integrated autocorrelation time τ int ∈ [0.5, 2] for the topological charge in units of measured configurations. The lowest autocorrelation is obtained obviously for the field theoretic definition without smearing, since one then basically observes uncorrelated ultraviolet fluctuations. All meaningful definitions yield compatible autocorrelation times with τ int ≈ 1.7 (3). The correlation between the values of Q from the two definitions in Fig. 6 is obvious even without com- 6 Mind that here we do not match the smoothing scales, but choose typical smoothing levels applied for the overlap definition and the gluonic definition. The smoothing in both cases serves very different purposes. For the gluonic definition, smoothing is mandatory, since the topological charge from unsmeared configurations is dominated by ultraviolet noise and a meaningful extraction of the charge is not possible. In the case of the overlap index definition, the topological charge can be meaningfully extracted even without any smoothing and the latter is applied only to reduce the computational time. puting the correlation coefficient (which is 88%; see the next subsection for a systematic analysis of correlations between different definitions). This is also illustrated in a scatterplot (Fig. 7). Although the correlation is evident in this plot, it demonstrates that the value of the topological susceptibility is larger for the index definition, indicating that cut-off effects affect the two definitions in a somewhat different manner. Next, we show typical histograms (Fig. 8) obtained with a fermionic definition (again, index HYP1 s = 0) and a gluonic one (Wilson flow at flow times t 0 , 2t 0 and 3t 0 ). For the index definition, we obtain a distribution that is compatible with a Gaussian 7 . For the field theoretic definition, we used an interval width of 0.1 to detect clustering around values close to an integer.
When the flow time is relatively small, of the order of t 0 , basically no clustering is observed. However, when increasing the flow time, at 2t 0 and 3t 0 , the filtering out of the ultraviolet noise is enough to discern peaks at positions close to 0.9, 1.8, 2.7 etc. (for the ensemble b40.16; for other lattice setups or other discretizations of the field stength tensor, the values can be different, but they will also be multiples of some number relatively close to 1). These non-integer positions of the peaks are sometimes "corrected" as mentioned in footnote 4, but this procedure is artificial. In particular, it is not needed to obtain the correct value of the topological susceptibility in the continuum limit. What is, however, relevant for a correct continuum limit is that the smearing procedure 7 It is worth to mention that the distribution is not expected to be ideally Gaussian. For an investigation of non-Gaussianities in the quenched case, see Ref. [70,71,90,91]. However, the detection of such non-Gaussianities requires very large statistics, at least an order of magnitude larger than in our present work. defines a proper smoothing scale, as discussed in the previous section. Such a scale is naturally defined in the gradient flow procedure and in the other smoothing schemes via the matching to gradient flow. If one applies e.g. APE smearing without proper matching to GF, one can not define a consistent procedure of extrapolating to the continuum limit. The traditional method of looking for a plateau in the smearing history is not enough, as it does not define a valid smoothing scale. However, if APE smearing (or any other non-GF type of smoothing) is matched to GF, such a smoothing scale is well-defined and one expects the proper continuum limit for the topological susceptibility.

Correlations between different definitions
For a complete comparison of as many definitions of the topological charge as possible, we concentrated on our ensemble b40.16, i.e. one with the smallest lattice volume and hence the smallest cost of the computations. For this ensemble, we took into account all of the definitions listed in Table 2.
We focus on the correlations between different definitions, expressed by the standard correlation coefficient, normalized to be in the interval [− 1, 1]. We used the bootstrap procedure (with blocking) to compute the error and the influence of autocorrelations.
To understand better the relations between all definitions, we discuss below correlations between different groups of definitions. We start with a general comparison, including the typical representatives of each family from Table 2. Then, we concentrate on the fermionic definitions. The comparison between the most abundant family of field theoretic definitions with several types of smearing that can be applied to the gauge fields to filter out the ultraviolet noise as well as different smoothing actions for the gradient flow is provided in the two sections in the Appendix.

Main comparison
In this subsection, we choose the following definitions as typical representatives: -index of the overlap Dirac operator applied to nonsmeared and smeared gauge fields (with one iteration of HYP smearing) (definitions 1, 3 from Table 2), -spectral flow of the Wilson-Dirac operator computed on gauge fields with one or five iterations of HYP smearing (4, 6), -spectral projectors with two values of the threshold parameter M (9, 10), -field theoretic without smearing (12), -field theoretic with GF at flow time t 0 , three types of smoothing action (16,22,25), -field theoretic with cooling matched to GF at flow time t 0 , three types of smoothing action (28,30,32), -field theoretic with stout smearing matched to GF at flow time t 0 (34), -field theoretic with APE smearing matched to GF at flow time t 0 (40), -field theoretic with HYP smearing matched to GF at flow time t 0 (44).
For the field theoretic definitions, we always use the clover discretization of the topological charge operator for this comparison. The effects of using other discretizations will be considered in one of the further subsections.
Our results are summarized in Fig. 9 and Table 4. In general, we observe very high correlations among different definitions of the topological charge, typically between 85% and 100% (the latter for equivalent definitions).
There are two exceptions to this feature. As expected, the field theoretic definition applied to non-smeared configuration measures basically only ultraviolet noise. The correlation coefficient with respect to other definitions is very small, although non-zero, which suggests that even on non-smeared configurations, some residual signal of the topological charge remains (the correlation coefficient as well as the topological susceptibility are non-zero with statistical significance; nevertheless, reliably extracting the susceptibility from nonsmeared gluonic definition is not possible). Nevertheless, smoothing of gauge fields is mandatory in the field theoretic definition to obtain a meaningful result. The second exception is the spectral projector method, which yields a 55-65% correlation with respect to other cases. One reason for this is obviously the stochastic ingredient in the estimation of Q with this method. However, with 12 stochastic sources that were used, this stochastic ingredient is largely, although not completely, eliminated. Apparently, there are other effects which result in the rather moderate correlation -it is very Table 4 Main comparison of selected definitions of the topological charge. The numbers correspond to the numbering given in Fig. 9. We give the correlation coefficient between different definitions and its error (0 means that the error is smaller than 0.005) likely that these are cut-off effects at the considered, relatively coarse, lattice spacing. Also, one should keep in mind that the spectral projector observable C was never intended to be used as a topological charge observable -it was rather introduced for computations of the topological susceptibility, for which the gauge ensemble average and the stochastic correction play an important role. Within the group of highly correlated definitions, we observe that the fermionic definitions are slightly more correlated with themselves than with the gluonic ones. Concerning the correlation of fermionic and gluonic definitions, it is interesting to note that the former are visibly better correlated with field theoretic ones with improved smoothing actionswhile the correlation with GF/cooling with the Wilson plaquette smoother is around 86-88%, the one with the Iwasaki smoothing action is up to 95%. If, however, one considers the spectral flow (or index) computed on configurations with 5 steps of HYP smearing applied, this effect is alleviated and actually the Iwasaki smoother gives a consistent result with the Wilson plaquette one, while tree-level Symanzik improved is slightly more correlated. We also observe that the correlation between fermionic definitions and GF (Wilson plaquette smoother), stout smearing and APE smearing (both matched to GF at flow time t 0 ) is basically the same. Interesting is the case of the correlation of the index/SF with the gluonic definition on HYP-smeared configurations (matched to GF at flow time t 0 , which for this ensemble implies 10 steps of HYP smearing). It is systematically higher than the one to stout/APE and follows the pattern of the correlation between fermionic and gluonic with GF and the tree-level Symanzik improved smoothing action. In particular, it yields a 97% correlation with SF HYP5. This suggests that HYP smearing has a somewhat similar effect both for fermionic and gluonic definitions. We have also checked the correlation of SF HYP5 and field theoretic with different numbers of HYP iterations and indeed the best correlation is achieved when this number is between 5 and 10 (no statistically significant difference between the latter). This suggests also that some matching between fermionic and gluonic definitions could be achieved.
In the next subsections, we analyze in detail the correlations between definitions inside some selected groups, with specific questions in mind, e.g. about the role of the used discretization of the topological charge operator or about the role of the used smoothing action.

Comparison of fermionic definitions
In this subsection, we make a comprehensive comparison of all fermionic definitions (1-11 from Table 2), see Fig. 10 and Table 5 for a summary. With respect to the main comparison of Sect. 4.4.1, we are able to conclude more about different parameter values that can be used in the definition of Q, i.e. the kernel parameter s for overlap and spectral flow, the number of HYP smearing steps applied before the calculation of Q and different values of the threshold parameter for spectral projectors.
We start by investigating the role of the s parameter of the Wilson-Dirac kernel operator (see Sect. 3.1). It is, in particular, responsible for the locality properties of the overlap Dirac operator [46][47][48] and thus is expected to be important for measurements of Q. For example, if the overlap operator is not local enough, then small topological objects may not be visible (they "fall through the lattice"). Indeed, one can notice sizable effects when s is changed from 0.4 (value that guarantees optimal locality for the non-smeared gauge fields case [48]) to 0 (very bad locality) -the correlation is only around 70%, which is much lower than the correlation with index/SF definitions with good locality properties (in particular, s = 0 for the case of 1 iteration of HYP smearing with 96% correlation). Similarly, the violation of locality in the HYP1 case (change from the optimal value s = 0 to s = 0.75 (more than twice smaller value of the decay rate of the norm of the overlap Dirac operator)) also leads to decreasing correlations. For the case HYP5, locality was not investigated in Ref. [48]. However, from the practically identical results at s = 0 and s = 0.5, one can infer that locality is similar for both of them and guarantees high correlation with respect to the index extracted with the optimally local s = 0 (HYP1) or s = 0.4 (no smearing) values. As stated in Sect. 3.2, the index and spectral flow definitions (with the same value of the s parameter) are exactly equivalent, i.e. should yield a 100% correlation. However, with the spectral flow at a coarse lattice spacing, it may be difficult to disentangle all the zero crossings that determine the value of Q. Similarly, with the index of the overlap operator, numerical precision issues may appear when using too relaxed (to decrease the cost) tolerance criterion for the solver in the procedure of finding zero modes. As a result, the obtained correlation was very close to, but not ideally 1, due to the occurrence of few cases where the value from overlap and from SF differed by ±1.
We now move on to discuss the role of the M parameter for spectral projectors. In Refs. [9,63], the renormalized M parameter (M R ) was set to around 100 MeV (MS scheme at the renormalization scale of 2 GeV). This is expected to be a reasonable choice, since it avoids both the region close to the renormalized quark mass and the region where a M ≈ 1. However, the above references considered relatively large lattices, while the small-volume lattice that we consider for this comparison can suffer from another effect that we discuss below. Namely, with a value of M R ≈ 100 MeV, the number of eigenmodes of the operator D † D is relatively small (typically 5-10), while the number of zero modes can be as high as 15. Hence, one can expect that the projector P M may not include ("count") all the zero modes in some cases, leading to a too small value of the topological susceptibility. To investigate how this feature can affect correlations, we performed the spectral projector calculations for four values of M. Interestingly, we found that the correlations between topological Table 5 Comparison of fermionic definitions of the topological charge. The numbers correspond to the numbering given in Fig. 10. We give the correlation coefficient between different definitions and its error (0 means that the error is smaller than 0.005) Concerning correlations between spectral projectors and index-type definitions, they are typically between 50% and 60%. However, there is no obvious tendency, like an improving correlation when increasing or decreasing M.
In the end, these results provide some warning about interpreting the spectral projector observable C as the topological charge. Although the topological susceptibility is welldefined with spectral projectors, it requires to correct C 2 with B /N when the number of stochastic sources is finite. Together with the presented results for correlations (in particular the ones for index-type definitions, which a priori could be expected to be highly correlated with C), this makes the interpretation of C on a single gauge field configuration difficult. It is, moreover, likely that the values of C extracted with different values of M are affected by different cut-off effects.

Correlation towards the continuum limit
Our next aim is to investigate how the correlations behave towards the continuum limit. The expectation is that very close to the continuum, all definitions agree -hence an increase of correlation coefficients should be observed when decreasing the lattice spacing. We chose one representative fermionic definition (index of the overlap operator evaluated on configurations with 1 step of HYP smearing applied) and one representative gluonic definition (with gradient flow, Wilson plaquette smoothing action at flow time t 0 ). In Fig. 11, we show that the correlation coefficient indeed increases towards the continuum limit, from around 84-88% for the two coarser lattice spacings to 92-93% for the two finer. It is therefore plausible that, as expected, the differences between results at a finite lattice spacing are cut-off effects.

Topological susceptibility
We also compared the topological susceptibility computed with different methods. We defined the topological susceptibility as for all cases, except spectral projectors, where one needs a correction for a finite number of stochastic sources and renormalization with (Z S /Z P ) 2 , see Eq. (13).  Table 2 The comparison is shown in Fig. 12. We find rather nice agreement between different definitions, with most of them giving a value in the range r 0 χ 1/4 ∈ [0.4, 0.5]. It is interesting to observe that the outlying values concern definitions for which certain theoretical doubts appear about their validity. In particular, the index definition for the non-smeared case and s = 0 gives a 20% smaller result than other index definitions. This may be due to the fact that with decreased locality, some topological structures are not counted properly and hence |Q| is too small on some configurations. Similarly, the lowest value of M 2 for spectral projectors might be too small to count all zero modes. However, these effects should go away in the continuum limit -for the index definition, strict locality is then recovered and for spectral projectors, all relevant modes become actual zero modes and hence are counted with any value of M 2 . The situation is somewhat different with the third outlier, the value of χ from the field theoretic definition without any smearing, since by using it one is basically averaging over ultraviolet noise and it is not clear whether any physical signal for the topological susceptibility is left.
The remaining differences between valid definitions, e.g. field theoretic ones with different kinds of smearing, are most likely due to cut-off effects. In particular, changing the smoothing action has a rather strong effect on the computed value of χ . It would certainly be desirable to perform the continuum limit extrapolation for the topological susceptibility from several definitions, but this is inconclusive with the current precision of our data (with typical error of χ at the level of 10%). Nevertheless, with current theoretical understanding, one can be rather certain that the continuum limit is correct for all the cases, excluding the ones for which obvious reservations can be made.
An interesting question regarding the field theoretic definition of Q is how the resulting topological susceptibility behaves as a function of the smoothing scale. Furthermore, we would like to test how the matching between the different smoothers affects χ . To this purpose, in Fig. 13, we present the topological susceptibility extracted using the clover definition of the corresponding charge density for the Wilson flow as a function of the flow time and compare it to other smoothing procedures with smoothing scales adjusted according to Table 3.
In the upper left panel of Fig. 13, we provide a comparison of r 0 χ 1/4 extracted via the Wilson flow and cooling. By rescaling n c → n c /3 according to Eq. (49) we observe that both smoothers give results which agree for the whole range of τ . The above picture resembles similar comparisons demonstrated in Ref. [19]. This suggests that cooling, which is much faster compared to the Wilson flow, provides very similar topological susceptibility as long as the number of cooling steps is rescaled appropriately.
The upper rightmost panel of Fig. 13 reveals an interesting feature of the HYP smoothing technique. Namely, a comparison between HYP smearing with the number of smearing steps rescaled according to Eq. (67) with the Wilson flow shows an approximate agreement (within the statistical accuracy). However, the topological susceptibility via the Wilson flow appears to manifest a plateau, and thus scale invariance, starting at small values of τ ∼ 3. On the contrary, when smoothing with HYP smearing, the plateau sets in at larger values of τ suggesting that this picture could be the outcome of larger cut-off effects in the topological charge. Without question, the non-local character of the HYP smearing introduces different instantonic properties, which could potentially explain this behaviour; this should be investigated in more detail.
In the lower left panel of Fig. 13, we demonstrate a comparison of r 0 χ 1/4 obtained via the Wilson flow as a function of the flow time with r 0 χ 1/4 calculated via stout smearing as a function of the rescaled smearing steps of ρ st n st . We consider three values of ρ st = 0.1, 0.05 and 0.01. Amazingly, the four sketched bands appear to fall on top of each other, yielding a message that the level of similarity between stout smearing and Wilson flow is indeed very high. In fact, the perfect matching of topological susceptibilities in addition to the correlation coefficient of 1.00 flags the exact numerical equivalence between the two smoothers. Once more, this result suggests that we could safely use stout smearing instead of the Wilson flow to measure topological observables and define physical reference scales according to Eq. (71).
Finally, in the lower right panel of Fig. 13, we provide a comparison of r 0 χ 1/4 measured with the Wilson flow as a function of τ with the value of r 0 χ 1/4 extracted with APE smearing vs. the rescaled smearing step α APE n APE /6. We did so for six different values of α APE , namely α APE = 0.6, 0.5, 0.4, 0.3, 0.2 and 0.1. Similarly to the stout smearing presented in the previous paragraph, the topological susceptibility for all six values of α APE appears to match exactly the result of the Wilson flow. Again, this is expected for APE by considering the high similarity of the topological charge revealed in Fig. 5 as well as the correlation coefficient of 1.00 noted in Table 15. Like for stout smearing, one can use APE smearing instead of the Wilson flow with a well defined physical reference scale given by Eq. (70).

Conclusions
In this paper, we have investigated several definitions of the topological charge. Our main conclusion is that all valid 8 definitions lead to a consistent behaviour and are highly correlated, meaning to give -to a large extent-the same topological charge for a given configuration. The progress in recent years, in particular the introduction of the gradient flow smoothing scheme, enabled good control over the topological charge extraction and made it possible to have a well-defined field theoretic definition of the renormalized topological susceptibility and hence a comparatively cheap way to compute the topological susceptibility, including its extrapolation to the continuum limit. This is possible, because the gradient flow provides a valid definition of a smoothing scale, which needs to be kept constant when approaching the continuum limit. The gradient flow also makes it possible to define such a smoothing scale for other kinds of smearing methods, via a well-defined matching procedure. Moreover, one can show that there are no short-distance singularities when the topological charge is defined at a finite flow time. This property, in conjunction with its computational efficiency, makes the gradient flow an excellent choice for computing the topological charge and the corresponding topological susceptibility. We see, however, no obstacle to also use other methods, such as cooling and smearing, for which the number of cooling or smearing steps can be related analytically to the gradient flow time. In this way, one can define a smoothing scale also for the other smoothing schemes and perform the continuum limit by keeping it constant in physical units, which is a prerequisite for a correct continuum limit.
A warning about the usage of field theoretic definitions is provided by our follow-up analysis [92], where we compared the GF definition with the spectral projector one on a wide range of ETMC's N f = 2 + 1 + 1 large-volume ensembles. We found that cut-off effects in the topological susceptibility from the GF are much larger than in the susceptibility from spectral projectors and can be up to 500% at the coarsest lattice spacing of around 0.09 fm. Nevertheless, the continuum limit is always compatible between GF and spectral projectors, and also independent of the flow time (GF) or the renormalized spectral threshold (spectral projectors), as long as these scales are fixed in physical units [93]. We refer the reader to this paper for more details.
In fact, the numerical equivalence introduced by the matching procedure between the gradient flow and other smoothing schemes suggests that in cases where high statistics are needed such as, for example, for the evaluation of higher moments of the topological charge in pure gauge theory, instead of using the gradient flow, one can opt for employing either cooling or APE or stout smearing. From our experience, we find that cooling, APE smearing or stout smearing are, respectively, 120, 20 and 30 times faster 9 than the gradient flow for typical parameter values.
Defining the topological charge and susceptibility using the spectral projectors is another relatively new method. They constitute another theoretically clean way of defining these quantities. However, the cost of this method is significantly larger than the one of the gradient flow. Nevertheless, it might be the method of choice for some applications, since it yields much smaller cut-off effects than the GF, at least in the setup of our follow-up work [92]. Concerning other fermionic definitions, such as the index of the overlap Dirac operator or the spectral flow of the Wilson-Dirac operator, they are theoretically very clean and provide integer values of the topological charge, but their cost is prohibitive for large-scale analyses.
In summary, we have shown in this paper that all valid definitions of the topological charge are highly correlated and, in principle, all of them can be used to analyze topological issues. Thus, the choice for using a certain definition of the topological charge will depend on the particular problem under consideration. under the name of BARYONS. KC was supported in part by the National Science Centre ( Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .

Appendix A: Comparison of gluonic definitions
In this section of the Appendix, we present a comprehensive comparison of the correlation of the topological charge between different field theoretic definitions. We always use the clover discretization, but we vary the type of smoothing procedure, smoothing action (where applicable) and/or other parameters entering the definition of smoothing and flow times: -GF at flow times t 0 , 3t 0 , three types of smoothing action (16,18,22,24,25,27), -cooling matched to GF at flow times t 0 , 3t 0 , three types of smoothing action (28-34), -stout smearing matched to GF at flow times t 0 , 3t 0 , two values of the stout parameter ρ st (34-37), -APE smearing matched to GF at flow times t 0 , 3t 0 , three values of the α APE parameter (38-43), -HYP smearing matched to GF at flow times t 0 , 3t 0 (44)(45).
A summary of our results in shown in Fig. 14 and table 6. We have already discussed the correlations within the class of gradient flow definitions. We now concentrate on comparisons for different groups of definitions and also between the groups.
If one uses cooling as the smoothing procedure, the relation between different cooling times (matched to different flow times) and different smoothing actions for the cooling procedure is very similar to the one for corresponding cases for GF (i.e. with the same smoothing actions). Using stout smearing, one observes that the smoothing parameter ρ st ∈ {0.01, 0.1} has no effect on the correlations. This is natural, since the stout smearing procedure is basically equivalent to the gradient flow with an appropriate step. If this step is small enough, one expects that the results are exact, i.e. the gauge fields evolve according to the continuous gradient flow equations, without any flow time discretization effects. For APE smearing, we also note that the APE parameter α APE has practically no effect on the resulting correlations (we also checked other values than the ones listed in table 2). However, if the number of APE steps is varied, keeping the α parameter fixed, the correlation decreases from 1 to around 0.95. This is qualitatively and quantitatively similar to the GF case and further demonstrates that the matching between GF and APE (or other kinds of smoothing procedure) is very robust. The slight decrease of correlation when going from flow time t 0 to 3t 0 demonstrates that at flow time t 0 , one still has not reached the plateau of Q, i.e. the values of Q, at least for some gauge field configurations, still change when increasing the flow time. However, this effect is much smaller for HYP smearing, where the correlation between the values of Q corresponding to numbers of HYP smearing steps matched to flow times t 0 and 3t 0 is 99%, as compared to typically 95% when using other kinds of smoothing.
Finally, we discuss correlations between different kinds of smoothers. We already argued that GF and stout smearing are equivalent, hence the correlation is perfect if the number of stout smearing steps is matched to the flow time. The correlation between GF and APE smearing is also close to 100% (for matched smoothing scales), while the one between GF/cooling (with the same smoothing action) and GF/HYP smearing is 97%. If one compares definitions at unmatched smoothing scales (e.g. GF at flow time t 0 with APE at a number of steps corresponding to 3t 0 ), one obviously observes a significant decrease of correlation. It is also worth to mention that while for stout, APE or HYP smearing, the notion of a smoothing action does not make sense, still taking into account the way they are constructed, they correspond more to GF with the Wilson plaquette smoother, rather than to GF with more complicated smoothing actions. This implies that the correlation between the values of Q for these smearing types is high with respect to GF with the Wilson plaquette smoothing action (97%-100%), but decreases to a large extent when comparing to GF (or cooling) with the tree-level Symanzik or Iwasaki smoother, to around 90% and 80%, respectively. The very lowest correlation (77%) is observed when comparing GF with the Iwasaki smoother to cooling with the tree-level Symanzik improved action, both at flow time 3t 0 . This correlation is actually even significantly smaller than the correlation of both these cases alone to indextype (fermionic) definitions.

Appendix B: Comparison of different smoothing actions and flow times for gradient flow
The aim of this section is to provide a comparison of the correlation of the topological charge extracted using, different smoothing actions for the gradient flow, as well as different flow times. The included cases are: -Wilson plaquette smoothing action, flow times t 0 , 2t 0 and 3t 0 (16,17,18), -tree-level Symanzik improved smoothing action, flow times t 0 , 2t 0 and 3t 0 (22,23,24), Table 6 Comparison of field theoretic definitions with different kinds of smoothing of UV fluctuations. The numbers correspond to the numbering given in Fig. 14. We give the correlation coefficient between different definitions and its error (0 means that the error is smaller than 0.005) -Iwasaki smoothing action, flow times t 0 , 2t 0 and 3t 0 (25,26,27).
A summary of our findings is given in Fig. 15 and Table 7.
As expected, when the smoothing action is fixed, correlations decrease with increasing difference between corresponding flow times. However, the decrease is very slight and practically invisible in the case of the Iwasaki smoother. This suggests that increasing the flow time has very small effect on the values of Q and the effect is almost absent for the Iwasaki case. When comparing different smoothing actions, one notices that while Wilson plaquette is still very much correlated with tree-level Symanzik improved (92%-97%), the correlation with respect to the Iwasaki smoothing action drops down significantly (to 80%-88%). The correlation of tree-level Symanzik improved to Iwasaki Fig. 16 Comparison of field theoretic definitions with GF smoothing and different discretizations of the topological charge operator. The correlation between different definitions is colour-coded (note the scale is different than in Figs. 9,10,15,14) than Wilson plaquette vs. Iwasaki, but still smaller than the one with respect to Wilson plaquette.
A summary of our findings is given in Fig. 16 and Table 8. Table 7 Comparison of field theoretic definitions with GF smoothing and different smoothing actions and flow times. The numbers correspond to the numbering given in Fig. 15. We give the correlation coeffi-cient between different definitions and its error (0 means that the error is smaller than 0.005)  (1) 1.00(0) 1.00(0) 1 Table 8 Comparison of field theoretic definitions with GF smoothing and different discretizations of the topological charge operator. The numbers correspond to the numbering given in Fig. 16. We give the cor-relation coefficient between different definitions and its error (0 means that the error is smaller than 0.005) The correlations between different discretizations at a fixed flow time are almost perfect (99%-100%) and decrease with an increase of the flow time difference. However, even the largest difference, the one between the simple plaquette discretization at flow time t 0 and the improved one (clover and rectangle terms) at flow time 3t 0 yields a very high correlation (95%). This behaviour is totally consistent with expectations.