Topological susceptibility and the sampling of field space in $N_f=2$ lattice QCD simulations

We present a measurement of the topological susceptibility in two flavor QCD. In this observable, large autocorrelations are present and also sizable cutoff effects have to be faced in the continuum extrapolation. Within the statistical accuracy of the computation, the result agrees with the expectation from leading order chiral perturbation theory.


Introduction
The topology of the gauge fields plays an important role in the understanding of the low energy structure of QCD. A prominent example is the Witten-Veneziano relation where the topological susceptibility of pure Yang-Mills theory is linked to the mass of the η meson [1,2]. The object of this paper is the topological susceptibility χ top in two flavor QCD, a quantity which vanishes in the chiral limit but whose value is non-zero for finite sea quark masses.
Precise measurements of χ top (m π ) are a challenge for the lattice. In particular significant effort went into finding a suitable definition of the topological charge on a discrete space time. Also large autocorrelations in the numerical simulations and sizable cutoff effects make the determination difficult.
These issues have a long history which we can only sketch briefly. The naive "field theoretical" definition of the topological charge density using some discretization of the field strength tensor F is known to lead to a topological susceptibility χ top which suffers from non-integrable short-distance singularities. This has led to a vast array of operational procedures based on cooling and link smearing which remove the ultra-violet fluctuations of the gauge fields, but whose behavior towards the continuum limit is not well understood. However, a number of practical definitions with a well-defined continuum limit are now available, notably based on the properties of chiral Dirac operators and the index theorem [3] or ratios of certain fermionic correlation functions [4][5][6].
In this study, we adopt the definition of the topological charge through the Wilson flow [7], which provides a smoothing of the gauge fields with a smearing radius that is kept constant in physical units as the continuum is approached. Numerical support for a common continuum limit of this definition of the susceptibility and the method of Refs. [4][5][6] in pure gauge theory has been given in Ref. [8].
While the computation of the topological charge using the Wilson flow requires only a moderate numerical effort, the simulations are still expensive due to large autocorrelation times in the topological charge. In particular if periodic boundary conditions are adopted, the global topological charge suffers from a severe critical slowing down [9] and a determination of the susceptibility becomes practically impossible. The particularly poor scaling of the autocorrelation times due to the topological charge can be solved by using open boundary conditions in time [10,11], however, the smooth fields in the construction of the topological charge are still slow in the Monte Carlo evolution.
Studying the effect of dynamical fermions on the susceptibility using lattice QCD has a long tradition, but even a decade ago, the situation was less than clear. While the Sesam collaboration did not find a convincing suppression of the susceptibility [12], the UKQCD collaboration did observe some effect [13]. Later, the importance of cut-off effects in this observable was discussed for simulations with Wilson fermions [14] and staggered quarks [15,16]. After taking the continuum limit, a behavior compatible with the continuum expectation from Chiral Perturbation Theory [17] was found. More recently, the suppression of the susceptibility with the quark mass in two flavor QCD was also demonstrated in Refs. [18] and [19] using Wilson fermions with a twisted mass and domain wall fermions, respectively. 1 Due to the limitations by the statistical error and the discretization effects, the scope of this study is twofold: On the one hand, we want to present the measurement of χ top (m π ) as well as it can be made with currently available ensembles and present evidence for the proper suppression of the topological charge fluctuations with the quark mass. On the other hand, studying the autocorrelations in a set of gauge field configurations which are used in many computations serves as a crucial input in the error analysis of these computations [9]. We will also have a first look at the impact of open boundary conditions in time in simulations including dynamical fermions.
The paper is therefore organized as follows: In Sect. 2 the observables are defined and an overview of the ensembles used in the current study is given. We examine the autocorrelations in observables constructed from smoothed gauge fields and give new estimates of the largest exponential autocorrelation times. In Sect. 3, after the separation of the gauge fields into topological sectors is presented for the two flavor theory, we then turn to the computation of the susceptibility, the discussion of its discretization effects and a comparison of the continuum limit of χ top (m π ) with Chiral Perturbation Theory.

