Properties and uses of the Wilson flow in lattice QCD

Theoretical and numerical studies of the Wilson flow in lattice QCD suggest that the gauge field obtained at flow time t>0 is a smooth renormalized field. The expectation values of local gauge-invariant expressions in this field are thus well-defined physical quantities that probe the theory at length scales on the order of sqrt(t). Moreover, by transforming the QCD functional integral to an integral over the gauge field at a specified flow time, the emergence of the topological (instanton) sectors in the continuum limit becomes transparent and is seen to be caused by a dynamical effect that rapidly separates the sectors when the lattice spacing is reduced from 0.1 fm to smaller values.


Introduction
Flows in field space are an interesting tool that may allow new insights to be gained into the physical mechanisms described by highly non-linear quantum field theories such as QCD. The flow B µ (t, x) of SU(3) gauge fields studied in this paper is defined by the equationṡ

JHEP08(2010)071
where A µ is the fundamental gauge field in QCD (see appendix A for unexplained notation; differentiation with respect to the "flow time" t is abbreviated by a dot). Evidently, for increasing t and as long as no singularities develop, the flow equation (1.1) drives the gauge field along the direction of steepest descent towards the stationary points of the Yang-Mills action.
In lattice QCD, the simplest choice of the action of the gauge field U (x, µ) is the Wilson action [1] S w (U ) = 1 g 2 0 p Re tr{1 − U (p)}, (1.3) where g 0 is the bare coupling, p runs over all oriented plaquettes on the lattice and U (p) denotes the product of the link variables around p. The associated flow V t (x, µ) of lattice gauge fields (the "Wilson flow") is defined by the equationṡ (1.4) in which ∂ x,µ stands for the natural su(3)-valued differential operator with respect to the link variable V t (x, µ) (see appendix A). The existence, uniqueness and smoothness of the Wilson flow at all positive and negative times t is rigorously guaranteed on a finite lattice [2]. Moreover, from eq. (1.4) one immediately concludes that the action S w (V t ) is a monotonically decreasing function of t. The flow therefore tends to have a smoothing effect on the field and it is, in fact, generated by infinitesimal "stout link smearing" steps [3]. The Wilson flow previously appeared in ref. [4] in the context of trivializing maps of field space. Familiarity with this paper is not assumed, but some mathematical results obtained there will be used here again. An important goal in the following is to find out whether the expectation values of local observables constructed from the gauge field at positive flow time can be expected to have a well-defined continuum limit. Evidence for the existence of the limit is provided by performing a sample calculation to one-loop order of perturbation theory directly in the continuum theory, using dimensional regularization, and through a numerical study of the SU(3) gauge theory at three values of the lattice spacing. Two applications of the Wilson flow are then discussed, one concerning the scalesetting in lattice QCD and the other the question of how exactly the topological (instanton) sectors emerge when the lattice spacing goes to zero.

Properties of the Wilson flow at small coupling
The aim in this section partly is to show that the Wilson flow can be studied straightforwardly in perturbation theory and partly to check that the expectation values of local gauge-invariant observables calculated at positive flow time are renormalized quantities.
For simplicity the perturbation expansion is discussed in the continuum theory using dimensional regularization. The gauge group is taken to be SU(N ) and it is assumed that there are N f flavours of massless quarks. As a representative case, the observable

JHEP08(2010)071
is considered and its expectation value is worked out to next-to-leading order in the gauge coupling.

Gauge fixing
The flow equation (1.1) is invariant under t-independent gauge transformations. E is therefore a gauge-invariant function of the fundamental field A µ and its expectation value can consequently be calculated in any gauge. In perturbation theory, the gauge invariance of the flow equation leads to some technical complications that are better avoided by considering the modified equatioṅ For any given value of the gauge parameter λ, the solution of eq. (2.2) is related to the one at λ = 0 through where the gauge transformation Λ(t, x) is determined bẏ The expectation value of E can thus be computed using the modified flow equation. Moreover, one is free to set λ = 1, which turns out to be a particularly convenient choice. Note that the use of the modified flow equation does not interfere with the fixing of the gauge of the fundamental field, since E is unchanged and therefore remains a gaugeinvariant function of the latter.

