Topological Susceptibility under Gradient Flow

We study the impact of the Gradient Flow on the topology in various models of lattice field theory. The topological susceptibility $\chi_{\rm t}$ is measured directly, and by the slab method, which is based on the topological content of sub-volumes ("slabs") and estimates $\chi_{\rm t}$ even when the system remains trapped in a fixed topological sector. The results obtained by both methods are essentially consistent, but the impact of the Gradient Flow on the characteristic quantity of the slab method seems to be different in 2-flavour QCD and in the 2d O(3) model. In the latter model, we further address the question whether or not the Gradient Flow leads to a finite continuum limit of the topological susceptibility (rescaled by the correlation length squared, $\xi^{2}$). This ongoing study is based on direct measurements of $\chi_{\rm t}$ in $L \times L$ lattices, at $L/\xi \simeq 6$.


Introduction
In some quantum field theories, the set of configurations is divided into topological sectors, labelled by a topological charge Q ∈ Z. This is the case in QCD, and in N-dimensional O(N + 1) models (with periodic boundary conditions for the gluon and spin fields), due to For usual lattice actions, all configurations can be continuously deformed into one another, at finite action, hence there are no topological sectors in a strict sense. Exceptions are topological lattice actions, with a sharp cutoff for the angles between nearest neighbour spin variables [1], or for each plaquette variable [2] (see also Ref. [3]) in spin models and gauge theories, respectively. However, even for conventional lattice actions there are established ways to divide the configurations into sectors, which turn into topological sectors in the continuum limit.
Here we consider 2-flavour QCD with twisted-mass quarks [4] (at full twist) and the Wilson gauge action. For the O(N) models we employ the standard lattice action, S [ e ] = β xy (1 − e x · e y ) , e x ∈ S N−1 ∀x , N = 2 or 3 , where the sum runs over all nearest neighbour lattice sites. For these O(N) models we apply the geometric definition of the topological charge density on the lattice [5], which leads to integer charges Q ∈ Z. In QCD we use a clover discretisation of F µνFµν , where F is the field strength tensor. In all cases under consideration, parity symmetry implies Q = 0, hence the topological susceptibility takes the form 2 The slab method to measure the topological susceptibility χ t Once we have fixed a formulation of the topological charge on the lattice, it is straightforward to measure χ t by means of Monte Carlo simulations, if the Markov chain frequently changes Q, such that the sectors are sampled correctly. In practice, however, such simulations are often confronted with the severe problem of "topological freezing": in particular, the algorithms, which proceed in small update steps, tend to get stuck in one topological sector for a huge number of steps, since the topological sectors are effectively separated by high potential barriers. The autocorrelation time with respect to Q increases with a high power of the inverse lattice spacing as we approach the continuum limit ("topological slowing down"), see e.g. Ref. [6]. A variety of approaches to handle this problem is reviewed in Ref. [7]. One strategy aims at extracting physical observables even from a Markov chain which is entirely trapped in a single topological sector. For general observables such a method was suggested in Ref. [8], and tested and extended in Refs. [9][10][11][12]. More specifically, a procedure to measure χ t within a fixed topological sector was proposed in Ref. [13] and tested in Refs. [9,12,14,15]. Here we consider the slab method as another way to evaluate χ t from data obtained at fixed Q (actually data from ±Q can be combined). The idea was mentioned in Ref. [16], implemented in Ref. [17], and further explored in Refs. [12,18,19]. A different variant was applied in Ref. [20], and there are similarities with the approach in Ref. [21].
We briefly review the simplest version of the slab method, which assumes the statistical distribution of the topological charges to be Gaussian [17], p(Q) ∝ exp(−Q 2 /(2χ t V)). We split the volume V into two sub-volumes ("slabs") xV and (1 − x)V (0 < x < 1). By summing up the topological charge density in each of them, in a configuration of topological charge Q, we obtain the slab charges q, Q − q ∈ R (they do not need to be integer, since the slabs do not have periodic boundaries). At fixed x, V and Q, the corresponding slab probability distributions p 1 and p 2 obey Measuring q 2 yields a value for q 2 = q 2 − x 2 Q 2 . A sequence of such measurements, at different parameters x, enables a fit to the prediction which provides a result for χ t . In practice, the most reliable fitting regime is around the center, x ≈ 0.5, because the size of the slabs should be large compared to that of topological excitations; one should not include x > ∼ 0 or x < ∼ 1 (these regions involve very small slabs, where the Gaussian distribution is not a good approximation).