Setup of the calculation
The Wilson flow in the space of lattice gauge fields U (x, µ) is defined by the equations where the action S W is the standard Wilson action and V t are the smoothed link variables at flow-time t. It was introduced in Ref. [7] and we follow the conventions adopted there. In particular we will look at the topological charge density q(x, t) with a clover-type discretization of the field strength tensor G µν (x) constructed from links V t . Furthermore the energy density constructed from these links will be considered The quantities defined through the flow need to be evaluated at a fixed, physical value of t. Since it provides a smoothing of the fields over a radius r = √ 8t, it is natural to use the reference scale t 0 introduced in the original paper [7] through Recently it has been argued that cutoff effects in some quantities can be reduced by considering the flow at a larger time t 1 , where the right hand side of eq. (2.4) is replaced by 2/3 [20]. We will give results in units of t 0 and t 1 . To convert to physical units the reader may use t 0 = 0.0236 fm 2 [21] and t 1 = 0.061 fm 2 , both obtained for N f = 2 from fixing the kaon decay constant and the pion mass, m π , to experiment. Their errors are negligible for our present purposes. While these numbers refer to the physical point, in the rest of this paper we always use the scales t 0 (m π ), t 1 (m π ) as obtained at the simulated values of m π . The topological susceptibility derived from the charge density eq. (2.2) does not require renormalization. While its continuum limit does not depend on the flow time t, discretization effects will be influenced by its choice. However, we will see below that choosing either t = t 0 or t = t 1 does not have a sizable effect on a 4 χ top .
Starting from the topological charge density defined in eq. (2.2), the topological charge can be directly obtained (2.5) In periodic boundary conditions, the susceptibility is then simply χ top = Q 2 /V .

Ensembles
The set of ensembles used for our study was generated within the CLS effort. 2 Two degenerate light quarks have been simulated using O(a)-improved Wilson fermions and Wilson gauge action. The ensembles are listed in Table 1: three lattice spacings and a range of pion masses are available with all lattices fulfilling m π L > 4.  Table 1: Overview of the ensembles used. We give the parameters of the action β and κ, the temporal and spatial extent of the lattice and the pion mass. For the simulations done with the DD-HMC algorithm, the ratio of active links, R, with which the τ int are scaled throughout the paper, is below one. The last two columns give the total statistics in molecular dynamics units (MDU) and how this compares to the τ exp estimated in eq. (2.9).
The majority of the ensembles was generated with the DD-HMC algorithm [22], while the HMC algorithm with Hasenbusch preconditioning [23] has been employed for those with smaller pion masses [24]. The ensemble where open boundary conditions in time have been imposed has been generated with the openQCD code [25]. The details of most of the simulations can be found in Ref. [26].
Since a significant part of this text deals with autocorrelations, it should be noted that algorithms and their parameters influence them. The different ways the fermion determinants are treated in the three algorithmic setups will certainly matter once accuracies are high. At the current level of precision, we only correct for the inactive links during the trajectories of the DD-HMC algorithm by scaling the molecular dynamics time with R, the fraction of active links. In the high statistics pure gauge theory study of Ref. [9] this has proven to correctly compensate for their effect. We will also scale autocorrelation times with the acceptance rate P acc in the final scaling formula.

