Chiral symmetry and the Yang-Mills gradient flow

In the last few years, the Yang-Mills gradient flow was shown to be an attractive tool for non-perturbative studies of non-Abelian gauge theories. Here a simple extension of the flow to the quark fields in QCD is considered. As in the case of the puregauge gradient flow, the renormalizability of correlation functions involving local fields at positive flow times can be established using a representation through a local field theory in 4+1 dimensions. Applications of the extended flow in lattice QCD include non-perturbative renormalization and O(a) improvement as well as accurate calculations of the chiral condensate and of the pseudo-scalar decay constant in the chiral limit.

1 Introduction In view of its renormalization properties [1,2], and since its application in lattice gauge theory is technically straightforward, the Yang-Mills gradient flow allows the dynamics of non-Abelian gauge theories to be probed in many interesting ways. The flow can be JHEP04(2013)123 used for accurate scale setting, for example, and it provides an understanding of how exactly the topological (instanton) sectors emerge in the continuum limit of lattice QCD [1]. 1 Moreover, observables at positive flow time are natural quantities to consider for nonperturbative renormalization and step scaling [3]- [6].
Matter fields may or may not be included in the flow equations. In this paper, a fairly trivial extension of the flow to the quark fields in QCD is considered, where the flow equation for the gauge field is unchanged, while the evolution of the quark fields as a function of the flow time is determined by a gauge-covariant heat equation (see section 2). The theoretical analysis and practical implementation of the flow is not significantly complicated by the inclusion of the quark fields. In particular, following ref. [2], the renormalization of correlation functions of gauge-invariant local fields at positive flow times can be shown to require no more than a multiplicative renormalization of the time-dependent quark fields once the parameters of QCD are renormalized as usual.
For illustration, two applications of the extended flow are worked out in this paper, one being a new strategy for the calculation of the axial-current renormalization constant in lattice QCD and the other a computation of the chiral condensate essentially through the evaluation of the expectation value of the scalar quark density at positive flow time. In both cases, the use of the flow proves to be technically attractive. The chiral condensate, for example, is easily obtained with high precision, because no additive renormalization is required.
In the next two sections, the extension of the gradient flow to the quark fields is discussed in the continuum theory. Since the flow equations respect chiral symmetry, the correlation functions of the time-dependent fields satisfy simple chiral Ward identities. As explained in section 4, these allow the correlation functions to be related to the physics of the light pseudo-scalar mesons. Then follows a more technical part of the paper (sections [5][6][7], where the flow is set up in the framework of lattice QCD. In particular, the associated field theory in 4+1 dimensions is shown to admit a well-defined local lattice regularization. The viability of the proposed applications of the flow in numerical lattice QCD is finally demonstrated through a sample calculation in 2+1 flavour QCD (sections 8,9).

The gradient flow in QCD
The theory considered in this paper is QCD with gauge group SU(N ) and a multiplet of two or more massive quarks in the fundamental representation of the gauge group. Many results are however of a fairly general nature and not limited to QCD. The theory is set up in Euclidean space and quantized through the functional integral as usual. In this and the following section, dimensional regularization is employed. The notational conventions are summarized in appendix A.

Flow equations
Let A µ (x) be the SU(N ) gauge potential, ψ(x) the quark field and ψ(x) the antiquark field integrated over in the functional integral. The latter carry Dirac, colour and flavour indices, which are usually suppressed for simplicity.
The Yang-Mills gradient flow evolves the gauge field as a function of a parameter t ≥ 0 that is referred to as the flow time. Starting from the fundamental gauge field, the time-dependent field B µ (t, x) is determined by the differential equation (see ref. [1] for an introduction to the subject). As already mentioned in section 1, the extension of the flow to the quark fields considered in this paper is a minimal one, where the evolution of the gauge field is left unchanged. The gauge field however appears in the flow equations 2 which, together with the initial conditions χ| t=0 = ψ, χ| t=0 = ψ, (2.6) define the time-dependent quark and antiquark fields χ(t, x) and χ(t, x). In the flow equations (2.4), the gauge-covariant Laplacian ∆ could be replaced by the square of the Dirac operator / D, for example, or be multiplied by a proportionality factor, but such modifications do not appear to offer any advantages.
The flow of quark fields introduced in this section is similar to the source smoothing operations proposed many years ago in refs. [7,8]. With respect to these popular methods, there are, however, a few important differences, one of them being the fact that the flow operates in four dimensions rather than on the fields on an equal-time hyper-plane. Moreover, the flow time varies continuously and the gauge field is not set to the fundamental gauge field, but is evolved together with the quark field.

Local composite fields
Since the flow equations are gauge-covariant, the time-dependent fields transform in the same way under gauge transformations as the fundamental fields. Examples of gaugeinvariant composite fields that are local both in space-time and in flow time are the densities E t (x) = − 1 2 tr{G µν (t, x)G µν (t, x)}, (2.7) 2 The left-action of any differential operator ∆ (or of a difference operator in lattice field theory) is defined through the requirement that the relation η † ∆ ← = (∆η) † holds for all complex-valued fields η.

JHEP04(2013)123
where r, s are flavour indices. Through the initial conditions (2.1), (2.6), the time-dependent fields depend on the fundamental fields. From the point of view of the QCD functional integral, composite fields like the densities (2.7), (2.8) are therefore observables, similar to Wilson loops or the ordinary pseudo-scalar densities, for example. In the following, the quantities of interest are the correlation functions of these fields, i.e. the expectation values of products of local fields composed from the basic fields at any flow time.

Perturbation theory
As explained in refs. [1,2], such correlation functions can be expanded in powers of the gauge coupling in a straightforward manner. First the flow equations are solved in powers of the fundamental fields. The substitution of these expansions in the local fields then leads to linear combinations of correlation functions of the fundamental fields that can be computed using the standard QCD Feynman rules.
In perturbation theory, the flow equations (2.2) and (2.4) are replaced by in order to avoid some technical subtleties. The terms proportional to the parameter α 0 > 0 serve to damp the gauge modes and could be removed by a time-dependent gauge transformation [1]. Correlation functions of gauge-invariant fields are therefore not affected by these extra terms.
The expansion of the time-dependent quark field χ(t, x) in powers of the fundamental fields is obtained by iterating the quark flow equation in its integral form, 12) and the corresponding integral equation for the time-dependent gauge field [1], where denotes the heat kernel of the Laplacian in D dimensions. At the lowest order in the gauge coupling, the interaction term in eq. (2.11) can be dropped and the expression is then obtained for the two-point function of the time-dependent quark field, where g 0 and M 0 are the bare coupling and (diagonal) quark mass matrix. The smoothing character of