Solution of the modified flow equation
In perturbation theory the gauge potential is scaled by the bare coupling, 5) and the functional integral is then expanded in powers of g 0 . The flow B µ (t, x) thus becomes a function of the coupling with an asymptotic expansion of the form When this series is inserted in eq. (2.2), and if one sets λ = 1, a tower of equationṡ is obtained, where the expressions on the right are given by

JHEP08(2010)071
and so on. In particular, in D dimensions the leading-order equation implies which shows explicitly that the flow is a smoothing operation. More precisely, the gauge potential is averaged over a spherical range in space whose mean-square radius in four dimensions is equal to √ 8t. The higher-order equations (2.7) can be solved one after another by noting that (2.13) Recalling eqs. (2.9),(2.10), it is clear that this formula generates tree-like expressions, where the fundamental field A µ is attached to the endpoints of the trees.

Expansion of E
When the series (2.6) is inserted in a sequence of terms of increasing order in g 0 is obtained. The lowest-order term is 15) and the terms at the next order are Each of these terms is a power series in the gauge coupling, which may be worked out by expressing the coefficients B µ,k (t, x) through the fundamental field A µ (x) and by expanding the correlation functions of the latter using the standard Feynman rules. In practice it is advantageous to pass to momentum space by inserting the Fourier representations

JHEP08(2010)071
and the corresponding expressions for the higher-order fields B µ,k (t, x). The shorthand has been introduced here in order to simplify the notation. Note that the flow-time integral in eq. (2.19) is similar to a Feynman parameter integral. In particular, if Feynman parameters are used for the diagrams contributing to the gluon correlation functions, one ends up with gaussian momentum integrals that can be easily evaluated in any dimension D. The integrals over the parameters then look very much the same as the ones usually encountered, except for the fact that the flow-time parameters are integrated up to t rather than infinity.

Computation of E 0
In the case of the lowest-order term (2.15), the steps sketched in the previous subsection lead to the formula where D(p) µν denotes the unrenormalized full gluon propagator. Setting D = 4 − 2ǫ and choosing the Feynman gauge, the propagator assumes the form from which one infers that To this order, the computation is then easily completed by quoting the known result for the gluon self-energy, γ E = 0.577 . . . being Euler's constant.

Computation of E to order g 4 0
At the next-to-leading order, the expectation value of E receives contributions from E 0 and the terms E 1 , E 2 , . . . up to order g 4 0 generated by the expansion of the flow equation (cf. subsection 2.3). Some of the latter involve the 3-point gluon vertex, but in all cases only the leading-order gluon correlation functions are required.

JHEP08(2010)071
In the case of the term E 2 , for example, the purely algebraic part of the calculation leads to the integral where p = q + r. This integral appears to be quite complicated, but the polynomial in the numerator of the integrand allows the expression to be simplified and eventually to be evaluated analytically (see appendix B for further details). The result of the computation, shows that these contributions are not free of ultra-violet singularities. The computation of the other terms follows the same pattern and does not present any additional difficulties. Collecting all contributions, the result is then obtained.

Renormalization
The bare coupling g 0 is related to the renormalized coupling g in the MS scheme and the associated normalization mass µ according to [5] Now if eq. (2.28) is written as an expansion in the renormalized coupling, the terms proportional to 1/ǫ cancel and one obtains To this order in the gauge coupling, E thus turns out to be a quantity that does not need to be renormalized and which therefore encodes some physical property of the theory. In terms of the running coupling α(q) at scale q = (8t) −1/2 , the expansion (2.32) assumes the form

JHEP08(2010)071
In particular, for N = 3 one obtains The next-to-leading order correction is reasonably small in this case and corresponds to a change in the momentum q by a factor 2 or so if N f ≤ 3.

Lattice studies of the Wilson flow
In QCD the perturbation expansion (2.36) is expected to be applicable at small flow times only, where the smoothing range √ 8t is at most 0.3 fm or so. The properties of the Wilson flow at larger values of t can however be studied straightforwardly using the lattice formulation of the theory and numerical simulations.
The simulations of the SU(3) gauge theory reported in this section mainly serve to clarify whether E scales to the continuum limit as suggested by perturbation theory. Along the way, a new scale-setting method will be proposed based on the observed properties of E .