Autocorrelations
Large autocorrelations are a well-known problem in lattice simulations. They are particularly pronounced in quantities constructed from the fields smoothed by the Wilson flow. In these quantities, a critical slowing down governed by a dynamical critical exponent of z = 2 is expected [27]. To compare simulations at different lattice spacings, we will therefore scale the Monte Carlo time with the dimension two quantity t 0 measured on the respective ensemble.
An additional problem is posed by the freezing of the topological charge as the continuum limit is approached. This occurs with a larger critical exponent -if it is not exponential -and makes simulation below a ≈ 0.05 fm exceedingly difficult. The problem is not restricted to the topological charge alone. In general, the associated slow eigenmodes of the Markov chain transition matrix contribute to the autocorrelation functions of all observables. For a correct estimate of the error a careful study of the slowest modes of the Markov chain is mandatory.
The basis of the error analysis of Markov Chain Monte Carlo data is the measurement of the normalized autocorrelation function ρ A with A(τ ) the Monte Carlo time history of the observable A and τ the Monte Carlo time.
From the analysis of the spectral decomposition a recipe to include the coupling to the slow modes of the transition matrix in the error analysis has been given in Ref. [9]: the normalized autocorrelation function ρ is summed as usual up to a window W and then the contribution of a single exponential with time constant τ exp is attached For this method an estimate of the largest time constant τ exp in the simulation is needed. Since the dynamics of parity even and parity odd observables decouples  Table 1. Points have been shifted for better legibility. ρ(τ ) is summed up to a window W = 100 MDU, neglecting a potential tail. The lines are fits of eq. (2.8) to the data with with c 0 ≈ 58 MDU for Q(t) 2 and c 0 ≈ 69 MDU for E(t).
[9], we can restrict ourselves to even ones. Then Q 2 (t) is an obvious candidate to look for the largest time constant, but we also investigate E(t).
For the autocorrelation analysis we have to pick a value of the flow time t. In Fig. 1 we present the dependence of τ int on this parameter and find a behavior similar to the one observed in pure gauge theory [25]: the data can be described by with the value at t = t 0 essentially capturing the asymptotic value. We therefore use this value of the flow time and drop the flow time argument in the following, using Q = Q(t 0 ) and E = E(t 0 ).

Results for autocorrelation times
Autocorrelation functions depend on the observable, the lattice spacing and the quark mass. To quantify the dependence on the physics parameters, we give examples of autocorrelation functions for two values of the quark mass in Fig. 2.
The data at β = 5.5 shows a suppression with the quark mass of ρ Q 2 , as also found in Ref. [28], whereas ρ E is not affected beyond statistics. This can be understood from the fact that the topological susceptibility goes to zero in the  Figure 3 shows the dependence of the autocorrelation functions on the lattice spacing. While the rescaling of the Monte Carlo time with t 0 confirms the scaling properties of ρ E already observed in pure gauge theory, the topological charge decorrelates much faster on the coarser lattices. On our finest lattice spacing, the autocorrelation functions of E and Q 2 are on top of each other and we expect the topological charge to dominate on finer lattices. This qualitative behavior is again similar to the one in pure gauge theory with periodic boundary conditions [10], only that it happens at smaller lattice spacing -not a surprise in light of the suppression of autocorrelations in Q 2 with the quark mass.
The results of the integrated autocorrelation times (summed up to Rτ ≈ 40t 0 /a 2 ) for Q 2 and E are given in Table 2. We also tried single exponential fits to the various ρ E in order to find τ exp , the largest time constant of the Monte Carlo chain. It turns out that the ρ(τ ) are well described by this ansatz over essentially the full range of τ .
Since E has the largest autocorrelations in our range of lattice spacings and its autocorrelation function is approximately an exponential function, we can identify τ int (E) with τ exp at the current level of accuracy. This leads to a parametrization  In an earlier publication [9], an estimate of τ exp (β) has been given based on the scaling of τ int (Q 2 ) in pure gauge theory normalized to the value of the E5 lattice. Since the topological charge is not the slowest quantity in this region, this is overly pessimistic for β = 5.5, where it gives twice the τ int (Q 2 ) observed and overly optimistic at β = 5.2.
Since the autocorrelations of the topological charge show a scaling compatible with what has been found in pure gauge theory, i.e. τ int (Q 2 ) ∝ a −5 , we expect that for our range of quark masses and β > 5.5 the charge will dominate the scaling of τ exp .

Open boundary conditions
The use of open boundary conditions is expected to accelerate the sampling of the topological charge, which at fine lattice spacings is expected to be the slowest observable. However, we observe that this is not the case for our lattices, in particular at β = 5.2 where oB6 has been simulated. Here the bulk tunneling dominates over the effect of the open boundaries.
Within the accuracy of the result for τ int (E) reported in  With open BC the notion of a global topological charge is lost, as will be discussed in the next section, therefore we do not compute τ int (Q 2 ) for oB6. For the alternative way to extract the susceptibility discussed below, we find comparable autocorrelation times for the two ensembles.