JHEP04(2013)123
the quark flow is evident from these equations. Moreover, the smoothing range at leading order, √ 8t, is seen to be the same as in the case of the gauge field. There are 8 self-energy diagrams that contribute to the two-point function (2.14) at one-loop order of perturbation theory. The vertices in these diagrams derive from the QCD action and from the iteration of the integral equation (2.11). Following the lines of ref. [2], it is straightforward to compute the ultraviolet divergent parts of the diagrams and one then finds that the two-point function can be renormalized by renormalizing the quark masses as usual and the fields according to In the MS scheme in D = 4 − 2ǫ dimensions, the calculation yields for the field renormalization constant, g being the renormalized coupling.

Quark condensate at non-zero flow time
At non-zero flow times, the large momenta in the integral (2.14) are exponentially suppressed and the quark two-point function consequently has no singularities as (t, x) → (s, y). The absence of short-distance singularities is a general feature of the correlation functions at positive flow times, which derives from the smoothing property of the flow equations. In particular, the "time-dependent quark condensates" do not require additive renormalization. In perturbation theory, the renormalized condensates can be worked out in powers of the renormalized coupling g, with coefficients that depend on the renormalized quark masses m R,r , m R,s , . . . and the flow time t. The leading-order term in four dimensions is given by and thus vanishes in the chiral limit, consistently with the fact that chiral symmetry can only be spontaneously broken at the non-perturbative level.

Field theory in D+1 dimensions
The discussion in the following sections heavily builds on the fact that the correlation functions of the time-dependent fields coincide with the correlation functions in a local

JHEP04(2013)123
field theory in D+1 dimensions, the extra dimension being the flow time. The observation dates back to the seminal work of Zinn-Justin and Zwanziger [9,10] on the Langevin equation and recently allowed the renormalizability of the (pure-gauge) gradient flow to be established to all orders of perturbation theory [2]. As will become clear below, the inclusion of the quark fields in the theory in D+1 dimensions is straightforward. The general setup is exactly the same as in ref. [2], where an introduction to the subject and further details can be found.

Fields and action
In addition to the fundamental and the time-dependent fields already encountered, the theory in D+1 dimensions involves the Lagrange-multiplier fields L µ (t, x), λ(t, x) and λ(t, x). The latter are fermion fields with the same indices as the quark fields, while iL µ (t, x) is a vector field that takes values in the Lie algebra of SU(N ).
As explained in ref. [2], gauge fixing in perturbation theory requires the introduction of the Faddeev-Popov ghost fields c(x),c(x) and the associated time-dependent fields No boundary condition is imposed on the fieldd(t, x), which, in many ways, plays a rôle similar to the Lagrange-multiplier fields. Apart from the boundary conditions (2.1), (2.6) and (3.1), all fields are assumed to be unconstrained at this point. The total action of the theory, includes the gauge-fixed QCD action S of the fields at flow time zero and the bulk actions With respect to the pure-gauge flow studied in ref. [2], the bulk quark action S F,fl is the only new term in the total action.

Correlation functions
The n-point correlation functions of the gauge, quark, ghost and Lagrange-multiplier fields are formally defined through the functional integral over all fields with "Boltzmann factor" exp{−S tot }.
In the case of the correlation functions of the gauge and quark fields, the functional integral can be worked out, to some extent, by performing the integral over the Lagrangemultiplier fields. Since the action is linear in these fields, the integration yields a product JHEP04(2013)123 of functional Dirac δ-functions. The bulk actions were chosen such that the δ-functions amount to imposing the flow equations on the time-dependent gauge and quark fields. Up to mathematical subtleties and possibly non-trivial Jacobian factors, the correlation functions thus coincide with the ones calculated directly, as in section 2, by solving the flow equations with the initial conditions (2.1), (2.6) and treating the calculated time-dependent fields as QCD observables.
The equivalence of the theory in D+1 dimensions and the direct computation of the correlation functions can be rigorously shown in perturbation theory [2] and in lattice QCD (section 5). In particular, all correlation functions that do not involve the time-dependent quark fields or the associated Lagrange-multiplier fields coincide with the ones studied in ref. [2], where the quark fields were not evolved in flow time. This "reduction property" of the theory in D+1 dimensions strongly constrains the form of the possible counterterms required for renormalization (see subsection 3.4).

Fermion integral
Since the action in D+1 dimensions is quadratic in the fermion fields, the functional integral over these fields can be evaluated (at fixed bosonic fields) by applying Wick's theorem. The basic Wick contractions are worth being given explicitly, because the expressions provide some insight into the structure of the theory in D+1 dimensions and will, in any case, be frequently needed later.

Renormalization
In perturbation theory, the renormalization properties of the theory in D+1 dimensions can be determined using power counting and symmetry arguments. The task of figuring out all possible local counterterms to the action is simplified by the fact that the renormalization of the theory where only the gauge field is evolved in flow time is already known to require no more than the usual QCD parameter and field renormalizations [2]. In particular, at flow time t > 0, the fields B µ (t, x) and L µ (t, x) need not be renormalized. When the time-dependent quark fields are included in the theory, further counterterms may need to be added to the total action. The reduction property mentioned at the end of subsection 3.2 however implies that any such term must be proportional to the quark fields or the associated Lagrange-multiplier fields. Moreover, as explained in subsection 7.2 of ref. [2], bulk counterterms (i.e. local terms integrated over all D+1 coordinates) are excluded, because the correlation functions at large flow times are given by tree diagrams.
Taking the symmetries of the unrenormalized theory into account, and the fact that the quark and Lagrange-multiplier fields have power-counting dimension 3/2 and 5/2, respectively, a moment of thought reveals that the only possible additional counterterm is proportional to The insertion of this term in the Feynman diagrams is equivalent to a renormalization of the fermion and, similarly, the anti-fermion fields at positive flow time. Apart from the usual renormalization of the QCD parameters and fields, the renormalization of the theory in D+1 dimensions thus requires the time-dependent quark and associated Lagrange multiplier fields to be renormalized as specified above.

JHEP04(2013)123 4 Chiral symmetry relations
Since the quark flow equations (2.4) preserve chiral symmetry, the time-dependent composite fields like the densities (2.8) transform in a simple way under chiral rotations and thus allow the spontaneous breaking of the symmetry to be studied in useful new ways. In this section, the principal goal is to relate the time-dependent condensates (2.17) to the ordinary chiral condensate and the pseudo-scalar decay constant in the chiral limit.

Quark field equation and PCAC relation
The space-time dimension is now set to 4 and the need for regularization and renormalization is ignored for the moment, i.e. the statements made in this subsection are, strictly speaking, only valid at tree level of perturbation theory. Perhaps somewhat unexpectedly, the quark field equation in the theory in 4+1 dimensions, differs from the one in QCD by a term that couples to the quark fields at non-zero flow time. Equation (4.1) holds for any local fields φ 1 , . . . , φ n and the contact terms vanish if all points (t 1 , x 1 ), . . . , (t n , x n ) in 4+1 dimensions are different from (0, x). The validity of the field equation (4.1) follows from the identities which are a direct consequence of the formulae quoted in subsection 3.3 for the Wick contractions of the fermion fields. Equations (4.2)-(4.4) moreover show that contact terms only arise from contractions of the fundamental quark fields.