Simulation parameters
Ensembles of representative gauge fields were generated on three lattices, using the Wilson gauge action and a combination of the well-known link-update algorithms. The values of the lattice spacing quoted in table 1 derive from the results in lattice units for the Sommer reference scale r 0 = 0.5 fm [6] published by Guagnelli et al. [7]. At the chosen couplings β = 6/g 2 0 , the spacings of the three lattices thus decrease from roughly 0.1 to 0.05 fm by factors of 1/ √ 2. Moreover, the lattice sizes in physical units are approximately constant. In order to safely suppress any residual statistical correlations of the N cnfg generated fields, the separation in simulation time of the fields was taken to be at least 10 times the integrated autocorrelation time of the topological charge. The definition of the latter on the lattice is ambiguous to some extent, but its autocorrelation time is largely independent JHEP08(2010)071 of the choices one makes and is known to increase very rapidly when the lattice spacing is reduced [8,9]. In particular, already at a = 0.07 fm all other usual quantities of interest tend to be far less correlated in simulation time.

Observables
For any given gauge field configuration U (x, µ), the flow equation (1.4) can be integrated numerically up to the desired flow time and one may then construct gauge-invariant local observables from the gauge field at this time. In particular, is a possible definition of the density E on the lattice, P x being the set of unoriented plaquettes with lower-left corner x.
Another more symmetric definition of E is obtained by introducing a lattice version of the field tensor G µν (x) (see figure 1). E is then simply given by the continuum formula (2.1). Both definitions are equally acceptable at this point, since they respect all formal requirements (locality and gauge invariance in particular) and since they converge to the correct expression in the classical continuum limit.
For the numerical integration of the flow equation (1.4), any of the widely known integration schemes can in principle be used (see ref. [10], for example). The Euler integrators previously discussed in ref. [4] are particularly simple but also the least efficient ones. Evidently, the integration errors should be much smaller than the statistical errors of the calculated quantities. A fairly simple and numerically stable integrator that allows this condition to be easily met is described in appendix C.

Time dependence of E
To leading-order perturbation theory, the dimensionless combination t 2 E is a constant proportional to the gauge coupling. At the next order, the scale invariance of the theory is broken and t 2 E develops a non-trivial dependence on the flow time t. Asymptotic freedom actually implies that the combination slowly goes to zero in the limit t → 0. for the Λ-parameter in the MS scheme obtained by the ALPHA collaboration [11], and using the four-loop evolution equation for the running coupling α(q) [12], the perturbation series (2.36) can be evaluated at any value of the flow time given in units of r 0 . The curve obtained in this way is shown in figure 2 together with the error band that derives from the error of Λ quoted in eq. (3.2). Perhaps somewhat fortuitously, the simulation results obtained at a = 0.05 fm accurately match the perturbative curve over a significant range of t. The symmetric definition of E has here been used and for clarity the data from only one lattice are shown (as discussed below, the lattice-spacing effects are small and the data from the other lattices would therefore lie nearly on top of the line plotted in figure 2).
Beyond the perturbative regime, t 2 E grows roughly linearly with t, at least so within the range covered by the simulation data. The slowdown of the density E from the perturbative 1/t 2 to a smoother 1/t behaviour may perhaps be explained by noting that the Wilson flow tends to drive the gauge field towards the stationary points of the gauge action. In the vicinity of these points of field space, the right-hand side of the flow equation (1.4) is small and E consequently changes only little with time.

Lattice-spacing effects
Perturbation theory suggests that the density E scales to the continuum limit like a physical quantity of dimension 4. The scaling behaviour of E can be checked by introducing a reference scale t 0 through the implicit equation (see figure 2). If E is physical, the dimensionless ratio t 0 /r 2 0 must be independent of the lattice spacing, up to corrections vanishing proportionally to a power of a.
The data plotted in figure 3 clearly show that the ratio of reference scales smoothly extrapolates to the continuum limit. As may have been suspected, the lattice effects appear to be of order a 2 . They are at most a few percent on the lattices considered and particularly small if the symmetric definition of E is employed. At other points in time, the scaling violations behave similarly, but tend to increase towards small t, where the smoothing range √ 8t is only 2 or 3 times larger than the lattice spacing.