Results by the slab method under Gradient Flow (GF)
As a smoothing procedure for lattice field configurations, the GF corresponds to a renormalisation group scheme. When the GF proceeds, the distinction between topological sectors becomes more marked, approaching the continuum feature of a total separation [22,23].
The slab method has been tested in 2-flavour QCD, with twisted-mass quarks and the Wilson gauge action, in a volume 16 3 × 32, at β = 3.9 and bare mass 0.015, which corresponds to a pion mass of m π 650 MeV and a lattice spacing a 0.079 fm [18,19].
The GF was implemented with the Runge-Kutta method; the results with time step dt = 0.01 and 0.001 agree. The GF flow time unit was fixed to t 0 /a 2 = 2.42, based on the criterion proposed in Refs. [22,23]: E t 0 = 0.3, where E is the mean energy density. The effect of the GF on the curves to be fitted in the slab method is shown in Figure 1 (left). As the GF proceeds, the fit has to be restricted to a narrower interval centered at x = 0.5. Moreover, the fitting function (4) has to be extended to where c is a constant (with respect to x), which increases roughly like 0.38 √ t [19]. The fits at different instances of the GF time, t = t 0 , 2t 0 . . . 8t 0 , yield very stable results for χ t , which are compatible both with a direct measurement (χ t a 4 = 7.8(2) · 10 −5 ), and with the method of Ref.
In the continuum, the Gradient Flow in O(N) models takes the form [24] ∂ t e(t, where t ≥ 0 is the GF time (of dimension [length] 2 ) and ∆ is the Laplace operator, which we handle by standard lattice discretisation. In order to proceed in discrete flow time steps, we apply the Runge-Kutta method. For a given configuration we first compute the gradients of all spin variables (at all In this manner, we considered L × L lattices, for instance with L = 120, β = 1.607, where the correlation length amounts to ξ = 19.1(2) [25]. Figure 2 (left) shows the (approximate) slab parabolae obtained for q 2 (x) in the sector |Q| = 1, in even multiples of the flow time unit t 0 = 0.0772 (which obeys E t 0 = 0.08, cf. Section 3). We see a qualitative difference from the QCD result in Figure 1: here the fits do not require any subtractive constant, i.e. the original formula (4) can be used, and χ t keeps on decreasing as the GF proceeds (the curvature of the parabola is reduced). This feature was also observed in all data sets to be reported in Section 3. The fitting results for χ t -obtained by the slab method, separately in the sectors |Q| = 0, 1, 2 -are close to the directly measured values, as we see in Figure 2 (right). In this case, the direct measurement is not problematic, since our simulations were carried out with the Wolff cluster algorithm, which proceeds in non-local update steps, thus suppressing the effect of topological freezing [26].
In order to investigate further these qualitatively different behaviours, we tested the 1d O(2) model, or quantum rotor, as a toy model. Here we refer to β = 2, as an example. In infinite volume we can compute analytically [27] (still for the standard lattice action, in lattice units) in agreement with our simulation results for size L = 100 (we used again the cluster algorithm).
Regarding the slab method under GF, Figure 3 shows the results for q 2 (x) in the sector |Q| = 2 at various flow times t, and χ t as a function of t. For the scaling quantity we obtain numerically at t = 0 : χ t ξ = 0.05388 (6), and at t = 10 : χ t ξ = 0.0496(1); thus we gradually approach the continuum value of χ t ξ = 1/(2π 2 Figure 3. The slab method in the 1d O(2) model: at short GF times, q 2 (x) (on the left) can be fitted to eq. (4), but long flow times require an extension to eq. (5), and the exclusion of x > ∼ 0 and x < ∼ 1 from the fit. After a long flow time, the sector Q = 0 provides the best results for χ t , as the plot on the right shows.

Topological scaling in the 2d O(3) model
In a d-dimensional quantum field theory with topological sectors, the dimensionless term χ t ξ d is supposed to be a scaling quantity, which converges to a finite value in the continuum limit. For the 1d O(2) model the continuum value χ t ξ = 1/(2π 2 ) is attained without problems [1,11,27]. In QCD, and in SU(N) Yang-Mills theories (N ≥ 2), a straight approach based on the expression χ t = x q 0 q x , where q x is the lattice topological charge density, faces problems. In conventional formulations one encounters a divergence due to the point x = 0, which prevents the continuum scaling. In these cases, there are known solutions to this problem, in particular the application of the GF [22,23].
In the 2d O(3) model, this question has been controversial in the 1980s and 1990s. The consensus is now that the quantity χ t ξ 2 seems to diverge in the continuum limit, as first observed in [5], and later underpinned by semi-classical studies, e.g. in Refs. [28,29]. Again the problem can be traced back to the topological density correlation at distance zero (see e.g. Ref. [1]), and the semi-classical picture suggests an abundance of very small topological windings, so-called dislocations, as the continuum limit is approached. Ref. [30] constructed and applied a sophisticated (truncated) classically perfect lattice action, with a host of couplings beyond nearest neighbour lattice sites, which suppress such dislocations. Nevertheless, the simulation results with this action still suggest a logarithmic divergence of the term χ t ξ d in the continuum limit (ξ → ∞ in lattice units). Hence the continuum limit of this popular model is generally assumed to be ill-defined, at least with respect to its topology (the reason is again the contribution q 0 q 0 , see e.g. Refs. [1,31,32]).
However, the question remains whether or not this divergence could be overcome by the GF; we gave preliminary results in Ref. [25]. In the following we summarise the status of this study. So far it involves nine L×L lattices, in the range of L = 24 . . . 404, where in each volume β has been tuned such that L/ξ 6. Therefore, increasing L corresponds to a controlled step towards the continuum limit, at a fixed and large physical box size. The statistics in each volume are 10 5 configurations (generated by the Wolff cluster algorithm, both in the single-cluster and the multi-cluster version). We repeat that we perform the GF with the Runge-Kutta 4-point method, with a time step of dt = 10 −4 , which is simultaneously applied to all spin variables.   value of the maximum of E t decreases, hence the reference value has to be sufficiently small, such that it is attained in all volumes under consideration. We are in the process of extending this study up to L = 494 and 606, where for instance 0.1 is not attained anymore. Hence we chose the reference value 0.08, which works up to L = 606 [33]. So we refer to the definition as we anticipated in Figure 2. Figure 4 (right) shows an example for the zero-momentum spin-spin correlation function C(r) at t = 0, 4t 0 , 8t 0 . At a fixed distance r between Euclidean time layers, C(r) increases as the GF proceeds, but the correlation length ξ remains unchanged within the errors. We measured ξ by a fit to a coshfunction in the range r ∈ [L/3, 2L/3]. We see that the GF -up to t = 8t 0 -hardly affects this long-range property. Figure 5 (left) illustrates the GF time evolution of the term χ t ξ 2 , which is supposed to be the scaling quantity. We see a rapid decrease at an early stage of the GF flow, in particular at large L and ξ. This observation is compatible with an increasing dominance of small dislocations on fine lattices: these topological windings are destroyed even by a short GF flow.
Finally, Figure 5 (right) shows χ t ξ 2 as a function of ξ, at various multiples of the flow time unit t 0 (the lines are drawn to guide the eye). At this stage, no trend towards a stabilisation is visible. The observed behaviour at t = 0 is well compatible with a logarithmically divergent function of the form χ t ξ 2 = c 1 ln(c 2 ξ + c 3 ); at t = t 0 . . . 6t 0 the quality of this fit is somewhat worse, but it still follows roughly this behaviour (although a power-law c 1 ξ c 2 + c 3 can be fitted with a similar quality) [25]. In any case, the present data do not provide a basis for revising the standard lore of a topologically ill-defined continuum limit in this model.

Summary and outlook
In Section 2 we have discussed the slab method, which enables a reliable measurement of the topological susceptibility within a fixed topological sector, i.e. from a Markov chain at fixed Q. This method does still work quite well when the GF is applied, but in some cases -in particular in QCD -we observed the necessity to subtract a constant from the expected fitting function.  Figure 5. Left: the GF time evolution of the term χ t ξ 2 in nine L × L lattices, in each case at L 6ξ. In large volumes we observe a rapid decrease at an early stage of the GF; this can be interpreted as the destruction of numerous small dislocations. Right: in the range t = 0 . . . 6t 0 , the quantity χ t ξ 2 does not seem to attain any finite continuum value as we approach the continuum (increasing ξ).
In the 2d O(3) model, the data presented in Section 3 do not suggest a continuum convergence of the term χ t ξ 2 after application of the GF. Further details, including a table with numerical results, are given in Ref. [25].
However, the impact rangex(t) of the GF is short in these examples. It can be estimated based on the heat kernel K(t, x); in d dimensions we obtain K(t, x) = 1 (2πt) d/2 exp − x 2 /(4t) ,x(t) = d d x x 2 K(t, x) In our 2d model it attains at mostx(6t 0 ) = √ 24t 0 1.62 at this stage of our study; this refers to our largest volume, L = 404, with ξ = 67.7 (3) x(6t 0 ). In order to arrive at conclusive results, we are now going to fix an extended GF time unit T 0 by a condition T 0 /ξ 2 = constant, and investigate flow times up to an impact range ofx(t) ≈ ξ/2. This study is in progress [33], and it should finally reveal whether or not the GF leads to a continuum scaling of the quantity χ t ξ 2 .