Now let
be the axial currents and densities in QCD. For r = s, the quark field equation then implies the PCAC relation where m 0,r , m 0,s are the bare masses of the quarks with flavour labels r, s and In particular, the identity

JHEP04(2013)123
holds at all flow times t > 0 and all points x including x = y. There are no contact terms in this case, because the fermion propagators that contribute to the correlation function on the left of eq. (4.8) all go from (0, x) to (t, y) in 4+1 dimensions or from (t, y) to (0, x). The contraction (4.2) thus does not occur. In view of the smoothing character of the flow equations, all parts of the correlation function are in fact regular functions of x.
Another purely algebraic consequence of the structure of the Wick contractions is the relation which, when combined with the PCAC relation (4.8), leads to an identity that allows the time-dependent quark condensates to be related to the physics of the pseudo-scalar mesons (see subsection 4.3). In the derivation of this equation, both the quark masses and the flow time were assumed to be positive, as otherwise the absence of boundary terms at large x and non-integrable singularities at x = y would not be guaranteed.

Renormalized PCAC relation
Beyond tree level of perturbation theory, the fields and quark masses in eqs. (4.8)-(4.10) require renormalization. The equations then become relations among renormalized correlation functions, whose exact form depends on the chosen normalization conventions. There are some natural choices one can make here and it is the goal in the following lines to specify these. The discussion implicitly assumes a sensible regularization of the theory, such as dimensional regularization or Wilson's lattice formulation [12], where most symmetries are preserved. As before, only the flavour non-diagonal channels are considered. The non-singlet axial currents and densities renormalize multiplicatively, 11) and the same is the case for the time-dependent densities As already mentioned, since there are no short-distance singularities at positive flow time, these fields (including the flavour diagonal ones) renormalize according to their quark content, i.e. like the product of a time-dependent quark and antiquark field at non-zero distance in 4+1 dimensions.
In the case of the fieldP rs , the multiplicative renormalizability, is less obvious, because the field has dimension 4 and could mix with other fields of dimension 3 and 4. The Wick contractions involving the Lagrange-multiplier fields λ and λ are however such that, up to contact terms, the two-point correlation functions ofP rs and any

JHEP04(2013)123
local field composed from the gauge and quark fields at flow time zero vanish. Divergent additive renormalizations by such fields are therefore excluded. One is then left with fields involving the Lagrange-multiplier fields, but sinceP rs is the only one among these with dimension 4 and the required symmetry properties, additive renormalizations are again not possible and the field must consequently renormalize multiplicatively. The normalizations of the renormalized fields and the renormalized quark masses m R,r may now be chosen such that the PCAC relation assumes the form Note that this convention only fixes the relative normalization of the axial current A rs R,µ , the product (m R,r + m R,s )P rs R and the densityP rs R . The normalization of the latter may however be fixed by requiring the equation and thus the identity Further normalization conditions then still need to be imposed on the renormalized quark masses and the time-dependent densities, but for the masses one can adopt one of the standard conventions, while the normalization of the time-dependent fields tends to cancel in the equations of interest and is left unspecified.
The normalization of the axial currents implicitly determined by eqs. (4.14), (4.15) coincides with the one usually chosen, where the axial charges assume integer values. There are different ways to show this, the perhaps most straightforward one starting from a formulation of lattice QCD that preserves chiral symmetry [13]- [23]. Using chiral Ward identities [24], the bare fields can be shown to satisfy eqs. (4.8), (4.9) in this case, up to terms vanishing proportionally to the lattice spacing. Equations (4.14), (4.15) then imply Z A =Z P = 1 and thus the canonical normalization of the renormalized axial currents. In view of the universality of the continuum limit, the canonical normalization is therefore guaranteed by these equations independently of the chosen regularization of the theory.

Pseudo-scalar matrix elements and the chiral condensate
The quark multiplet is now assumed to contain two light quarks, referred to as the u and d quarks, with the same renormalized mass m R . If there are not too many further quarks, chiral symmetry in the (u, d)-channel is expected to be spontaneously broken and there is a triplet of pseudo-scalar mesons, the "pions", whose mass vanishes as m R goes to zero.
At large distances, the light-quark pseudo-scalar correlation functions are dominated by the intermediate one-pion states. In particular, as x 0 → ∞ where ∆E denotes the energy gap to the higher states and M π , F π and G π are the pion mass, decay constant and vacuum-to-pion matrix element of the axial density. From the point of view of the theory in 4 dimensions, P du R,t (x) is not a local field, but it is still a well localized expression in the fundamental fields since the smoothing kernel K(t, x; 0, y) falls off approximately like a Gaussian at large distances |x − y|. At times x 0 much larger than the smoothing range √ 8t, the two-point function therefore decays in the same way as the ordinary pseudo-scalar correlation functions. The coefficient G π,t coincides with the vacuum-to-pion matrix element of a certain operator, but is here considered to be defined through eq. (4.19). The low-energy effective constants characterizing the (two-flavour) chiral limit, are the chiral condensate Σ and the value F of the pion decay constant at m R = 0. As already mentioned, these can be related to the time-dependent light-quark condensate Σ R,t = Σ uu R,t via the chiral Ward identity (4.16). Using the PCAC relation the latter may be written in the form The integral on the right of this equation diverges when m R → 0 and its asymptotic behaviour is determined by the pion pole of the correlation function in momentum space. Recalling eq. (4.19), one is then led to conclude that In the chiral limit, the time-dependent condensate Σ R,t thus provides an order parameter for the spontaneous breaking of chiral symmetry, which coincides with the standard condensate Σ up to a proportionality constant. A remarkable feature of eqs. (4.24), (4.25) is the fact that there is no reference to the axial currents, not even implicitly via the normalization conditions, and that the equations hold for any (positive) flow time t. Moreover, the renormalization constant Z χ drops out so that the decay constant F , for example, is directly given through the bare time-dependent condensate and pseudo-scalar vacuum-to-pion matrix element.

Chiral perturbation theory
Close to the chiral limit, the asymptotic behaviour of many quantities in QCD can be described by chiral perturbation theory. If the Compton wavelength of the pion is much larger than the smoothing range of the gradient flow, i.e. if the time-dependent densities S rs R,t and P rs R,t are local fields from the point of view of pion physics, with the same symmetry properties as the densities at vanishing flow time. At low energies, the correlation functions of the time-dependent fields and the ordinary axial currents and densities are therefore expected to be dominated by the pion intermediate states and should hence be accessible to a quantitative description in the framework of chiral perturbation theory.
The time-dependent densities can in fact easily be included in chiral perturbation theory by adding the appropriate source terms to the chiral Lagrangian [25]. At each order of the chiral expansion, these terms require a number of effective couplings to be introduced, which will, in general, depend on the flow time. The expansion of the chiral condensate Σ R,t and the matrix element G π,t , for example, can then be worked out and provides insight into how exactly the chiral limits (4.24), (4.25) are reached.

