1 Introduction

QCD gauge fields can have non-trivial topological properties, manifested in non-zero and integer values of the so-called topological charge. Such topological properties are believed to have important phenomenological implications, e.g. the fluctuations of the topological charge are related to the mass of the flavour-singlet pseudoscalar \(\eta '\) 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 technique 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 \(\tau \) 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 \(\tau =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)Footnote 1. The included definitions are:Footnote 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\),

    • stout smearing with three different values of the stout parameter \(\rho _{\mathrm{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 \(\alpha _{\mathrm{APE}}\) and smearing steps matched to GF time for \(t_0\), \(2t_0\) and \(3t_0\),

    • HYP smearing for a given set of parameters \(\alpha _\mathrm{HYP1}\), \(\alpha _{\mathrm{HYP2}}\), \(\alpha _{\mathrm{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.

Table 1 Parameters of the employed ETMC \(N_f=2\) gauge field configuration ensembles [20,21,22]. The columns contain: the inverse bare coupling \(\beta \), the approximate values of the lattice spacing a [28,29,30], \(r_0/a\) [28, 30], the scheme- and scale-independent renormalization constants ratio \(Z_P/Z_S\) [31,32,33], the lattice size \((L/a)^3\times (T/a)\), the bare twisted light quark mass \(a\mu _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_\pi L\)
Table 2 The relevant characteristics of each topological charge definition. For each definition, we give a number, full name, type of smearing of gauge fields (– = no smearing, HYPn = n iterations of HYP smearing, GF (action,t) = gradient flow with a given smoothing action (Wplaq = Wilson plaquette, tlSym = tree-level Symanzik improved, Iwa = Iwasaki) and at flow time t, cool (action,t) = cooling (smoothing action as for GF) and a number of steps corresponding to GF at flow time t, stout (\(\rho _{\mathrm{st}}\),t) = stout smearing with a given \(\rho _{\mathrm{st}}\) parameter and a number of steps corresponding to GF at flow time t, APE (\(\alpha _\mathrm{APE}\),t) = APE smearing with a given \(\alpha _{\mathrm{APE}}\) parameter and a number of steps corresponding to GF at flow time t, HYP (t) = HYP smearing with a number of steps corresponding to GF at flow time t), short name (used in plots) and definition type (G=gluonic, F=fermionic)

2 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

$$\begin{aligned} S_G[U]= & {} \frac{\beta }{3}\sum _x\Big (b_0 \sum _{\begin{array}{c} \mu ,\nu =1\\ 1\le \mu <\nu \end{array}}^4 \text {Re\,Tr} \big ( 1 - P^{1\times 1}_{x;\mu ,\nu } \big ) \nonumber \\&+ b_1 \sum _{\begin{array}{c} \mu ,\nu =1\\ \mu \ne \nu \end{array}}^4 \text {Re\,Tr}\big ( 1 - P^{1 \times 2}_{x; \mu , \nu } \big ) \Big )\,, \end{aligned}$$
(1)

with \(\beta =6/g_0^2\), \(g_0\) being the bare coupling and \(P^{1\times 1}\), \(P^{1\times 2}\) the plaquette and rectangular Wilson loops respectively. The configurations were generated with the tree-level Symanzik improved action [23], i.e. \(b_1=-\frac{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 basis by

$$\begin{aligned} S_l[\psi , {\bar{\psi }}, U] = a^4 \sum _x {\bar{\chi }}_l(x) \big ( D_W + m_0 + i \mu _l \gamma _5 \tau _3 \big ) \chi _l(x)\,, \end{aligned}$$
(2)

where \(\tau ^3\) acts in flavour space and \(\chi _l=(\chi _u,\,\chi _d)\) is a two-component vector in flavour space, related to the one in the physical basis (\(\psi \)) by a chiral rotation, \(\psi = \exp (i\omega \gamma _5\tau _3/2)\chi \), with \(\omega \) being the twist angle (\(\omega =\pi /2\) at maximal twist). The bare untwisted and twisted quark masses are, respectively, \(m_0\) and \(\mu _l\), while the multiplicatively renormalized light quark mass is \(\mu _R=Z_P^{-1}\mu _l\). \(D_W\) is given by:

$$\begin{aligned} D_W = \frac{1}{2} \big ( \gamma _{\mu } (\nabla _{\mu } + \nabla ^*_{\mu }) - a \nabla ^*_{\mu } \nabla _{\mu } \big ), \end{aligned}$$
(3)

where \(\nabla _{\mu }\) and \(\nabla ^*_{\mu }\) represent the forward and backward covariant derivatives, respectively.

Twisted mass fermions are providing an automatically \({\mathcal {O}}(a)\)-improvement if the twist angle is set to \(\pi /2\) (maximal twist). This can be achieved by non-perturbative tuning of the hopping parameter \(\kappa = (8+2 a m_0)^{-1}\) to its critical value, i.e. such that the PCAC quark mass vanishes [25, 34,35,36,37,38].

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

3 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:

$$\begin{aligned} Q = n_{-} - n_{+}, \end{aligned}$$
(4)

where \(n_\pm \) 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

$$\begin{aligned} D=\frac{1}{a}\left( 1-\frac{A}{\sqrt{A^\dagger A}} \right) ,\qquad A=1+s-aD_W\,, \end{aligned}$$
(5)

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, 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].

3.2 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) = \gamma _5 D(m_0)\) with bare mass \(m_0\ne 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) = \gamma _5 (D_W + m_0)\) can occur for any value of \(m_0\) in the region \(-8 \le m_0 \le 0\) [52, 53] and counting the net number of crossings in the region \(-(1+s) \le m_0 \le 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

$$\begin{aligned} H_W(m_0) = \gamma _5 (D_W + m_0), \end{aligned}$$
(6)

and its eigenvalues \(\lambda _k^{H_W}(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\),

$$\begin{aligned} H_W(m_0) \psi = 0 \qquad \Longleftrightarrow \qquad D_W \psi = -m_0 \psi . \end{aligned}$$
(7)

It follows from this equation that the real modes \(\lambda _k^W \in {\mathbb {R}}\) of \(D_W\) correspond to zero modes of \(H_W(m_0=-\lambda _k^W)\), while the chirality of the modes is given through first order perturbation theory by the derivative of the spectral flow at \(m_0 =-\lambda _k^W\) [50],

$$\begin{aligned} \left. \frac{d \lambda _k^{H_W}}{dm_0} \right| _{m_0=-\lambda _k^W} = \langle k |\gamma _5| k\rangle . \end{aligned}$$
(8)

Finally, summing up the chiralities of the real modes \(\lambda _k^W \in {\mathbb {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,

$$\begin{aligned} Q = \sum _{\lambda _k^W \in {\mathbb {R}}} \text {sign}(\langle k |\gamma _5| k\rangle ), \end{aligned}$$
(9)

where the sum is over \(\lambda _k^W < 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].

3.3 Spectral projectors

Another fermionic definition of the topological charge was introduced by Giusti and Lüscher [9, 59]. One introduces the projector \({\mathbb {P}}_M\) to the subspace of eigenmodes of the squared Hermitian Dirac operator \(D^\dagger D\) with eigenvalues below \(M^2\). The projectors \({\mathbb {P}}_M\) can be evaluated stochastically and for chirally symmetric fermions, the topological charge can be defined in terms of it as \(Q=\mathrm{Tr}\,\{\gamma _5{\mathbb {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 \(\pm 1\), but it can be schematically written as \(\pm 1+{\mathcal {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 \(\eta _1\), \(\ldots \), \(\eta _{N_{src}}\) of \(N_{src}\) pseudofermion fields with the action \(S_\eta =\sum _{j=1}^{N_{src}} (\eta _j,\eta _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

$$\begin{aligned} {{{\mathcal {C}}}}=\frac{1}{N_{src}}\sum _{j=1}^{N_{src}} \left( {\mathbb {P}}_M\eta _j,\gamma _{5}{\mathbb {P}}_M\eta _j\right) . \end{aligned}$$
(10)

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

$$\begin{aligned} {{{\mathcal {B}}}}= & {} \frac{1}{N_{src}}\sum _{j=1}^{N_{src}} \left( {\mathbb {P}}_M\gamma _5{\mathbb {P}}_M\eta _j,{\mathbb {P}}_M\gamma _5{\mathbb {P}}_M\eta _j\right) , \end{aligned}$$
(11)
$$\begin{aligned} {{{\mathcal {A}}}}= & {} \frac{1}{N_{src}}\sum _{j=1}^{N_{src}} \left( {\mathbb {P}}_M^2\eta _j,{\mathbb {P}}_M^2\eta _j\right) . \end{aligned}$$
(12)

Then, the topological susceptibility, \(\chi \), is defined as

$$\begin{aligned} \chi =\frac{1}{V}\frac{\langle {{{\mathcal {A}}}}\rangle ^2}{\langle \mathcal{B}\rangle ^2} \left( {\langle {{{\mathcal {C}}}}^2\rangle -\frac{\langle \mathcal{B}\rangle }{N_{src}}} \right) . \end{aligned}$$
(13)

If the ratio \(Z_P/Z_S\) is known from another computation, one can replace \(\langle {{{\mathcal {A}}}}\rangle ^2/\langle {{{\mathcal {B}}}}\rangle ^2\) in the above equation with \(Z_S^2/Z_P^2\). We can, therefore, define as a proxy of the topological charge the quantity

$$\begin{aligned} Q_{\mathrm{eff}} = \frac{Z_S}{Z_P} {{{\mathcal {C}}}}, \end{aligned}$$
(14)

with \(Q = \lim _{N_{src} \rightarrow \infty } Q_{\mathrm{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} \approx 6\), while for large \(N_{src} \rightarrow \infty \), 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 \({\mathbb {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.

3.4 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

$$\begin{aligned} Q=\int d^4 x \,q(x), \end{aligned}$$
(15)

where q(x) denotes the topological charge density defined as

$$\begin{aligned} q(x) = \frac{1}{32 \pi ^2} \epsilon _{\mu \nu \rho \sigma } \mathrm{Tr} \left\{ F_{\mu \nu } F_{\rho \sigma } \right\} . \end{aligned}$$
(16)

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

$$\begin{aligned} Q=a^4 \sum _{x} q_{L}(x). \end{aligned}$$
(17)

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=\mathrm{integer} \pm 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_{\mu \nu }\) 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

$$\begin{aligned} q_L^{\mathrm{plaq}}(x) = \frac{1}{32 \pi ^2} \epsilon _{\mu \nu \rho \sigma } \mathrm{Tr} \left( C^{\mathrm{plaq}}_{\mu \nu } C^{\mathrm{plaq}}_{\rho \sigma } \right) \,, \end{aligned}$$
(18)

with

(19)

where the square pictorializes the path ordered product of the links lying along plaquette sides in the directions \({{\hat{\mu }}}\) and \({{\hat{\nu }}}\). This definition of \(q_{L}(x)\) has a low computational cost and leads to lattice artifacts of order \({{{\mathcal {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

$$\begin{aligned} q_L^{\mathrm{clov}}(x) = \frac{1}{32 \pi ^2} \epsilon _{\mu \nu \rho \sigma } \mathrm{Tr} \left( C^{\mathrm{clov}}_{\mu \nu } C^{\mathrm{clov}}_{\rho \sigma } \right) , \end{aligned}$$
(20)

with

(21)

Like in the plaquette definition, clover includes lattice artifacts of order \({{{\mathcal {O}}}}\left( a^2 \right) \). This can be viewed easily by perturbatively expanding \(C^{\mathrm{plaq}}_{\mu \nu }(x)\) and \( C^\mathrm{clov}_{\mu \nu }(x)\) and obtaining \(1+a^4 F_{\mu \nu }(x)+ \mathcal{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}\left( x \right) \) and then perturbatively canceling the terms which contribute to higher powers of a. Such a definition is the Symanzik tree-level improved expressed as

$$\begin{aligned} q^{\mathrm{imp}}_{L}(x)=b_0 q_L^{\mathrm{clov}}(x) + b_1 q_L^{\mathrm{rect}} (x), \end{aligned}$$
(22)

with

$$\begin{aligned} q_L^{\mathrm{rect}}(x) = \frac{2}{32 \pi ^2} \epsilon _{\mu \nu \rho \sigma } \mathrm{Tr} \left( C^{\mathrm{rect}}_{\mu \nu } C^{\mathrm{rect}}_{\rho \sigma } \right) , \end{aligned}$$
(23)

and

(24)

\(q^{\mathrm{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 \times 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 \({{{\mathcal {O}}}}(a^4)\) in the continuum limitFootnote 3. Hence, a way to obtain topological quantities with small lattice artifact contributions is by using improved topological density operators.Footnote 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 \(\theta \)-dependence of quantities extracted on configurations produced with the topological density weighted by an imaginary theta term \(i \theta \) 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].

3.5 Smoothing procedures

Smoothing a gauge link \(U_{\mu }(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

$$\begin{aligned} S_G= & {} \frac{\beta }{3} \mathrm{Re} \mathrm{Tr}\{ { X}^\dagger _{\mu }(x) { U}_{\mu }(x) \} \nonumber \\&\ \ \ \ \ \ \ \ \ \ \ \ + \{ \mathrm{terms \ independent \ of} \ { U}_{\mu }(x) \}, \end{aligned}$$
(25)

where \(X_{\mu }(x)\) is the sum of all the path ordered products of link matrices, called the “staples”, which interact with the link \(U_{\mu }(x)\). If we consider the Wilson gauge action, the main components in \(X_{\mu }(x)\) are the staples extending over \(1 \times 1\) squares (in lattice units). We can, therefore, write \(X_{\mu }(x)\) as

$$\begin{aligned} X_{\mu }(x)= & {} \sum _{\nu \ge 0, \nu \ne \mu }\Big [ U_{\nu }(x)U_{\mu }(x+a\hat{\nu })U_{\nu }^{\dag }(x+a\hat{\mu }) \nonumber \\&+ \!U_{\nu }^{\dag }(x-a\hat{\nu })U_{\mu }(x-a\hat{\nu })U_{\nu }(x-a\hat{\nu }+a\hat{\mu })\Big ]. \nonumber \\ \end{aligned}$$
(26)

According to the above equation, for a given link \(U_{\mu }(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_{\mu }(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.

3.5.1 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]

$$\begin{aligned} {\dot{V}}_{\mu } \left( x, \tau \right)= & {} -g_0^2 \left[ \partial _{x, \mu } S_G(V(\tau )) \right] V_{\mu } \left( x, \tau \right) , \nonumber \\ V_{\mu } \left( x, 0 \right)= & {} U_{\mu } \left( x \right) , \end{aligned}$$
(27)

where \(\tau \) is the dimensionless gradient flow time. In the above equation, the link derivative is defined as

$$\begin{aligned} \partial _{x, \mu } S_G(U)= & {} i \sum _{a} T^a \frac{\mathrm{d}}{{\mathrm{d}} s} S_G \left( e^{is Y^a} U \right) \Bigg |_{s=0} \, \nonumber \\\equiv & {} i \sum _{a} T^a \partial ^{(a)}_{x,\mu } S_G(U), \end{aligned}$$
(28)

with

$$\begin{aligned} Y^a(y,\nu )=\left\{ \begin{array}{ll} T^a &{}\quad \mathrm {if}\ (y,\nu )=(x,\mu ), \\ 0 &{}\quad \mathrm {if}\ (y,\nu ) \ne (x,\mu ), \end{array} \right. \ \end{aligned}$$
(29)

and \(T^a\) (\(a=1,\cdots , 8\)) the Hermitian generators of the SU(3) group. If we now set \(\varOmega _{\mu }(x) = U_{\mu }(x) X^{\dagger }_{\mu } (x)\), we obtain

$$\begin{aligned} g_0^2 \partial _{x,\mu } S_G(U) (x) =&\frac{1}{2} \left( \varOmega _{\mu } (x) - \varOmega ^{\dagger }_{\mu } (x) \right) \nonumber \\&- \frac{1}{6} \mathrm{Tr} \left( \varOmega _{\mu } (x) - \varOmega ^{\dagger }_{\mu } (x) \right) . \end{aligned}$$
(30)

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 \(\epsilon \). This is performed using the third order Runge-Kutta scheme, as explained in Ref. [16]. We set \(\epsilon =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 \tau \) in physical units for which

$$\begin{aligned} t^2 \langle E(t) \rangle |_{t=t_0} = 0.3, \end{aligned}$$
(31)

where E(t) is the action density

$$\begin{aligned} E(t)=-\frac{1}{2V} \sum _{x} \mathrm{Tr} \left\{ F_{\mu \nu }(x,t) F_{\mu \nu }(x,t) \right\} . \end{aligned}$$
(32)

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\), \(2 t_0\) and \(3 t_0\).

3.5.2 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_{\mu }(x) \in SU(3)\) by updating it, from an old value \(U^{\mathrm{old}}_{\mu }(x)\) to \(U^{\mathrm{new}}_{\mu }(x)\), according to the probability density

$$\begin{aligned} P(U) \propto \mathrm{exp}\left\{ -\lim _{\beta \rightarrow \infty } \beta \frac{1}{3} \mathrm{Re} \mathrm{Tr}{ X_{\mu }}^\dagger (x) { U_{\mu }}(x) \right\} . \end{aligned}$$
(33)

The basic step of the cooling algorithm is to replace the given link \(U^{\mathrm{old}}_{\mu }(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^{\mathrm{new}}_{\mu }(x) \in SU(3)\) that maximizes

$$\begin{aligned} \mathrm{Re} \mathrm{Tr}\{ { U}^{\mathrm{new}}_{\mu }(x) { X}^{\dagger }_{\mu }(x) \}. \end{aligned}$$
(34)

In the case of an SU(2) gauge theory, the maximization is achieved by

$$\begin{aligned} U^{\mathrm{new}}_{\mu }(x) = \frac{X_{\mu } (x)}{ \sqrt{{\mathrm{det}} X_{\mu } (x)}}. \end{aligned}$$
(35)

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.

3.5.3 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 \(\alpha _{\mathrm{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 \(\mathrm{det}\,U =1\), thus, we should project back to SU(3). The above operation is noted as

$$\begin{aligned} U_{\mu }^{\left( n_{\mathrm{APE}} +1 \right) } \left( x \right)= & {} {Proj}_{SU(3)} \left[ \left( 1 - \alpha _{\mathrm{APE}} \right) U_{\mu }^{\left( n_{\mathrm{APE}} \right) } \left( x \right) \nonumber \right. \\&\left. + \frac{\alpha _{\mathrm{APE}}}{6} X_{\mu }^{\left( n_{\mathrm{APE}} \right) } (x)\right] \,. \end{aligned}$$
(36)

The APE smearing scheme can be iterated \(n_{\mathrm{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 \(\alpha _{\mathrm{APE}}\) parameter:

$$\begin{aligned} \alpha _{\mathrm{APE}} = 0.1, \; 0.2, \; 0.3, \; 0.4, \; 0.5 \;\; \mathrm{and} \;\; 0.6. \end{aligned}$$
(37)

One can project back to SU(3) by maximizing the expression \(\mathrm{Re} \mathrm{Tr}\{ {{{\tilde{U}}}}^{\left( n_{\mathrm{APE}} +1 \right) }_{\mu }(x) { X}^{\dagger }_{\mu }(x) \}\) with \({{{\tilde{U}}}}^{\left( n_{\mathrm{APE}} +1 \right) }_{\mu }(x)\) 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.

3.5.4 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_{\mu } (x)\) denote the weighted sum of the perpendicular staples which begin at lattice site x and terminate at neighboring site \(x+a{{\hat{\mu }}}\). Now, we give a weight \(\rho _{\mathrm{st}}\) to the staples according to

$$\begin{aligned} C_{\mu } (x) = \rho _{\mathrm{st}} X_{\mu } (x). \end{aligned}$$
(38)

The weight \(\rho _{\mathrm{st}}\) is a tunable real parameter. Then, the matrix \(Q_{\mu } (x)\), defined in SU(3) by

$$\begin{aligned} Q_{\mu } (x)= & {} \frac{i}{2} \left( \varXi ^{\dagger }_{\mu }(x) - \varXi _{\mu }(x) \right) \nonumber \\&\quad \quad - \frac{i}{6} \mathrm{Tr} \left( \varXi ^{\dagger }_{\mu }(x) - \varXi _{\mu }(x) \right) \, \end{aligned}$$
(39)

is Hermitian and traceless, where

$$\begin{aligned} \varXi _{\mu }(x) = C_{\mu } (x) U^{\dagger }_{\mu }(x). \end{aligned}$$
(40)

Thus, we define an iterative, analytic link smearing algorithm in which the links \(U^{(n_{\mathrm{st}})}_{\mu }(x)\) at stout smearing step \(n_{\mathrm{st}}\) are mapped into links \(U^{(n_{\mathrm{st}}+1)}_{\mu }(x)\) at stout smearing step \(n_{\mathrm{st}}+1\) using

$$\begin{aligned} U^{(n_{\mathrm{st}}+1)}_{\mu }(x) = \mathrm{exp} \left( i Q^{n_\mathrm{st}}_{\mu } (x) \right) U^{(n_{\mathrm{st}})}_{\mu }(x). \end{aligned}$$
(41)

This step can be iterated \(n_{\mathrm{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.

3.5.5 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.

  1. 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:

    $$\begin{aligned} U_{\mu }^{\left( n_{\mathrm{HYP}} +1 \right) } \left( x \right) = {Proj}_{SU(3)}&[ \left( 1 - \alpha _3\right) U_{\mu }^{\left( n_\mathrm{HYP} \right) } \left( x \right) \nonumber \\&+ \frac{\alpha _3}{6} {{{\tilde{X}}}}_{\mu }(x) ], \end{aligned}$$
    (42)

    with

    $$\begin{aligned} {{{\tilde{X}}}}_{\mu }(x)= & {} \!\!\!\!\!\! \sum _{\nu \ge 0, \nu \ne \mu }\!\!\!\!\!\Big [ {{{\tilde{U}}}}_{\nu ; \mu }(x){{{\tilde{U}}}}_{\mu ; \nu }(x+a\hat{\nu }) {{{\tilde{U}}}}_{\nu ; \mu }^{\dag }(x+a\hat{\mu }) \nonumber \\&+ {{{\tilde{U}}}}_{\nu ; \mu }^{\dag }(x-a\hat{\nu }){{{\tilde{U}}}}_{\mu ; \nu }\nonumber \\&(x-a\hat{\nu }) {{{\tilde{U}}}}_{\nu ; \mu }(x-a\hat{\nu }+a\hat{\mu })\Big ]. \end{aligned}$$
    (43)

    The link \(U_{\mu }^{\left( n_{\mathrm{HYP}} \right) } \left( x \right) \) is the original link in \({{\hat{\mu }}}\) direction which has been smoothed \(n_{\mathrm{HYP}}\) times. \({{{\tilde{U}}}}_{\mu ; \nu }(x)\) are fat links along the direction \({{\hat{\mu }}}\) resulting from the second step of the HYP smearing procedure with the staples extending over direction \({\hat{\nu }}\) being not smeared. The parameter \(\alpha _3\) is tunable and real.

  2. 2.

    The second step of HYP smearing creates the decorated links \({{{\tilde{U}}}}_{\mu ; \nu }(x)\) by applying a modified APE smearing procedure according to

    $$\begin{aligned} {{{\tilde{U}}}}_{\mu ;\nu } \left( x \right)= & {} {Proj}_{SU(3)} [ ( 1 - \alpha _2) U_{\mu }^{\left( n_{\mathrm{HYP}} \right) } \left( x \right) \nonumber \\&\quad \quad \quad \quad \quad + \frac{\alpha _2}{4} {{{\overline{X}}}}_{\mu ; \nu } (x) ], \end{aligned}$$
    (44)

    with

    $$\begin{aligned}&{{{\overline{X}}}}_{\mu ;\nu }(x) \nonumber \\&\quad =\sum _{\rho \ge 0, \rho \ne \nu , \mu }\!\!\!\!\!\!\!\Big [ {{{\overline{U}}}}_{\rho ; \nu , \mu }(x){{{\overline{U}}}}_{\mu ; \rho , \nu }(x+a\hat{\rho }) {{{\overline{U}}}}_{\rho ; \nu , \mu }^{\dag }(x+a\hat{\mu }) \nonumber \\&\qquad + {{{\overline{U}}}}_{\rho ; \nu , \mu }^{\dag }(x-a\hat{\rho })\nonumber \\&\qquad {\overline{U}}_{\mu ; \rho , \nu }(x-a\hat{\rho }) {{{\overline{U}}}}_{\rho ; \mu , \nu }(x-a\hat{\rho }+a\hat{\mu })\Big ]\,.\nonumber \\ \end{aligned}$$
    (45)

    Once more, the link \(U_{\mu }^{\left( n_{\mathrm{HYP}} \right) } \left( x \right) \) is the link in \({{\hat{\mu }}}\) direction which has been smoothed \(n_{\mathrm{HYP}}\) times. Now, \({{{\overline{U}}}}_{\mu ; \rho , \nu }(x)\) are fat links along direction \({{\hat{\mu }}}\) resulting from the first step of the HYP smearing procedure with staples extending in the directions \({{\hat{\rho }}}\), \({{\hat{\nu }}}\) being non smeared. The parameter \(\alpha _2\) is again tunable and real.

  3. 1.

    The decorated links \({{{\overline{U}}}}_{\rho ; \nu , \mu }(x)\) are built from the links which have been smeared \(n_{\mathrm{HYP}}\) using the modified APE smearing step:

    $$\begin{aligned} {{{\overline{U}}}}_{\mu ;\nu , \rho } \left( x \right)= & {} {Proj}_{SU(3)} \left[ \left( 1 - \alpha _1 \right) U_{\mu }^{\left( n_{\mathrm{HYP}} \right) } \left( x \right) \nonumber \right. \\&\left. \quad \quad \quad + \frac{\alpha _2}{2} {\breve{X}}_{\mu ; \nu \rho } (x) \right] , \end{aligned}$$
    (46)

    with

    $$\begin{aligned}&{\breve{X}}_{\mu ;\nu , \rho }(x) \nonumber \\&\quad =\sum _{\sigma \ge 0, \sigma \ne \rho , \nu , \mu }\!\!\!\!\!\!\!\!\!\!\Big [ {U}^{\left( n_{\mathrm{HYP}} \right) }_{\sigma ; \rho , \nu , \mu }(x){U}^{\left( n_{\mathrm{HYP}} \right) }_{\mu ; \sigma , \rho , \nu }(x+a\hat{\sigma })\nonumber \\&\qquad {U^{\dagger }}_{\sigma ; \rho , \nu , \mu }^{\left( n_{\mathrm{HYP}} \right) }(x+a\hat{\mu }) \nonumber \\&\qquad +{U^{\dagger }}_{\sigma ; \rho , \nu , \mu }^{\left( n_\mathrm{HYP} \right) }(x-a\hat{\sigma }){U}_{\mu ; \sigma , \rho , \nu }^{\left( n_{\mathrm{HYP}} \right) }(x-a\hat{\sigma }) \nonumber \\&\qquad {U}_{\sigma ; \rho , \mu , \nu }^{\left( n_{\mathrm{HYP}} \right) }(x-a\hat{\sigma }+a\hat{\mu })\Big ]. \end{aligned}$$
    (47)

    In the above expression, only the two staples orthogonal to directions \({{\hat{\mu }}}\), \({{\hat{\nu }}}\), \({{\hat{\rho }}}\) are included in the smearing procedure.

For the purpose of our work, we choose the values [88]

$$\begin{aligned} \alpha _1 = 0.75, \; \alpha _2 = 0.6, \; \alpha _3 = 0.3. \end{aligned}$$
(48)

3.5.6 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

$$\begin{aligned} \tau \simeq \frac{n_c}{3}. \end{aligned}$$
(49)

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 flowFootnote 5 for a finite flow time \(\tau \) can be expanded as

(50)

with \(u^a_{\mu }(x, \tau ) \in {\mathbb {R}}\) assumed to be infinitesimal. Using Eq. (26), the plaquette staples are written as

(51)

where \(w^a_{\mu }(x, \tau )\) 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 \(\varOmega _{\mu }(x, \tau )\) as

(52)

Hence, Eq. (30) becomes

$$\begin{aligned} g_0^2 \partial _{x,\mu } S_G(U) \simeq i \sum _{a} \left[ 6 u^a_{\mu }(x,\tau ) - w_{\mu }^a(x,\tau ) \right] T^a. \end{aligned}$$
(53)

Using the above expression, the evolution of the gradient flow by an infinitesimally small flow time \(\epsilon \) can be approximated as

$$\begin{aligned} u_{\mu }^a (x,\tau +\epsilon ) \simeq u_{\mu }^{a} (x,\tau ) - \epsilon \left[ 6 u^a_{\mu }(x,\tau ) \!-\! w_{\mu }^a(x,\tau ) \right] . \end{aligned}$$
(54)

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_{\mu }(x,n_c)\) is substituted with the projection of \(X_{\mu }(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_{\mu }(x,n_c)\) by Eq. (51). In the perturbative approximation, this leads to

(55)

The above update corresponds to the substitution

$$\begin{aligned} u^a_{\mu } (x, n_c+1) = \frac{w_{\mu }^a(x,n_c)}{6}. \end{aligned}$$
(56)

Comparing Eqs. (54) and (56), one sees that the flow would evolve the same as cooling if one chooses a step of \(\epsilon = 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 \(\tau \) 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 \(b_1=(1-b_0)/8\). This equivalence in manifested by the formula

$$\begin{aligned} \tau \simeq \frac{n_c}{3-15 b_1}. \end{aligned}$$
(57)

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_{\mu } (x, n_{\mathrm{APE}})\) in terms of elements of its Lie algebra

(58)

Subsequently, we apply this expansion to Eq. (36) and obtain the evolution equation

$$\begin{aligned}&u_{\mu }^a (x,n_{\mathrm{APE}} + 1) \simeq u_{\mu }^{a} (x,n_\mathrm{APE}) \nonumber \\&\quad - \frac{\alpha _{\mathrm{APE}}}{6} \left[ 6 u^a_{\mu }(x,n_{\mathrm{APE}}) - w_{\mu }^a(x,n_{\mathrm{APE}}) \right] . \end{aligned}$$
(59)

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 \(\epsilon = {\alpha _{\mathrm{APE}}}/{6}\). Hence, the perturbative relation between the Wilson flow time \(\tau \) and the number of APE smearing steps \(n_{\mathrm{APE}}\) so that both smoothers have the same effect on the gauge field is

$$\begin{aligned} \tau = \frac{\alpha _{\mathrm{APE}}}{6} n_{\mathrm{APE}}. \end{aligned}$$
(60)

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_{\mu } (x, n_{\mathrm{st}})\) can be written as

(61)

Hence, \(\varXi _{\mu } (x, n_{\mathrm{st}})\) is written as

(62)

The above leads to

$$\begin{aligned} Q_{\mu }(x, n_{\mathrm{st}}) \!\simeq \! - \rho _{\mathrm{st}} \sum _{a} \!\left[ 6 u^a_{\mu }(x, n_{\mathrm{st}}) \!-\! w_{\mu }^a(x, n_{\mathrm{st}}) \right] T^a. \end{aligned}$$
(63)

If we now apply the exponantiation and multiplication according to Eq. (41), we obtain that

(64)

and in terms of \(u_{\mu }^{a} (x,n_{\mathrm{st}})\) that

$$\begin{aligned} u_{\mu }^a (x,n_{\mathrm{st}}+1)\simeq & {} u_{\mu }^{a} (x,n_{\mathrm{st}}) \nonumber \\&- \rho _{\mathrm{st}} \left[ 6 u^a_{\mu }(x,n_{\mathrm{st}}) - w_{\mu }^a(x,n_{\mathrm{st}}) \right] . \end{aligned}$$
(65)

Comparing Eqs. (54) and (65), we observe that the Wilson flow would evolve the same as stout smearing if one chooses a step of \(\epsilon = \rho _{\mathrm{st}}\). Therefore, the predicted perturbative relation between the flow time \(\tau \) and the number of stout smearing steps \(n_{\mathrm{st}}\) so that both smoothers have the same effect on the gauge field is

$$\begin{aligned} \tau = \rho _{\mathrm{st}} n_{\mathrm{st}}. \end{aligned}$$
(66)

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 \(\tau (n_{\mathrm{HYP}})\) and interpolating it with an ansatz. The function \(\tau (n_{\mathrm{HYP}})\), once more, is defined as the Wilson flow time \(\tau \) for which the average action density changes by the same amount as when \(n_{\mathrm{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

$$\begin{aligned} \tau (n_{\mathrm{HYP}}) = A n_{\mathrm{HYP}} + B n^2_{\mathrm{HYP}} + C n^3_{\mathrm{HYP}}. \end{aligned}$$
(67)

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 \(\lambda _{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 \(\lambda _{S}\). This can be successfully applied by fixing the length scale \(\lambda _{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

$$\begin{aligned} \lambda _{S} = \sqrt{8 t} = a \sqrt{8 \tau }. \end{aligned}$$
(68)

The above equation enables us to translate the length scale \(\lambda _{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

$$\begin{aligned} \lambda _S = a \sqrt{\frac{8 n_c}{3} }. \end{aligned}$$
(69)

For APE smearing, we can write \(\lambda _S\) as a function of \(n_\mathrm{APE}\) as

$$\begin{aligned} \lambda _S = a \sqrt{\frac{4 \alpha _{\mathrm{APE}} n_{\mathrm{APE}}}{3} }. \end{aligned}$$
(70)

In the case of stout smearing, \(\lambda _S\) takes the form

$$\begin{aligned} \lambda _S = a \sqrt{{8 \rho _{\mathrm{st}}} n_{\mathrm{st}} }. \end{aligned}$$
(71)

Finally, for HYP smearing, we can use the numerically extracted \(\tau _{\mathrm{HYP}} (n_{\mathrm{HYP}})\) in Eq. (67) to define \(\lambda _S\) as a function of \(n_{\mathrm{HYP}}\) according to

$$\begin{aligned} \lambda _S = a \sqrt{{8 \tau _{\mathrm{HYP}} (n_{\mathrm{HYP}}) }}. \end{aligned}$$
(72)

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) of \(\lambda _S =\sqrt{8 t} = {\mathcal {O}}(0.1 \mathrm{fm})\). The value of \(\lambda _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, \(\lambda _S\) corresponds to a plateau for the topological susceptibility reflecting the scale invariance of the observable. We, therefore, extract the topological susceptibility at fixed \(\lambda _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 \tau (n)\) where \(\tau (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 \tau _0 = a^2 \times 2.5\); in other words the dimensionless flow time \(\tau _{0}=2.5\). For each individual smoother, we can extract a smoothing scale which matches this flow time.

Table 3 The matching prefactors between the smoothing schemes of the Wilson flow with time \(\tau \), cooling at level \(n_c\), APE smearing with level \(n_{\mathrm{APE}}\) and finally stout smearing with level \(n_{\mathrm{st}}\). The leftmost column corresponds to the left hand side of the matching equation while the uppermost row to the scale of the right hand side i.e. \(n_{\mathrm{APE}} \simeq \frac{2}{\alpha _{\mathrm{APE}}} n_{c}\)

4 Results

4.1 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 divergences. 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 \(\tau \) and \(n_c\). In the left panel of Fig. 1 we show the function \(\tau (n_c)\) defined as the Wilson flow time \(\tau \) for which the average action density

$$\begin{aligned} \left\langle {\bar{S}}_{G} \right\rangle= & {} 1 - \left\langle \frac{\sum _{x} \sum _{\begin{array}{c} \mu ,\nu =1\\ 1\le \mu <\nu \end{array}}^4 \mathrm{Re} \mathrm{Tr}U^{1\times 1}_{x,\mu ,\nu } (\tau )}{6 Va^{-4}N} \right\rangle \end{aligned}$$
(73)

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 \(\tau =n_c/3\). Furthermore, in the right panel of Fig. 1, we present the average action density as a function of \(\tau \) 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 \(\tau =n_c/3\) leads to equivalent results for the average action density between the Wilson flow and cooling for small values of \(\tau \) and \(n_c\).

Fig. 1
figure 1

Left Panel: The behavior of \(\tau (n_c)\) as a function of \(n_c\) for the Wilson smoothing action. The line corresponds to \(\tau = n_c/3\). Right Panel: The average action density as a function of the Wilson flow time \(\tau \) or the corresponding cooling step \(n_c/3\)

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 \(\alpha _{\mathrm{APE}}\), namely \(\alpha _{\mathrm{APE}}=0.4\), 0.5 and 0.6. Subsequently, we calculated the function \(\tau (\alpha _{\mathrm{APE}}, n_{\mathrm{APE}})\) defined as the Wilson flow time \(\tau \) for which the average action density reduces by the same amount as when \(n_\mathrm{APE}\) smearing steps, for a fixed value of \(\alpha _{\mathrm{APE}}\), are performed. In the left panel of Fig. 2, we demonstrate \(\tau (\alpha _{\mathrm{APE}}, n_{\mathrm{APE}})\) for the three different values of \(\alpha _{\mathrm{APE}}\). The data points for the three values of \(\alpha _{\mathrm{APE}}\) appear to agree with the lines \(\tau =(\alpha _{\mathrm{APE}}/6) \times n_{\mathrm{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 \(\tau \) and \((\alpha _{\mathrm{APE}}/6) \times n_{\mathrm{APE}}\) demonstrating that the four sets of data perfectly agree with each other.

Fig. 2
figure 2

Left Panel: The behavior of \(\tau (\alpha _{\mathrm{APE}}, n_{\mathrm{APE}})\) as a function of \(n_{\mathrm{APE}}\) for \(\alpha _\mathrm{APE}=0.4\), 0.5, 0.6. The red, blue and black lines correspond to \(\tau = (0.4/6) n_{\mathrm{APE}}\), \((0.5/6) n_{\mathrm{APE}}\) and \((0.6/6) n_{\mathrm{APE}}\) respectively. Right Panel: The average action density as a function of the Wilson flow time \(\tau \) or the corresponding rescaled APE smearing step \(\frac{\alpha _{\mathrm{APE}}}{6} n_{\mathrm{APE}}\)

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 equivalence, we smoothed the gauge configurations using stout smearing and three different values of the \(\rho _{\mathrm{st}}\), namely \(\rho _{\mathrm{st}}=0.01\), 0.05 and 0.1. Now we define the function \(\tau (\rho _{\mathrm{st}}, n_{\mathrm{st}})\), like in the previous cases, as the Wilson flow time \(\tau \) for which the average action density alters by the same amount as when \(n_{\mathrm{st}}\) stout smearing steps with a given \(\rho _{\mathrm{st}}\) are performed. In the left panel of Fig. 3, we present the function \(\tau (\rho _{\mathrm{st}}, n_{\mathrm{st}})\) for the three different values of the parameter \(\rho _{\mathrm{st}}\). The data points for the three values of \(\rho _{\mathrm{st}}\) appear to agree with the lines \(\tau =\rho _\mathrm{st} \times n_{\mathrm{st}}\) providing strong evidence that Eq. (66) provides the right rescaling for which Wilson flow and stout smearing become numerically equivalent. This agreement sets in for \(\tau \simeq 5\), 2, 1 for \(\rho _{\mathrm{st}} = 0.1, \ 0.05, \ 0.01\), respectively, demonstrating that the smaller the value of \(\rho _{\mathrm{st}}\) the closer we approach the Wilson flow. Furthermore, in the right panel of Fig. 3, we provide the average action density as a function of \(\tau \) or \(\rho _{\mathrm{st}} \times n_{\mathrm{st}}\) for the Wilson flow or stout smearing, respectively, demonstrating that the four sets of data perfectly agree with each other. Similar comparisons of the topological charge and susceptibility are provided in Section 4.2.

Fig. 3
figure 3

Left Panel: The behavior of \(\tau (\rho _{st}, n_\mathrm{st})\) as a function of \(n_{\mathrm{st}}\) for \(\rho =0.01\), 0.05, 0.1. The red, blue and black lines corresponds to \(\tau = 0.01 n_\mathrm{st}\), \(0.05 n_{\mathrm{st}}\) and \(0.1 n_{\mathrm{st}}\), respectively. Right Panel: The average action density as a function of the Wilson flow time \(\tau \) or the corresponding rescaled stout smearing step \(\rho _{\mathrm{st}} \times n_{\mathrm{st}}\)

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 perturbative rescaling between the Wilson flow time \(\tau \) and the number of HYP smearing steps \(n_{\mathrm{HYP}}\). Thus, instead, we attempted a numerical fit using a parametrization of \(\tau (n_\mathrm{HYP})\) in \(n_{\mathrm{HYP}}\) according to Eq. (67). In the left panel of Fig. 4, we provide the function \(\tau (n_{\mathrm{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)\times 10^{-5}\). Of course, these numbers depend on the parameters \(\alpha _{\mathrm{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).

Fig. 4
figure 4

Left Panel: The behavior of \(\tau (n_{\mathrm{HYP}})\) as a function of \(n_{\mathrm{HYP}}\). The green line corresponds to a linear approximation which is valid up \(t_0\) while the blue line corresponds to the numerical fit \(\tau (n_{\mathrm{HYP}})\). Right Panel: The average action density as a function of the Wilson flow time \(\tau \) and the corresponding numerical matching \(\tau ({n_\mathrm{HYP}})\)

4.2 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 \(\tau \) and for matched smoothing scales for cooling, APE, stout as well as HYP smearing has been investigated. In Fig. 5, we present the clover definition of the topological charge as a function of \(\tau \), \(n_c/3\), \(\alpha _{\mathrm{APE}} n_{\mathrm{APE}}/6\) for \(\alpha _\mathrm{APE}=0.5\), \(\rho _{\mathrm{st}} n_{\mathrm{st}}\) for \(\rho _{\mathrm{st}}=0.05\) and \(\tau _{\mathrm{HYP}}(n_{\mathrm{HYP}})\), for four randomly chosen configurations in the ensemble b40.16; each panel corresponds to each configuration.

Fig. 5
figure 5

For four different gauge field configurations, we show the clover definition of the topological charge as a function of the gradient flow time \(\tau \) for Wilson flow (line in red), \(n_c/3\) for cooling (), \(0.5 \times n_{\mathrm{APE}}/6\) for APE smearing (), \(0.05 \times n_{\mathrm{stout}}\) for stout smearing () and \(\tau _{\mathrm{HYP}} (n_{\mathrm{HYP}}))\) for HYP smearing (). For this ensemble \(t_0\simeq 2.5a^2\)

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_{\mathrm{APE}}\) and \(n_{\mathrm{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 \(\tau =6-8\), the approximate plateau shifts from \(Q\simeq 2\) to \(Q \simeq 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 \(\rho _{\mathrm{st}}\) and \(\alpha _{\mathrm{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 the topological charge is not necessarily the same; this is demonstrated also in Section 4.6 of this manuscript.

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_{\mathrm{HYP}}\) with the numerically-extracted 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 \(\tau \); however, for this range of \(\tau \) 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.

4.3 Monte Carlo histories and distribution histograms

In this subsection, we show some typical features of the topological charge evaluated with one representative fermionic definitionFootnote 6 (index of the overlap operator evaluated on configurations with 1 step of HYP smearing applied) and one representative gluonic definition (with the Wilson flow at flow time \(t_0\)). This serves two purposes. First, we show that correlations between the two classes of definitions are apparent even from a visual inspection of Monte Carlo histories. Second, we want to investigate whether the distribution of the topological charge is approximately Gaussian and, for the gluonic case, whether clustering of the values around integers occurs.

Fig. 6
figure 6

Monte Carlo histories for one representative fermionic definition (index of the overlap operator evaluated on configurations with 1 step of HYP smearing applied) and one representative gluonic definition (Wilson flow at flow time \(t_0\))

The Monte Carlo histories are shown in 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 \(\tau _{\mathrm{int}}\in [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 \(\tau _{\mathrm{int}}\approx 1.7(3)\). The correlation between the values of Q from the two definitions in Fig. 6 is obvious even without computing 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.

Fig. 7
figure 7

Scatterplot of the topological charge for one representative fermionic definition (index of the overlap operator evaluated on configurations with 1 step of HYP smearing applied) and one representative gluonic definition (Wilson flow at flow time \(t_0\)). For better visibility, the integer values of the index were randomly shifted by a small non-integer value

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 GaussianFootnote 7. For the field theoretic definition, we used an interval width of 0.1 to detect clustering around values close to an integer.

Fig. 8
figure 8

Histograms of the topological charge for ensemble b40.16. The employed definitions are: index of the overlap operator evaluated on configurations with 1 step of HYP smearing applied (top left) and the Wilson flow at flow times \(t_0\) (top right), \(2t_0\) (bottom left) and \(3t_0\) (bottom right)

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 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.

4.4 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.

4.4.1 Main comparison

In this subsection, we choose the following definitions as typical representatives:

  • index of the overlap Dirac operator applied to non-smeared 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).

Fig. 9
figure 9

Main comparison of selected definitions of the topological charge. The correlation between different definitions is colour-coded (note the scale is different than in Figs. 10151614)

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)

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 non-smeared 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 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 \({{{\mathcal {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 actions – while 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.

4.4.2 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.

Fig. 10
figure 10

Comparison of fermionic definitions of the topological charge. The correlation between different definitions is colour-coded (note the scale is different than in Figs. 9151614)

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)

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 \(\pm 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 (\(\overline{\text {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 \(aM\approx 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\approx 100\) MeV, the number of eigenmodes of the operator \(D^\dagger 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 \({\mathbb {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 charges extracted from different M values are very small – the ones with the smallest M are only 25%-30% correlated with the ones at higher values of M. The correlations among higher values of M (\(M^2\ge 0.0004\)) are somewhat higher (45–60%), but still only moderate. This suggests that the stochastic noise is not entirely suppressed and that cut-off effects are possibly very different for different values of M.

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 \({{{\mathcal {C}}}}\) as the topological charge. Although the topological susceptibility is well-defined with spectral projectors, it requires to correct \(\langle \mathcal{C}^2\rangle \) with \(\langle {{{\mathcal {B}}}}\rangle /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 \({{{\mathcal {C}}}}\)), this makes the interpretation of \({{{\mathcal {C}}}}\) on a single gauge field configuration difficult. It is, moreover, likely that the values of \({{{\mathcal {C}}}}\) extracted with different values of M are affected by different cut-off effects.

4.5 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.

Fig. 11
figure 11

Increase of correlation towards the continuum limit between 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\)). The dashed red line is a guide to the eye, showing that it is plausible that the correlation will become 1 in the continuum limit

4.6 Topological susceptibility

Fig. 12
figure 12

Topological susceptibility for ensemble b40.16 (with pion mass around 340 MeV), using representatives of different kinds of definitions. The definition numbers corresponds to the ones in Table 2

We also compared the topological susceptibility computed with different methods. We defined the topological susceptibility as

$$\begin{aligned} \chi =\frac{\langle Q^2\rangle }{V}, \end{aligned}$$
(74)

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).

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\chi ^{1/4}\in [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 \(\chi \) 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 \(\chi \). 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 \(\chi \) 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.

Fig. 13
figure 13

A comparison of the topological suscpeptibility \(r_0 \chi ^{1/4}\) as a function of the gradient flow time \(\tau \) for the Wilson flow, or, in the upper left panel, the rescaled cooling step \(n_c/3\) for cooling, or, in the upper right panel, the rescaled HYP smearing \(\tau _{\mathrm{HYP}}(n_{\mathrm{HYP}})\), or, in the lower left panel, the rescaled \(n_{\mathrm{stout}}\) smearing \(\rho _{\mathrm{st}} \times n_{\mathrm{stout}}\) for three values of \(\rho _{\mathrm{st}}\), or, finally in the lower right panel, the rescaled APE smearing \(\alpha _{\mathrm{APE}} \times n_{\mathrm{APE}} /6\) for six different values of the \(\alpha _\mathrm{APE}\) parameter. For this ensemble \(t_0\simeq 2.5a^2\)

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 \(\chi \). 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 \chi ^{1/4}\) extracted via the Wilson flow and cooling. By rescaling \(n_c\) \(\rightarrow \) \(n_c/3\) according to Eq. (49) we observe that both smoothers give results which agree for the whole range of \(\tau \). 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 \(\tau \sim 3\). On the contrary, when smoothing with HYP smearing, the plateau sets in at larger values of \(\tau \) 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 \chi ^{1/4}\) obtained via the Wilson flow as a function of the flow time with \(r_0 \chi ^{1/4}\) calculated via stout smearing as a function of the rescaled smearing steps of \(\rho _{\mathrm{st}} n_{\mathrm{st}}\). We consider three values of \(\rho _\mathrm{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 \chi ^{1/4}\) measured with the Wilson flow as a function of \(\tau \) with the value of \(r_0 \chi ^{1/4}\) extracted with APE smearing vs. the rescaled smearing step \(\alpha _{\mathrm{APE}} n_{\mathrm{APE}}/6\). We did so for six different values of \(\alpha _{\mathrm{APE}}\), namely \(\alpha _{\mathrm{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 \(\alpha _\mathrm{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).

5 Conclusions

In this paper, we have investigated several definitions of the topological charge. Our main conclusion is that all validFootnote 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 fasterFootnote 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.