Results
From the discussion in the previous section it should be clear that while we are confident that the statistics is sufficient for having control over the statistical errors, the accuracy of the susceptibility which can be reached from these ensembles only allows for qualitative statements. It is certainly not on par with the precision reached in pure gauge theory.
We will start by discussing the separation of topological sectors as the continuum is approached. In some sense, this is a measure for the physics associated with topology being realized already at a finite lattice spacing. It also is the reason for the Hybrid Monte Carlo algorithm having difficulty to decorrelate the topological charge.
Once this is established, we turn to the topological susceptibility, studying systematic effects between different ways of constructing it. This will be important for its determination imposing open boundary conditions in time, where the global topological charge is no longer a suitable observable.

Separation between the sectors
With periodic boundary conditions field space is disconnected in the continuum gauge theory. This property is obfuscated by ultra-violet fluctuations of the gauge fields and for a long time whether and in which way these sectors emerge as the lattice spacing decreases has been unclear. Using the smoothing provided by the Wilson flow, the emergence of the topological sectors in the continuum limit can finally be understood [10]. For convenience we here repeat the arguments in [10] in a rather explicit manner.
First one notices that it was shown long ago, that the space of lattice gauge fields U separates into disconnected sectors provided the gauge fields U are smooth enough, which means that for all plaquettes p on the lattice [29,30]. 3 Secondly, the Wilson flow defines a map U ≡ V 0 ↔ V t of the gauge fields in our path integral to smoothed fields V t , t > 0. Since this map is one-to-one, a classification of V t in sectors is also one of the original fields U . Fields V t are smooth since the total (global) action S W (V ) = p 1 g 2 0 s p (V ) is decreased by the flow, S W (V t ) < S W (V t ) for t > t . Of course, this does not mean that the local bound eq. (3.2) will be satisfied for V t at a given t. However, if the number of plaquettes which violate s p (V t ) < s cut decreases strongly as one approaches a → 0, we can say that sectors develop dynamically in the continuum limit.
Thirdly, in Ref. [10] violations of the bound eq. (3.2) are investigated numerically using the probability that one plaquette violates the bound 3) It finds that in pure gauge theory with the Wilson action P t (s cut ) ∼ a 10 , which also suggests that in a fixed volume the probability for any plaquette violating the bound decreases as a 6 . As can be seen in Figure 4, we find the same power law behavior in our two flavor simulations. While two values of s cut are shown, the scaling does not depend significantly on s cut in the investigated range up to s cut = 0.1 and the exponent is compatible with the pure gauge theory result of Ref. [10]. It is remarkable that the effect of the quark mass on P t 0 (s cut ) is completely compensated by its effect on t 0 and we observe universal scaling curves for each value of s cut . At a fixed β, smaller quark masses lead to a stronger separation of the sectors. In particular the presence of fermions does not alter the picture of Ref. [10].