Lattice regularization
The starting point in this section is a formulation of lattice QCD, where the quark fields are put on the lattice as proposed by Wilson [12]. No particular assumptions are made on the lattice Dirac operator or the gauge action, but both should be local and respect all symmetries preserved by the standard Wilson theory. For notational convenience, the lattice spacing a is set to unity. The lattice is taken to be infinitely extended in all directions unless specified otherwise.

Flow equations
On the lattice, the time-dependent gauge field V (t, x, µ) is defined by where U (x, µ) is the fundamental lattice gauge field and (∂S w )(V, x, µ) the gradient of the Wilson plaquette action (see appendix A and ref. [1] for further explanations).
In the case of the quark fields, the boundary conditions and the form of the flow equations are the same as in the continuum theory, but the operator ∆ in eq. (2.4) gets replaced by the lattice Laplacian ∇ µ and ∇ * µ being the gauge-covariant forward and backward difference operators in presence of the time-dependent gauge field.

JHEP04(2013)123
As already noted in ref. [1], the global existence, uniqueness and smoothness of the solution of the flow equation (5.2) is rigorously guaranteed. Since ∆ is a bounded operator, and since it depends smoothly on the flow time, the same can be shown to be true in the case of quark flow equations [26]. Local fields at non-zero flow time, such as those considered in the previous sections, are thus completely well defined on the lattice.

Lattice theory in 4+1 dimensions
Similarly to the continuum theory, the lattice QCD correlation functions of the timedependent local fields can be obtained from a lattice field theory in 4+1 dimensions. To be able to fully control this construction, the latter is first set up for a discretized version of the flow equations, where the flow time is restricted to a finite set of values separated by a time step ǫ > 0. The integration of the discretized evolution equation for the gauge field amounts to an approximate integration of the continuous equation (5.2) using the Euler method (or, equivalently, through the application of a sequence of "stout smearing" steps [27]). In the case of the quark fields, the discretized flow equation has the same form as the continuous one, but the time derivative of the quark field is replaced by the forward difference Note that the discrete equations are such that they can be solved in steps, starting from the fundamental fields at time t = 0 and recursively proceeding from time t to t + ǫ. In particular, the solutions are uniquely determined by the initial values of the fields.
In the lattice theory in 4+1 dimensions, all fields integrated over in the functional integral are defined at the times (5.4) only. Since gauge fixing is not needed in the present context, there are no ghost fields and no gauge terms in the action, but the boundary conditions (2.6) and (5.1) are imposed as before. The components of the Lagrange-multiplier fields at time T decouple from the other fields and can therefore be set to zero (or be omitted from the beginning).
A possible choice of the bulk actions is then maps any complex N × N matrix to an element of the Lie algebra of SU(N ). Apart from these actions and the lattice QCD action, the measure term must be included in the total action of the lattice theory (see appendix B for the definition of the Jacobian matrix J(W ) ab of the mapping (5.9) and the neighbourhood N ⊂ SU(N ) of unity).
Having specified the fields and their action, the functional integral of the lattice theory in 4+1 dimensions is defined in the standard manner. Note that the measure term (5.10), (5.11) effectively restricts the link variables in the functional integral to the domain where Since the action S G,fl is purely imaginary, the integral over the Lagrange-multiplier field L(t, x, µ) is not absolutely convergent, but it is understood that the integration ambiguity this may entail is removed through the addition of an infinitesimal quadratic term in the Lagrange-multiplier field to the action.

Recovering the theory in four dimensions
The correlation functions defined in this way of local fields composed from the fields V , χ and χ can be evaluated, to some extent, by first integrating over the Lagrange-multiplier fields. Since the total action is linear in the latter, the integration yields a product of Dirac δ-functions. The arguments of these δ-functions vanish if V , χ and χ satisfy the discrete flow equations, and for sufficiently small step sizes ǫ there is in fact no other way the arguments can be made to vanish. 3 In the case of the gauge field, the important point to note here is that the mapping P in the δ-functions operates on SU(N ) matrices in a small neighbourhood of unity. The properties of the mapping quoted in appendix B then

JHEP04(2013)123
imply that the constraints imposed by the δ-functions are equivalent to imposing the flow equation (5.5) on the gauge field. The integral over the field variables at all flow times t > 0 can now be performed using the δ-functions. Starting with the fields at time t = T and proceeding from there to the smaller times, the δ-functions allow all fields to be eliminated recursively until the only integration variables left are the fundamental fields U , ψ and ψ. In this process, the integrals over the δ-functions for the gauge field variables give rise to a product of Jacobians of the mapping P, but these factors are canceled by the measure term in the action (see appendix B).
The calculation thus shows that the correlation functions in 4+1 dimensions of local fields that do not involve the Lagrange-multiplier fields are exactly equal to the correlation functions of the same fields defined in 4 dimensions through the discrete flow equations and the QCD functional integral.

Continuous time limit
At this point, the time-dependent fields in the correlation functions are still the ones obtained by the Euler integration of the flow equations. Since the associated integration errors are of order ǫ, the continuous-time correlation functions are recovered when ǫ → 0 (and, trivially, T → ∞). Moreover, the measure term vanishes in this limit. In the following, the continuous-time formulation of the lattice theory will normally be considered, but the functional integral in 4+1 dimensions is always assumed to be defined through the limit ǫ → 0 of the discrete-time functional integral.

On-shell O(a) improvement
The lattice-spacing dependence of correlation functions involving the fields at non-zero flow time is best discussed in the theory in 4+1 dimensions. In this framework, the standard argumentation based on locality, power counting and symmetry can be applied. The relevant Symanzik local effective theory [28,29] then coincides with the continuum theory in 4+1 dimensions, with O(a) and higher-order terms in the lattice spacing a added to the action and the local fields. Moreover, on-shell improvement [30] means that the cancellation of latticespacing effects is limited to correlation functions at non-zero distances in 4+1 dimensions.
Although the theoretical discussion is more widely applicable, the O(a)-improved Wilson formulation of lattice QCD [31,32] is assumed from now on. Where possible, the same conventions and notation as in ref. [32] are used.

Effective action
Since the flow equations in the lattice theory are classically O(a) improved, and since loop diagrams do not contribute to the correlation functions at asymptotically large flow times,

JHEP04(2013)123
all O(a) terms in the Symanzik effective action must be boundary terms, i.e. local terms at flow time zero. Furthermore, if the effective theory is to describe on-shell quantities only, as is the case here, many terms can be eliminated using the quark field equations, the flow equations and the boundary conditions (2.1), (2.6). In the chiral limit, a possible choice of the O(a) term in the effective action is then where F µν (x) is the field tensor of the fundamental gauge potential (here and below, c 1 , c 2 , . . . denote coefficients of the effective theory). Quite many more terms are required to describe the O(a) lattice-spacing effects at non-zero quark masses. In the theory in 4 dimensions, these were completely classified by Bhattacharya et al. [33]. The Symanzik effective action in 4+1 dimensions in addition includes a term that contributes to correlation functions involving the bulk fields, M being the quark mass matrix in the effective theory. Its effect on the correlation functions is equivalent to a renormalization 3) of the bulk fermion (and anti-fermion) fields.