Synthesis
(a) Continuum limit. The numerical studies of the Wilson flow strongly support the conjecture that the gauge fields generated at positive flow time are smooth renormalized fields (except for their gauge degrees of freedom). In the case of QCD with a non-zero number of sea quarks, similar studies however still need to be performed. Additional confirmation in perturbation theory and through the consideration of further observables is evidently desirable.
(b) Scale setting. The time t 0 defined through eq. (3.3) may serve as a reference scale similar to the Sommer radius r 0 . With respect to the latter, t 0 has the advantage that its computation does not require any fits or extrapolations. Moreover, an ensemble of only 100 independent representative field configurations allows t 0 to be obtained with a statistical precision of a small fraction of a percent (for illustration, the values of t 0 /a 2 computed using the symmetric definition of E are listed in table 1).
(c) Universality. On the lattice the definition of local gauge-invariant quantities like E is not unique, but the results obtained in this section indicate that the differences become irrelevant in the continuum limit (see figure 3). This kind of universality is a consequence of JHEP08(2010)071 the fact, further elucidated in section 4, that the gauge fields generated by the Wilson flow are smooth on the scale of the lattice spacing. The universality classes of the local fields composed from the gauge field at fixed t/t 0 are therefore determined by the asymptotic behaviour of the fields in the classical continuum limit.

Functional integral and topological sectors
In lattice gauge theory, the space of gauge fields is connected and the concept of a topological sector has therefore no a priori well-defined meaning. However, one also knows since a long time that any classical continuous gauge field (including multi-instanton configurations) can be approximated arbitrarily well by lattice fields. The sectors are hence included in the field space but are not separated from one another. An understanding of how exactly the sectors get divided when the lattice spacing is taken to zero can now be achieved using the Wilson flow. The existence of the topological sectors thus turns out to be a dynamical property of the theory rather than being a consequence of an assumed or imposed continuity of the fields.

Field transformation
On a finite lattice, however large, the transformation is invertible and actually a diffeomorphism of the field space [4]. One can therefore perform a change of integration variables in the QCD functional integral from the fundamental field to the field at flow time t 0 . As already noted in ref. [4], the associated Jacobian can be worked out analytically and be expressed through the Wilson action along the flow. The expectation value of any observable O(V ) is then given by where S(U ) denotes the total action (including the quark determinants) of the theory before the transformation. Evidently, the functional integral may be rewritten in this way as an integral over the gauge field at any flow time t. Setting t to the reference time t 0 is however an interesting and natural choice. In particular, the fields that dominate the transformed integral then have a characteristic wavelength on the order of the fundamental low-energy scales of the theory. Since all terms in the action (4.3) have the same sign and favour smooth configurations, it is certainly plausible that large plaquette values s p are strongly suppressed. The probability for s p on a given plaquette p to be larger than some specified value s in fact decreases roughly like a 10 when the lattice spacing is reduced (see figure 4). In a finite volume of fixed physical size, the statistical weight of the field configurations with values of h beyond a given threshold is therefore rapidly going to zero in the continuum limit, while the fields that dominate the transformed functional integral become uniformly smooth in space.

Dynamical separation of the topological sectors
Many years ago, the space of lattice gauge fields was shown to divide into disconnected sectors if only fields satisfying a certain smoothness condition are admitted [13,14]. The proof is based on a geometrical construction of a local lattice expression for the topological charge which assumes integer values and which has the correct classical continuum limit.
In the case of QCD, the fields satisfying 1 h < 0.067 (4.5) are included in the subspace covered by the geometrical construction and thus fall into topological sectors very much like the continuous gauge fields in the continuum the-JHEP08(2010)071 ory. The space of fields "between the sectors" may accordingly be characterized by the inequality h ≥ 0.067. Recalling the discussion in subsection 4.2 of the smoothness properties of the fields that dominate the functional integral (4.2), the topological sectors are now seen to emerge as a consequence of the fact that the fields with large values of h have a rapidly decreasing weight in the integral when the lattice spacing is taken zero. In the case of the simulations reported in section 3, for example, the fraction of representative fields satisfying the bound (4.5) increases from 0% at a = 0.1 fm to 8% and 70% at a = 0.07 and 0.05 fm, respectively. The asymptotic behaviour of the weight of the fields between the sectors cannot be safely determined from these data, but the curves shown in figure 4 suggest that the probability of finding such a configuration on lattices of a given physical size goes to zero proportionally to a 6 .