Computation of the topological susceptibility
With periodic boundary conditions for the gauge fields, is the natural estimator for the topological susceptibility. With open boundary conditions in time, the same quantity would show large finite volume effects. 4 Therefore we rewrite eq. (1.2) restricting the time separation between the two densities to an upper bound r with χ corr (T /2 − a) + a 2 T L 3 q(0)q(T /2) giving the full formula for periodic boundaries with even T /a. The correlator q(z 0 )q(x 0 + z 0 ) will be studied below. As we will see, it falls off quickly for large x 0 such that the sum can be truncated at moderate r.
To avoid finite T effects, all temporal arguments need to be sufficiently far away from the boundaries, i.e., x 0 , T − x 0 , x 0 + z 0 , T − (x 0 + z 0 ) need to be large in units of the inverse mass of the lightest state with vacuum quantum numbers.
In order to study the effect of the truncation in eq. (3.5), we plot the correlator between the slice charges for several values of the lattice spacings on the periodic lattices in the left panel of Fig. 5. This correlation function is positive for small distances, then turns negative and approaches zero for x 0 → ∞. The sum over x 0 as a function of the upper bound is depicted on the right panel. As we can see, for our statistics a plateau is reached around r/ √ t 1 = 5. This holds for all our ensembles. We therefore cut the summation at this point; the results are collected in Table 3.
In these figures, it is interesting to note that while the large x 0 region shows almost perfect scaling, a significant cut-off effect is visible at x 0 = 0. Indeed, if this point is excluded from the summation in eq. (3.5), the sum does not show any scaling violation beyond its statistical accuracy.
In Table 3 we also report the result of our ensemble with open BC, namely oB6. We observe that for this lattice spacing and pion mass the influence of the boundaries in q(x 0 ) 2 is visible with the present accuray up to 1.5 fm. Therefore, we restrict the computation of the susceptibility given in eq. (3.5) to time slices in the range [20,171]; the summation window remains fixed to r = 5 √ t 1 . No differences to B6 are observed in the final result and the errors are about the same. The reason is that the statistics is roughly the same. Furthermore, the reweighting factors of the twisted mass reweighting employed in the oB6 ensemble decouple almost completely from pure gluonic observables such as the charge.

Numerical estimate of the tail
We then fitted the correlation function to a single exponential as shown in Fig. 5, and estimated the systematic error of χ corr (r) from the integral over the fit func-  Figure 5: Determination of the susceptibility from the correlation function of the time slice summed charge at t = t 0 (open symbols) and t = t 0 /2 (filled symbols); the scales t 0 and t 1 are taken at the quark mass of the ensemble. In both plots the ensembles are A5, F6 and N6. For the time-slice correlation functions in the left panel, the long distance behavior can be described by a single exponential decay as discussed in the text. In the right panel, the sum up to a window r is represented, with the vertical line the one used in our analysis.  Results of the topological susceptibility from eq. (3.4), with Q measured also at t = t 1 , and eq. (3.5). The scales t 0 and t 1 do always take the value at the finite quark mass of the corresponding ensemble. All the values have been multiplied by a factor 10 3 . For most of the ensembles, χ corr turns out to be more precise.
tion. As argued below this ansatz is justified in this situation despite the fourdimensional smoothing of the flow. The fit was done to the data with t = t 0 /2 shown in the figure, namely together to the data at the three values of a, neglecting scaling violations and in a range 2.5 ≤ x 0 t −1/2 1 ≤ 6. The thus determined systematic error turns out to be below 1%, much smaller than our statistical one. Even if our model underestimated the integral by a factor 4, which seems unlikely given the reasonable agreement with the data, the tail would still be negligible with our present accuracy.

Effects of the smoothing footprint in the large time behavior
The neglected tail gives a systematic error of this measurement. In order to justify the above estimate of its contribution, we need to be sure that the single exponential ansatz is expected to give a satisfactory description of the data in the range where it is applied.
The exponential decay with the mass of the flavor singlet pseudoscalar is modified due to two effects: the contribution of excited states and the fact that we apply a four dimensional smoothing to our observables. Since the flow generates a footprint of size √ 8t with a Gaussian profile, one expects that for sufficiently large x 0 this modification is negligible. With a lowest mass m 0 in the considered channel, dimensional analysis leads one to expect that it can be neglected once x 0 /(8t m 0 ) 1. To gain more insight into the influence of the smoothing on the shape of the two-point function, we consider the model of a free smoothed scalar field, noting that in lowest order of perturbation theory, the gluon field is smoothed in exactly the same way. Its time-slice correlation function at zero spatial momentum is given by where m is the the mass of the scalar field. The large time asymptotics is given by where we used dimensionless variables µ = √ 8t m , ξ = x 0 / √ 8t. The corrections to the asymptotic exponential decay due to the footprint of the smoothing are of order ξ −2 e −(ξ−µ/2) 2 . They are small once ξ − µ/2 > c, or equivalently x 0 > c √ 8t + 4tm, and c = 1.5 or c = 3, with the effect of the smoothing in the per cent and far below the per mille range, respectively.
While our explicit computation is in a model, we believe that it provides a good general guideline. In particular the structure with the two different terms is expected to be a general feature.
For our case, the relevant mass is m 0 ≈ 1 GeV. Choosing t = t 0 /2, we have µ ≈ 1.5 and ξ > 2 translates into a footprint effect of less than 3%. Since the discussion in Sect. 3.2.1 is about a systematic effect on the neglected tail, i.e., on a small systematic error itself, such an accuracy is more than sufficient in particular in view of our statistical uncertainties.