Effective local fields
The form of the effective fields representing the local lattice fields in the Symanzik theory is constrained by the same general principles that determine the form of the effective action. In particular, the effective fields representing fields at positive flow time have no O(a) terms. Such terms however occur in the case of the local fields at flow time zero. Considering again the theory with only massless quarks, the effective fields representing the axial currents with flavour indices r = s, for example, are given by where the ellipsis stands for terms of order a 2 . The associated effective pseudo-scalar densities, have a similar expansion.

JHEP04(2013)123
The mass-dependent O(a) corrections to these fields coincide with the ones in the theory in 4 dimensions [33]. In general, the determination of the O(a) terms however requires a careful consideration of a possible mixing of fields and may lead to fairly complicated expressions.

O(a) improved lattice action
Following common practice, the mass-dependent O(a) counterterms are omitted in the improved action and are instead included in the parameter and field renormalization factors [32]. The O(a) counterterms in 4+1 dimensions are then given by where F µν (x) denotes the standard clover expression for the gauge field tensor F µν (x) on the lattice. While the first term in eq. (6.7) (the Sheikholeslami-Wohlert term [31,32]) is already required for the improvement of the theory in 4 dimensions, the one proportional to the improvement coefficient c fl (g 0 ) is needed to cancel the O(a) lattice effects in correlation functions involving the bulk fields.
Having specified the action, the Wick contractions of the basic fermion fields in the O(a) improved theory can be worked out straightforwardly (see appendix C). In the continuous time limit, the contractions are practically the same as in the continuum theory with the obvious modifications (integrals over the space-time coordinates are replaced by sums over lattice points, the quark propagator S(x, y) by the inverse of the massive O(a) improved lattice Dirac operator and the kernel K(t, x; s, y) by the fundamental solution of the quark flow equation on the lattice). A special case is the contraction of the time-dependent quark fields, which is the only one depending on the improvement coefficient c fl . In particular, the contraction and thus the correlation functions of the fields ψ and ψ are independent of c fl , as has to be the case, since the theory in 4 dimensions is already improved through the inclusion of the Sheikholeslami-Wohlert term. A perhaps puzzling aspect of eq. (6.8) is the fact that the limit lim t→0 χ(t, x)χ(s, y) = ψ(x)χ(s, y) − c fl K(s, y; 0, x) † (6.9) differs from ψ(x)χ(s, y) by a term of order a, which is not just a contact term. The subtractions required for O(a) improvement thus depend on whether the quark fields in the correlation functions are at zero and non-zero flow time. This is actually not too surprising since these fields also renormalize differently. Both differences merely reflect the fact that the continuum limit of correlation functions involving the time-dependent fields must be taken at fixed flow times given in physical units and that the Symanzik effective theory in 4+1 dimensions correctly describes the deviation of the lattice theory from its continuum limit only when the latter is approached in this way.

Renormalized improved fields
As usual, the additively renormalized bare quark masses m q,r are defined by m q,r = m 0,r − m c , (6.10) where m c denotes the critical bare mass in the theory with mass-degenerate quarks. It is also helpful to introduce the associated subtracted bare mass matrix M q and the combinations m q,rs = 1 2 (m q,r + m q,s ) (6.11) of quark masses. The lattice fermion fields at positive flow time require multiplicative renormalization and O(a) improvement by a mass-dependent factor (cf. subsection 6.1). Quark and antiquark fields are scaled with the same factor, while the Lagrange-multiplier fields λ, λ are renormalized with the inverse of the same factor. 4 Time-dependent composite fields like the pseudo-scalar density P rs R,t = Z χ (1 + b χ m q,rs +b χ trM q )P rs t (6.13) then renormalize according to their quark content. Starting from the bare fields and these are renormalized according to [32,33] A rs R,µ = Z A (1 + b A m q,rs +b A trM q )A rs I,µ , (6.20) P rs R = Z P (1 + b P m q,rs +b P trM q )P rs I . (6.21) The new improvement coefficients,c A andc P , are required for the O(a) improvement of correlation functions involving the fermion fields at non-zero flow time, but all other coefficients already occur in the theory in four dimensions [31]- [33]. At tree-level of perturbation theory, All coefficientsb X are, incidentally, of order g 4 0 [33]. JHEP04(2013)123

Improvement and renormalization of the fieldP rs
In the chiral symmetry relations derived in section 4, the fieldP rs plays a prominent rôle. The improvement and renormalization is a bit complicated in this case and is therefore discussed separately from the other fields. As already noted in subsection 4.2, the field can only mix with local composite fields that include at least one Lagrange-multiplier field. Charge conjugation, the lattice symmetries and the flavour symmetry then imply that the field is multiplicatively renormalizable. There are, however, several fields that can mix with the density at order a, among them two fieldsP that have not appeared before. Inspection then shows that a possible choice of the improved and renormalized densities is P rs I =P rs +ĉ A∂µÃ rs µ +ĉ PP rs , (6.25) P rs R =Z P (1 + bP m q,rs +bP trM q )P rs I +bP (m q,r − m q,s )Q rs , (6.26) where, as usual, the field equations were used to reduce the number of terms. The fields on the right of these equations are distinguished by their symmetry properties or content in Lagrange-multiplier fields. They are all multiplicatively renormalizable and the contributions of the counterterms proportional toÃ rs µ ,P rs and Q rs in on-shell correlation functions are therefore of order a.
The expression (6.25), (6.26) for the renormalized improved density can be simplified by noting that the relations x P rs (x)P sr t (y) = S rr t (y) + S ss t (y) + 2c fl x P rs (x)P sr t (y) , (6.27) x Q rs (x)P sr t (y) = S ss t (y) − S rr t (y) , (6.28) hold exactly, for any t > 0, as a consequence of the form of the Wick contractions of the fermion fields. On the other hand, as discussed in subsection 4.2, the renormalization constants can be (and are to be) chosen so that the normalization condition (4.15) is satisfied in the continuum limit. Since the correlation functions in eq. (4.15) converge to the continuum limit with a rate proportional to a 2 , the comparison with the unrenormalized identities, eqs. (6.27), (6.28), then shows that The renormalized densitỹ thus assumes a fairly simple form, in whichĉ A is the only new improvement coefficient.

JHEP04(2013)123
6.6 PCAC relation on the lattice In the continuum limit, the renormalized PCAC relation (4.14) holds provided the renormalization constants Z A , Z P and the renormalized quark masses m R,r are chosen appropriately. If also all improvement coefficients are properly tuned, this implies (6.32) where, following the tradition, the divergence of the improved axial current is taken to be A technical detail worth emphasizing here again is the fact that the PCAC relation (6.32) holds at all x, including x = y, as long as the flow time t is set to a positive value given in physical units.