Topological susceptibility
Once the functional integral (4.2) splits into a sum of integrals over effectively disconnected sectors, the assignment of the topological charge Q to the field configurations in a representative ensemble of fields becomes unambiguous. The geometrical construction mentioned before could be used to compute the charge, but in view of the universality property discussed at the end of section 3, a straightforward discretization of the topological density, using the symmetric lattice expression for the field tensor G µν , is expected to give the same results for the moments Q n in the continuum limit.
On the lattices listed in table 1, the second moment Q 2 computed along these lines is about 50 and the topological susceptibility χ t turns out to be independent of the lattice spacing within statistical errors. A fit by a constant then gives the result χ 1/4 t = 187.4(3.9) MeV, (4.6) which differs by less than two standard deviations from the value 194.5(2.4) MeV [16] obtained at flow time t = 0 using a chiral lattice Dirac operator and the index theorem [15]. Moreover, as shown in a forthcoming publication, the ratios of the fourth root of the second moments computed here and through a variant [18] of the universal formula proposed in [17] extrapolate to 1 in the continuum limit to within a statistical uncertainty of 3%. All these empirical results support the conjecture that the moments of the charge distribution in the functional integral (4.2) coincide with those obtained from the index theorem [15]- [18] and thus the ones appearing in the chiral Ward identities [19,20]. However, while highly plausible, the equality remains to be theoretically established.

Concluding remarks
The results reported in this paper shed some new light on the nature of non-abelian gauge theories and the continuum limit in lattice QCD. They are obviously incomplete in several respects and give rise to some interesting questions. In particular, an all-order analysis of the Wilson flow in perturbation theory will probably be required in order to achieve a structural understanding of its renormalization properties.

JHEP08(2010)071
The one-loop calculation in section 2 suggests that the fields obtained at positive flow time are renormalized fields for any number of quark flavours. So far the question was studied numerically only in the pure gauge theory, but it may be encouraging to note that the conjecture is true in QED, for any charged matter multiplet, as a consequence of the gauge Ward identity.
An intriguing aspect of the transformed functional integral (4.2) is the fact that it is dominated by smooth fields. Since the field transformation (4.1) is invertible, the integral nevertheless encodes all the physics described by the theory. Some properties (the division into topological sectors, for example) however become particularly transparent in this formulation of the theory, while its behaviour at high energies is better discussed in terms of the fundamental gauge field.

JHEP08(2010)071
acting on functions f (U ) of the gauge field are While these depend on the choice of the generators T a , the combination can be shown to be basis-independent.

B One-loop integrals
The Feynman integrals other than E 0 which contribute to E at order g 4 0 involve an integration over two momenta, q and r, and over at most two flow-time parameters. In all cases the integrand is of the form where P (q, r) is a Lorentz-invariant polynomial and s, u, v are linear combinations of the flow-time parameters.
As already mentioned, these integrals can in principle be evaluated by substituting the Feynman parameter representation for the propagators 1/q 2 , 1/r 2 and 1/(q + r) 2 . A more economic computation is, however, always possible using the basic integrals The second formula, for example, allows the integral q,r e −t(q 2 +r 2 +(q+r) 2 ) q 2 r 2 = 1 (4π) 4 (2t) 2 4) is an ordinary first-order differential equation of the forṁ where V t ∈ G and Z(V t ) ∈ g. The Runge-Kutta scheme described in this appendix obtains the solution of the flow equation at times t = nǫ, n = 1, 2, 3, . . ., recursively, starting from the initial configuration at t = 0. The rule for the integration from time t to t + ǫ is where Z i = ǫZ(W i ), i = 0, 1, 2. (C.3) Note that this rule is fully explicit. Moreover, since the gauge field can be overwritten from one equation to the next, and since Z 0 can be overwritten by 8 9 Z 1 − 17 36 Z 0 , intermediate storage space for only one of these latter fields is required.
A straightforward calculation shows that the integration scheme (C.2) is accurate up to errors of order ǫ 4 per step. The total error of the integration up to a specified flow time thus scales like ǫ 3 . Empirically one finds that the integration is numerically stable in the direction of positive flow time, the integration errors in the link variables being on the order of 10 −6 if ǫ = 0.01.
Open Access. This article is distributed under the terms of the Creative Commons Attribution Noncommercial License which permits any noncommercial use, distribution, and reproduction in any medium, provided the original author(s) and source are credited.