Distribution of the topological charge
In the large volume limit, the topological charge is expected to follow a Gaussian distribution, with the width given by the topological susceptibility. In Fig. 6 we show the cumulative distribution function C(Q), which gives the probability of a configuration having charge smaller than Q, along with the measurement on three different ensembles. The theoretical lines use as width of the Gaussian distribution the measured value of χ corr of Table 3. The agreement with a Gaussian distribution is very reasonable, with only small deviations visible in the tails at a = 0.048 fm, which we attribute to our limited sampling of the tails.

Continuum limit and dependence on the light quark mass
In leading order in Chiral Perturbation Theory, the topological susceptibility is linear in the quark mass [17], or using the Gell-Mann-Oakes-Renner relation We therefore display the results for χ corr from Table 3 in Fig. 7 in terms of the dimensionless variable t 1 m 2 π , where we always take t 1 = t 1 (m π ), the value of the respective ensemble. With the lattice spacing determined from the kaon decay constant [26] and the pion decay constant set to 130.4 MeV, we obtain the line represented in the figure. Also shown is the pure gauge result of Ref. [34]. The comparison confirms the strong effect of the presence of the fermions found in previous studies. 5 In the figure, a strong cut-off effect is visible, with the topological susceptibility significantly enhanced for the coarser lattices. Because of the breaking of  Table 3 along with the leading order expectation from ChPT and the pure gauge value. chiral symmetry, at finite lattice spacing there is no reason for the susceptibility to vanish at vanishing pion mass. Assuming leading effects of O(a 2 ), Wilson Chiral Perturbation Theory suggests the following ansatz 6 t 2 1 χ top = c t 1 m 2 π + a 2 b t 1 . (3.10) The result is also shown in Fig. 7. As can be seen, the description of the data is satisfactory given the large error bands and the continuum limit is compatible with the LO ChPT expectation. We find for the slope c = 2.8(5) · 10 −3 , for the intercept b = 5.1(7) · 10 −3 and a χ 2 /d.o.f of the global fit of 1.17.

Conclusions
Measurements of observables with long autocorrelation times are difficult and the topological susceptibility is no exception. However, we found that for the action used in the CLS simulation, the topological charge is not yet the slowest observable in the range of lattice spacings under investigation, though the rapid growth in τ int indicates that on slightly finer lattices the charge would effectively freeze in typical simulations. We show that light quarks mitigate the problem of the slow topological charge to a certain extent, whereas the other slow modes, as seen, e.g., in E(t 0 ), do not profit.
Larger quark masses and also cut-off effects (for the action employed in this paper) both increase the susceptibility, however, even on our coarsest lattices and largest pion mass of 630 MeV, we observe a suppression of the susceptibility with respect to the pure gauge case by more than 50%. The three lattice spacings with a range of pion masses at our disposal allow us to take the continuum limit with theoretical support from Wilson Chiral Perturbation Theory. As we have shown, for fitting the pion mass dependence of χ top at finite lattice spacing, it is crucial to include a term that does not vanish at zero pion mass and finite cutoff. Only with chirally symmetric sea quarks such a term can be excluded.
We stress that discretization effects seem large, because the quantity of interest is small for small quark masses. In such a situation, large relative discretization errors are generic, however, the scale to judge them is the deviation from the pure gauge theory result and on this scale we reach a percent level accuracy in the continuum for the whole quark mass region.
Finally let us point out that the discussion in Sect. 3.2.2 on the effects of the footprint is particularly relevant when one wants to extract masses from correlation functions at positive flow time. For example the glueball mass determination in [35] may be affected.