Calculation of correlation functions
In numerical lattice QCD, correlation functions involving fields at non-zero flow time can be computed following the steps usually taken in the case of ordinary hadronic correlation functions. The fact that the flow equations must be solved at some point nevertheless represents a complication that needs to be carefully considered. For illustration, the details are worked out in this section for two representative cases.

Pseudo-scalar two-point function
In practice one is interested in the correlation function at vanishing spatial momentum and computes its average over a lattice Γ of source points, using random source fields [34], in order to reduce the statistical errors. If x is taken to be the source point, the averaging amounts to the substitution (7.2) where N Γ is the number of points in Γ and are randomly chosen complex spinor fields on Γ with mean zero and variance

JHEP04(2013)123
In a QCD simulation, the computation of the averaged correlation function at zero spatial momentum, 1 N Γ x∈Γ y P rs (x)P sr t (y) = − 1 N Γ N s Ns k=1 y φ k,s (t, y) † φ k,r (t, y) , (7.5) then amounts to calculating the functions for a representative ensemble of gauge fields, all k = 1, . . . , N s and h = r, s. For a given gauge field, the computation of φ k,h (t, y) proceeds in two steps, where one first calculates φ k,h (0, w) by solving the lattice Dirac equation with mass m 0,h and source field η k (x). The calculated field must then be evolved in flow time from time 0 to time t. In the direction of increasing flow time, the flow equations are numerically stable and the solution can easily be obtained, with negligible integration errors, using a Runge-Kutta integrator (appendix D).

Chiral condensate
In the case of the time-dependent quark condensate, an averaging over the position x is again desirable. Introducing random source fields as above, one is then led to the representation in terms of the fields ξ k (t; s, w) = x K(t; x; s, w) † η k (x) (7.9) at flow time s = 0. The computation of these fields requires the adjoint flow equation to be solved from time s = t (where ξ k coincides with η k ) to s = 0. A straightforward Runge-Kutta integration should not be used here, because the Laplacian ∆ is the one in presence of the gauge field determined by the gradient flow at time s, which would therefore have to be evolved backward in time, i.e. in the unstable direction. As explained in appendix E, this difficulty can be overcome through a hierarchical scheme that avoids the backward integration of the gauge field.

JHEP04(2013)123
8 Strategy for the computation of Z A and c fl The PCAC relation (6.32) holds provided the relevant renormalization constants and improvement coefficients are chosen appropriately. As explained below, the constant Z A and the coefficient c fl can conversely be determined, up to terms of order a 2 and a, respectively, by requiring the relation to be satisfied exactly at two points in the space of kinematical parameters. The computation through numerical simulation of Z A and c fl along these lines assumes that the improvement coefficients c sw and c A have already been calculated. Periodic boundary conditions should be imposed in the space directions, but apart from this no further assumptions are made. In particular, the boundary conditions in time do not need to be specified, since the computed values of Z A and c fl depend on them only at the level of the lattice-spacing effects.

Integrated form of the PCAC relation
The O(a) improvement unfortunately leads to an inflation of terms when the renormalized PCAC relation is expanded in the correlation functions of the unimproved bare fields. A slight simplification can however be achieved by summing the equation over an interval of time x 0 around y 0 and over all points x on the spatial lattice. One of the bare correlation functions entering the relation is then and there are 6 further such correlation functions in which P rs is replaced bỹ P rs ,P rs , Q rs ,∂ µ A rs µ ,∂ µÃ rs µ or ∂ * µ ∂ µ P rs .
In the case of the last three fields in this list, the sum over x 0 reduces to a sum of terms at the boundary of the summation range as in for example. For notational convenience, it is now helpful to introduce the abbreviationŝ

4)
Z rs P = Z P 1 + b P m q,rs +b P trM q . the integrated PCAC relation then becomeŝ Z rs A C rs ∂A +c A C rs ∂Ã + c A C rs ∂∂P − m rs C rs P +c P C rs This equation still looks fairly complicated, but one should keep in mind that it is a linear relation among computable correlation functions, which must hold for all flow times t > 0, time ranges d ≥ 0 and quark flavours r = s. The unknown coefficients are thus strongly constrained by the identity. At values of d larger than the smoothing range √ 8t of the gradient flow, C rs ∂Ã falls off like a Gaussian and rapidly becomes negligible with respect to the other terms in eq. (8.7). Moreover, the sum of the terms in the curly bracket as well as all other terms in the equation are then practically independent of d. Without further notice, only this range of d is considered in the following and the two terms proportional to C rs ∂Ã are dropped.

Up-down flavour channel
In QCD with mass-degenerate up and down quarks, the PCAC relation in the up-down channel assumes the form Clearly, as long as the chosen flow times are in the scaling range and are held fixed in physical units, the calculated ratios do not depend on them, up to terms of order a 2 and a, respectively. At quark masses close to their physical values, and lattice spacings a ≤ 0.1 fm, the factor 1 −Ẑ ud Ac P m ud differs from unity by a fraction of a percent. An estimate of the size of this correction term can be obtained by inserting the known value of m ud and the tree-level value (6.22) ofc P .Ẑ ud A and c fl can thus be computed up to this small correction factor. Note that the bare pion decay constant determined through the vacuum-to-pion matrix element of the improved axial current is to be renormalized with the mass-dependent fac-torẐ ud A rather than Z A . The difference is however small and can be estimated using the tree-level value b A = 1.

Up-strange flavour channel
The additional term that appears in the PCAC relation Ac P m us C us P (8.10)

JHEP04(2013)123
in this channel is proportional to the square of the mass difference m q,u − m q,s and is therefore likely to be negligible in practice. To be able to determine the ratios (8.9) (with ud replaced by us), eq. (8.10) must otherwise be evaluated at three values of the flow time t. For a ≤ 0.1 fm and physical quark masses, the correction factor 1 −Ẑ us Ac P m us may deviate from unity by up to one percent or so. The comparison with the results obtained in the ud channel, unfortunately only allows a certain linear combination of the coefficients b A andc P to be estimated.

Sample calculation in 2+1 flavour QCD
The aim in the following paragraphs is to demonstrate the feasibility of the proposed strategies for the computation of the renormalization constant Z A , the improvement coefficient c fl and the low-energy constants Σ and F . A single simulation run suffices for this test and no attempt is made to actually compute the low-energy constants, as this would require the data to be extrapolated to the infinite-volume, chiral and continuum limit.

Simulation details
The calculations reported below are based on a recent simulation of QCD with 2+1 flavours of quarks [35]. In this run (labeled I 1 in ref. [35]), a 64 × 32 3 lattice with open boundary conditions in time [36] was simulated at a point in parameter space, where the lattice spacing is estimated to be 0.090 fm [37,38] and where the pion and kaon masses are about 203 and 520 MeV, respectively [35]. The physical size of the lattice is thus such that the finite-volume effects on the calculated meson masses and matrix elements are expected to be small. For the Iwasaki gauge action [39] and the standard O(a)-improved Wilson-Dirac operator [31,32] used in the simulation, the coefficients c sw and c A are known nonperturbatively [40,41]. A representative ensemble of 150 gauge-field configurations, separated by 9.6 units of molecular-dynamics time, was produced in the run. In order to suppress any residual autocorrelation effects, the statistical errors quoted in this section were estimated by dividing the measurement data series for the primary observables into bins of 5 consecutive measurements and by applying the jackknife method to the binned data. All correlation functions in the PCAC relation can easily be obtained with small statistical errors over the whole range of the summation window d. As expected, the sum of correlation functions multiplying the axial-current renormalization constant and the other correlation functions reach a plateau at values of d larger than 2 or 3 times the smoothing range of the gradient flow (see figure 1 for a representative case; the slight bending down of the function multiplyingẐ ud A close to the boundary of the lattice is a lattice-spacing effect). Following the strategy outlined in section 8, the ratios (8.9) are now determined by requiring the PCAC relation (8.8) to hold exactly at two values t 1 , t 2 of t. The choice of the summation window d has very little influence on the calculated ratios as long as one stays within the plateau range of the correlation functions. Setting d = 18 then leads to the results quoted in the second and third column of table 1, where the factor 1 −Ẑ ud Ac P m ud was estimated using the tree-level value (6.22) of the coefficientc P . As anticipated, this correction is very small (about 0.2%).

Computation ofẐ rs
In the up-strange flavour channel, there are two sources of O(a) mass corrections. Proceeding as above, the correction coming from the factor 1 −Ẑ us Ac P m us is estimated to be approximately 1.3%, while the one deriving from the term in eq. (8.10) proportional to b χ changes the results by about 0.2%. As can be seen from table 1, there is a remarkable consistency among the values of c fl obtained in the two flavour channels. There is also little dependence on the flow times and the difference of the renormalization constantŝ  The results for the renormalization constantẐ ud A quoted in table 1 are in agreement with the value Z A = 0.781 (20) previously obtained in ref. [42] using a method based on the Schrödinger functional [43,44]. Note that the improvement coefficient c fl turns out to be close to its tree-level value (6.22), although there is no very good reason for this to be so at the gauge coupling considered.

Computation of the chiral condensate
Having determined the improvement coefficient c fl , the unrenormalized time-dependent condensate Σ uu t can be computed straightforwardly as described in subsection 7.2. The source time x 0 is taken to be in the middle of the lattice in order to minimize the excited-states contributions to the expectation value of the scalar density. Using 12 random source fields, and evaluating c fl at t 1 , t 2 = 3.86, 5.56, the calculation yields the results quoted in the second column of table 2. In the range of flow time considered, the time-dependent condensate is thus seen to decrease, roughly like 1/t, when t increases. The statistical error follows this evolution and is, in any case, encouragingly small.
To be able to relate the time-dependent condensate to the quark condensate Σ in the chiral limit, the vacuum-to-pion matrix elements G π = Z P (1 + b P m q,ud +b P trM q )G ud , (9.1) need to be computed as well. Actually, since only the ratio appears in eqs. (4.24), (4.25), it suffices to calculate G π and the unrenormalized matrix element G ud t . The latter can be determined from the pseudo-scalar correlation function (7.1) at large time separations |x 0 − y 0 | in the same way as the bare matrix element G ud at vanishing flow time is usually extracted from the ordinary pseudo-scalar two-point function. Setting x 0 to a value next to the boundaries of the lattice, there is in all cases a wide range in time y 0 , where the one-pion intermediate states dominate and the desired matrix elements can be easily determined. Similarly to the time-dependent condensate, the matrix element G ud t is a monotonically decreasing function of t. The ratios listed in the last two columns of table 2 however vary much less with t, consistently with the fact that the flow-time dependence of the ratios must disappear in the chiral limit.

Remarks
Conversion to physical units. The spacing of the simulated lattice, a = 0.08995(40) fm, was determined by PACS-CS through a computation of the mass of the Ω baryon [38]. Using step scaling and the Schrödinger functional, PACS-CS also calculated the renormalization constant Z P = 0.580 (21) [42] that relates the bare matrix element G ud to G π in the MS scheme at 2 GeV. 5 At t = 3.86, for example, the results Σ R,t G π,t = 91(2) MeV, (9.5) are then obtained, where the O(a) mass correction in eq. (9.1) was neglected (the error in eq. (9.4) is anyway dominated by the error of Z P ). Clearly, while these ratios may be close to the condensate Σ and decay constant F in the chiral limit, simulations of a range of lattices will be required to be able to determine Σ and F with controlled systematic errors.
Pseudo-scalar decay constants. Since the calculation of the axial-current renormalization constant does not require separate simulations, the renormalized pseudo-scalar meson decay constants become directly accessible on each simulated lattice. Depending on the flavour channel and the precision that is to be attained, the mass-dependent O(a) corrections in the PCAC relation used to determine the renormalization constantẐ rs A may however need to be estimated non-perturbatively.
Scaling to the continuum limit. The results of the computations reported in this section depend on a choice of flow times. When lattices with different spacings are simulated, these parameters should be held fixed in units of a suitable reference scale such as the reference flow time t 0 [1]. Only then can the calculated renormalized quantities be expected to converge to the continuum limit with a rate proportional to a 2 .
Through its extension to the quark fields, the gradient flow becomes a powerful tool for studies of the spontaneous breaking of chiral symmetry in QCD and the physics of the light pseudo-scalar mesons. The focus in this paper has been on the theoretical framework and its implementation in one of the widely used formulations of lattice QCD. In view of the renormalizability properties of the flow, the choice of the lattice regularization is however expected to be irrelevant in the continuum limit as long as the locality and gauge invariance of the lattice theory are guaranteed. The theoretical analysis obviously applies to field theories with other gauge groups and fermion multiplets as well.
In the cases worked out so far, the application of the gradient flow in numerical lattice QCD offers important technical advantages with respect to the established computational strategies. Apart from the ones already mentioned in section 1, there are potentially many further uses of the flow. The flavour singlet channel, for example, has not been considered yet and the time-dependent condensate is likely to be an easily accessible chiral order parameter in QCD at non-zero temperatures. Moreover, the application of the flow in QCD-like theories near the conformal window may conceivably lead to interesting qualitative insights.

JHEP04(2013)123
Scalar products γ µ p µ are abbreviated by / p, as usual, and σ µν = i 2 [γ µ , γ ν ]. In any dimension D close to 4, Note that the axial currents require renormalization if dimensional regularization is used with these conventions for the Dirac matrices.

A.3 Lattice theory
The lattice theories considered in this paper are set up as usual on a four-dimensional hypercubic lattice with spacing a = 1. For the quark fields, Wilson's formulation is chosen, where the quark fields reside on the sites of the lattice and carry the same indices as the quark fields in the continuum theory. The gauge covariant forward and backward difference operators acting on a quark field ψ(x) in presence of the gauge field U (x, µ) are defined by whereμ denotes the unit vector in direction µ. On the lattice, ∂ µ and ∂ * µ stand for the ordinary forward and backward difference operators and for their symmetric combination. The differential operators ∂ a x,µ act on differentiable functions f (U ) of the gauge field U according to While these operators depend on the choice of the generators T a , the gradient field can be shown to be basis-independent.
B Properties of the mapping (5.9) The mapping (5.9) is differentiable and satisfies P(e sX ) = is nowhere singular.

JHEP04(2013)123
Using simple norm estimates, the interior of the domain can be shown to be such a neighbourhood of unity. Moreover, the Jacobian det J(U ) is strictly positive on N. As a consequence, the identity holds for all integrable functions f (U ) and all V ∈ N, the Dirac δ-function being the one associated to the SU(N ) invariant measure on the Lie algebra. The proportionality constant k in this equation is positive and could be worked out explicitly, but its value is not needed in this paper.

C Wick contractions
The Wick contractions of the basic fermion fields in the O(a) improved lattice theory are first worked out for the case of the discretized flow equations with time step ǫ. Time-ordering ambiguities are avoided in this way and the limit ǫ → 0 can then easily be taken at the end of the calculation.

C.1 Fermion action and functional integral
To get started, it may be helpful to recall that the fermion part of the total action of the O(a) improved lattice theory in 4+1 dimensions is the sum of a boundary term, and of the bulk action (5.8). The Wilson-Dirac operator D in eq. (C.1) includes the Sheikholeslami-Wohlert term [31,32] and the mass term with bare mass matrix M 0 . In the bulk action (5.8), it is understood that the time-dependent fields satisfy the boundary conditions (2.6). The unconstrained fermion (and similarly the anti-fermion) fields integrated over in the functional integral can thus be taken to be ψ(x), χ(ǫ, x), . . . , χ(T, x), λ(0, x), . . . , λ(T − ǫ, x). (C.2) Since the action is quadratic in these fields, the fermion integral is given through Wick's theorem.

C.2 Field transformation
In the following, an important rôle is played by the kernel K ǫ (t, x; s, y) that satisfies at all t in the range 0 ≤ t ≤ T and the discretized quark flow equation x; s, y) = 0 (C.4)

JHEP04(2013)123
if 0 ≤ s ≤ t < T . Clearly, since the flow equation (C.4) can be solved step by step starting from time t = s, these equations uniquely determine the kernel at all t ≥ s. The calculation of the fermion propagators can now be simplified by performing a field transformation, where χ(t, x) is expressed through ψ(x) and another field φ(t, x), defined at t = ǫ, . . . , T , according to x; s, y)φ(s, y). (C.5) For any given field ψ(x), the mapping from the field variables χ(ǫ, x), . . . , χ(T, x) to the new variables φ(ǫ, x), . . . , φ(T, x) is invertible in view of the identity Moreover, the associated Jacobian matrix is lower triangular and its determinant does not depend on the gauge field. After the field transformation and the corresponding transformation of the anti-fermion fields, the fermion action assumes the form (C.7) The transformation thus achieves a nearly perfect decoupling of the field variables.

C.4 Contractions in the continuous time limit
When the time step ǫ is taken to zero, the kernel K ǫ (t, x; s, y) smoothly converges to the fundamental solution K(t, x; s, y) of the quark flow equation in the continuous time formulation of the lattice theory. The complete set of the non-zero contractions of the fermion fields, given by eqs. (C.8) and (C.12)−(C.16), therefore has a well-defined limit.

D Integration of the lattice flow equations
The numerical integration of the pure-gauge gradient flow has already been discussed in appendix C of ref. [1]. Here the 3rd order Runge-Kutta integrator described there is extended to the quark flow equation.

D.1 Integration of the gradient flow
Following ref. [1], the flow equation (5.2) is written in an abstract form where the gauge field V t is considered to be an element of a (high-dimensional) Lie group and the gradient Z(V t ) of the Wilson action an element of the associated Lie algebra. The Runge-Kutta integrator mentioned above proceeds in time steps of size ǫ. Assuming V t is known at some flow time t, an approximation to the exact solution V t+ǫ at time t + ǫ is obtained by computing the fields and W 3 thus provides the desired approximation to V t+ǫ .

D.2 Integration of the quark flow equation
It is again helpful to rewrite the flow equation in an abstract form, where the quark field χ t is considered to be an element of a complex vector space. Given V t and χ t at some flow time t, and assuming the gauge fields W 0 , W 1 and W 2 are as above, the quark fields φ 0 = χ t , It is not difficult to show that φ 3 = χ t+ǫ + O(ǫ 4 ), (D. 8) and φ 3 thus approximates χ t+ǫ to an accuracy that matches the one attained by the integrator of the pure-gauge flow.
E Integration of the adjoint flow equation (7.10) The numerical integration of the adjoint flow equation is complicated by the fact that the backward integration of the flow equation for the gauge field is exponentially unstable. In this appendix, the problem is first reformulated and its solution is then discussed. For simplicity, the index k in eq. (7.10) is omitted.

E.3 Improved scheme
The adjoint flow equation can be integrated more efficiently by dividing the sequence 0, ǫ, . . . , t of flow times into n b consecutive blocks of about equal size. One may then first compute the field V ǫ r at the time r at the beginning of the last block by forward integration from time 0 and keep this field in the memory of the computer. After that the adjoint flow equation is integrated from time t to r by applying the update step (E.8) and using the safe reconstruction method described above for the gauge field, where the reconstruction now starts from the stored field rather than the field at time 0.
In the next step the gauge field is calculated at the first time r in the next-to-last block and stored in memory. The integration of the adjoint flow equation can then proceed to time r as before. Clearly, the blocks can be processed in this way one after another until the integration reaches time 0.
The total number of gauge update steps in this scheme is minimized if m = n 2 b . While m will rarely be the square of an integer, one can always set n b = ⌈m 1/2 ⌉ and divide m into blocks of size m b = ⌊m/n b ⌋ and m b + 1. The computational effort then scales approximately like m 3/2 and is much smaller than the one required for the safe scheme described in the previous subsection already at low values of m.

E.4 Hierarchical scheme
A further acceleration of the computation can be achieved by dividing each block into n b smaller blocks, and these into n b even smaller blocks, and so on, with a saved gauge field at each block level. Given m and the number n s of saved gauge fields, a nearly optimal choice of the block number n b is For m ≤ 10 3 and n s = 4, for example, the total number of gauge-field updates that need to be performed is then not more than about 5m. Moreover, when several quark fields are evolved simultaneously, the intermediate gauge fields need to be computed only once.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.