The gradient flow coupling in the Schr\"odinger Functional

We study the perturbative behavior of the Yang-Mills gradient flow in the Schr\"odinger Functional, both in the continuum and on the lattice. The energy density of the flow field is used to define a running coupling at a scale given by the size of the finite volume box. From our perturbative computation we estimate the size of cutoff effects of this coupling to leading order in perturbation theory. On a set of Nf=2 gauge field ensembles in a physical volume of L ~ 0.4 fm we finally demonstrate the suitability of the coupling for a precise continuum limit due to modest cutoff effects and high statistical precision.


Introduction
Finite-volume renormalization schemes have now a long history in lattice field theory (see [1][2][3] or the pedagogical reviews [4,5]). Asymptotic freedom tells us that at small distances QCD is well described by perturbation theory, while at large scales QCD is a strongly interacting theory. Instead of trying to accommodate these two scales in a single lattice simulation, the idea of finite-size scaling exploits the size of a finite volume world as renormalization scale. A single lattice simulation can resolve only a limited range of scales, but one can match different lattices and adopt a recursive procedure to cover a large range of scales. In this way one can connect the perturbative and non-perturbative regimes of QCD.
Beside a successful application of the finite-size scaling technique in the case of pure Yang-Mills theory [2,3], the running coupling [6][7][8][9][10] and quark mass [11][12][13] have been computed non-perturbatively in QCD with different flavour content by ALPHA and other collaborations. Since the general idea of finite-size scaling is a very powerful tool to solve scale dependent renormalization problems, it is not surprising that it is broadly used also in other strongly interacting theories, even in effective theories such as HQET [14][15][16]. There has been a growing interest in other than QCD strongly interacting gauge theories, especially in connection with electroweak symmetry breaking and quasi-conformal behavior (see for example [17] and references therein). Finite-size scaling techniques are also a powerful tool to study these systems.
Basically there are two things that are needed to perform the previously sketched program. First one needs to define exactly what is meant by a finite-volume scheme, i.e., one has to specify the boundary conditions of the fields. Second, one needs a nonperturbative definition of the coupling. In principle there are many valid possibilities, but practical considerations have to be taken into account. Good options should allow for an easy evaluation of the coupling constant both in perturbation theory and in a numerical, non-perturbative (lattice) simulation.
The rest of this section is mainly dedicated to explain why we choose the Schrödinger functional (SF) scheme [2] as our finite-volume setup and the Wilson flow for a nonperturbative definition of the coupling [18]. To simplify the following discussion we will argue about a pure SU (N ) gauge theory in 4-dimensional Euclidean space-time.
In the Schrödinger functional [2,19] one embeds the fields in a finite volume box of dimensions L 3 × T . Gauge fields in the SF are periodic in the three spatial directions and have Dirichlet boundary conditions in time direction (i.e. one fixes the value of the gauge fields at x 0 = 0, T ). The value of the gauge fields at the time boundaries are called boundary fields. One can interpret the partition function of the theory as the transition amplitude of the gauge field to propagate from the boundary value at x 0 = 0 to the boundary value at x 0 = T . Such a setup has nice properties in perturbation theory. In particular with a smart choice of the boundary fields one can guarantee that there is a unique gauge field configuration (up to gauge transformations) that is a global minimum of the action. This avoids some difficulties with perturbation theory [20][21][22]. The reader interested in this issue will appreciate the original literature, as well as the nice discussion in [23].
Recently, the gradient flow has been used in different contexts [24][25][26], but it is the proposal made in [18] to define a renormalized coupling through the gradient flow in nonabelian gauge theories what inspires this work. The gradient flow defines a family of gauge fields parametrized by a continuous flow time t. The flow equation brings the gauge field towards the minimum of the Yang-Mills action, and therefore represents a smoothing process. The key point is that correlation functions of the smoothed gauge field defined at t > 0 are automatically finite [27]. One can use the expectation value of the energy density, where G µν (t) is the field strength of the gauge field at flow time t, to give a non-perturbative definition of the gauge coupling. This idea was applied to set the scale in lattice simulations [18,28], to tune anisotropic lattices [29] and more recently in a similar context of this work (finite-size scaling, but using a box with periodic boundary conditions) to compute the step scaling function in SU (3) with four fermion species [30].
In this paper we investigate the perturbative behavior of the Wilson flow in the Schrödinger functional. This motivates us to propose a gradient flow coupling with a normalization factor N to be determined later, valid for an arbitrary SU (N ) gauge field coupled (or not) to fermions. Relating t and L the coupling depends only on one scale, the size of the finite volume box, and therefore can be used for a finite-size scaling procedure in the same way as the traditional SF coupling. The paper is organized as follows: in the next section we investigate the perturbative behavior of E(t) in the SF, both in the continuum and on the lattice. Section 3 uses this information to define the gradient flow coupling in the SF, and to discuss some practical issues: cutoff effects, boundary fields and fermions. In section 4 we investigate this coupling numerically on a set of lattices in a physical volume of L ∼ 0.4 fm and finally conclude in section 5. Details needed for the computation have been summarized in form of appendices: a summary with some useful notation A, heat kernels B, propagators in the SF C and finally some practical details on how to integrate the Wilson flow in numerical simulations D.

Perturbative behavior of the Wilson flow in the SF
We would like to start this section by recalling the original proposal of using the Wilson flow and the energy density as a definition for a coupling in gauge theories [18]. Later it will become clear what role the SF setup plays.

Generalities
By considering the gauge fields to be functions of an extra flow time t, not to be confused with Euclidean time, denoted x 0 , the Wilson flow is defined by the non-linear equation δBµ gauge fields along the flow become smoother, eventually reaching a local minimum of the Yang Mills action: the flow smooths the fields over a region of radius √ 8t. The somewhat surprising result of [18,27] is that correlation functions made of this smoothed field have a well-defined continuum limit. In particular the energy density in SU (N ) Yang-Mills theory in infinite volume has the perturbative behavior At a scale µ = 1/ √ 8t, c 1 is a numerical constant and g MS (µ) is the renormalized coupling in the MS scheme. Therefore one can define a running coupling constant α(µ) from These expressions are valid in infinite volume. What about the Schrödinger Functional?
The computation is completely analogous, but we have to impose the correct boundary conditions to the gauge fields. As we have mentioned in the SF gauge fields are restricted to a box of dimensions L 3 × T . They are periodic in the three spatial directions and the spatial components have Dirichlet boundary conditions at x 0 = 0 and x 0 = T . We are going to work exclusively with zero boundary fields, which means The flow equation (2.1) has to be solved maintaining these boundary conditions at all flow times t. To apply the idea of finite-size scaling, as has previously been done in [23] in a periodic box, one simply has to run the renormalization scale with the size of the finite volume box given by L via Here c is a dimensionless constant that represents the fraction of the smoothing range over the total size of the box. In this way the flow coupling will not depend on any scale other than L. The renormalization scheme will depend on the values of c, ρ = T /L and 1 x 0 /T 8) where N −1 (c, ρ, x 0 /T ) will be computed in the next section in order to ensure

Continuum
Our computation follows the lines of [27]. First we consider the modified flow equation One can transform a solution of the last equation into a solution of the canonical flow equation (2.1) (corresponding to α = 0) by a flow-time dependent gauge transformation. In particular, if B µ is a solution of (2.10) one can construct a solution of (2.1) via as long as Λ obeys the equation This shows that gauge-invariant quantities are independent of α. For instance, setting α = 1 turns out to be a very convenient choice for perturbative computations. Due to the periodicity in the spatial directions it is natural to expand the gauge fields as As already mentioned, in the SF the gauge field is periodic in the three spatial directions and its spatial components have Dirichlet boundary conditions in time, eq. (2.5) and (2.6) respectively. On the other hand the boundary conditions of the time component of the gauge field are not fixed but naturally emerge through the gauge fixing condition. 2 To properly derive the boundary conditions for B 0 it is convenient to work in the lattice formulation and derive the boundary conditions by taking the continuum limit. We will postpone this derivation to the next section and simply state the result here: B 0 obeys Neumann boundary conditions at non-vanishing spatial momentum, while for zero momentum B 0 obeys mixed boundary conditions. Thus in the present set-up the full set of boundary conditions reads The modified Wilson flow equation with α = 1 is given by After rescaling the gauge potential with the bare coupling A µ → g 0 A µ , the flow becomes a function of the couplingB Inserting this expression in the modified flow equation, we find that to leading order in g 0 the flow equation is just the heat equation with initial condition A µ : i.e., to leading order the Wilson flow is the heat flow. We also observe that different momentum modes do not couple to each other at this order. Together with the fact that the zero momentum mode B 0 (0, x 0 , t) does not contribute to the observable of interest, E(t) = 1 4 G µν G µν , we can safely neglect the special treatment that the boundary conditions of the zero momentum mode B 0 (0, x 0 , t) would otherwise require in the following discussion.
We have to solve the heat equation respecting the boundary conditions (2.14). This is easily done by using appropriate heat kernels Since the boundary conditions of the fieldB µ,1 (p, x 0 , t) are inherited from the boundary conditions of the heat kernels, we have to choose them with the correct boundary conditions. Heat kernels with either Dirichlet (K D (x, x , t)) or Neumann (K N (x, x , t)) boundary conditions can be constructed from the basic periodic (K P (x, x , t)) heat kernel in [0, L] given by Explicit expressions are given in appendix B. Our observable, the energy density E(t, x 0 ) , has an expansion in powers of g 0 . The leading contribution is given by We are going to split the computation in two parts, one involving only the spatial components of G µν , and the other involving the mixed time-space components of G µν The final result is obtained inserting the SF gluon propagator [31,32]. Since our observable is invariant under gauge transformations of the A µ (x) field we will use the Feynman gauge, where the expression for the gluon propagator turns out to be more easy (for additional details see appendix C) 3 .
To shorten notation we use After some algebraic work one arrives at the expression and a very similar computation leads to

Lattice
On the lattice one defines the Wilson flow as is an arbitrary function of the link variable U µ (x), the components of its Liealgebra valued derivative ∂ a x,µ are defined as In a neighborhood of the classical vacuum configuration the lattice fields U µ (x) and V µ (x, t) are parametrized as follows:

Gauge fixing
To simplify our perturbative computations it is useful to study a modified equation with a gauge damping term. It is easy to check that the lattice flow equation (2.29) is invariant under flow-time independent gauge transformations. On the other hand one can consider the modified equation with V Λ µ (x, 0) = U µ (x) and the forward lattice covariant derivativeD Λ µ acting on Lie-algebra valued functions according tô With∂,∂ * we denote the forward/backward finite differences respectively as defined in appendix A.
The solutions of the modified equation (2.32) and the original flow equation (2.29) are related by a gauge transformation and therefore one can freely choose the function Λ(x, t). To fix the gauge the most natural choice is to use the same functional that is used for the conventional gauge fixing. As is detailed in appendix C, we choose Note Λ(x, t) does not depend on x at x 0 = 0 and x 0 = T , as a decent gauge transformation should be in the Schrödinger functional according to our conventions (see appendix C for details). We observe that on the lattice the time component of the gauge field B 0 (x, t) is completely free and does not obey any particular boundary conditions. To understand how the boundary conditions for B 0 (x, t) arise in the continuum theory, one can extend the domain of definition of B 0 (x, t) to −a ≤ x 0 ≤ T and choose to fix the additional variables with the condition∂ * This equation can be interpreted as a boundary condition for the B 0 (x, t) field. In particular B 0 (x, t) has Neumann boundary conditions at x 0 = 0, T , except for its spatial momentum zero mode that has a mixture of Neumann boundary conditions at x 0 = T and Dirichlet boundary conditions at x 0 = −a.
This justifies our previous choice of boundary conditions in the continuum, eq. (2.14). With this useful convention in mind eq. (2.35) simply reads

Behaviour of E(t) in lattice perturbation theory
We again note that the value of any gauge invariant observable is independent of our choice of α in equation (2.39). In particular, with the choice α = 1 the modified flow equation reads 40) and to first order in g 0 Using periodicity in the spatial directions, we expand and the flow equation becomes wherep is the usual spatial lattice momentum, see appendix A. Now we have to solve a special type of heat equation in which the Laplacian is substituted by a discrete version, but the flow time remains a continuous variable. The strategy is very similar: We find the fundamental solutions of this equation, i.e., the discrete heat kernels given in appendix B, and writẽ Then we have to insert this in our lattice observable E . We use the clover definition for G µν that to leading order in g 0 reads . The computation is completed by using the lattice gluon propagator For the spatial part of the contribution to the energy density we arrive at while for the mixed part we obtain (2.50) The definitions of the lattice momentap,p,p and the functionsŝ p (x),ĉ p (x) are summarized in appendix A.

Tests
There are several tests that can be performed to check the previous computations. Let us first concentrate on the continuum computation. At fixed t, the infinite volume limit L → ∞ (with ρ kept constant) is taken through c → 0. For this case the continuum expression transforms into an integral via and we obtain thus recovering the infinite volume result of [27].
Another rather obvious check is that one should recover the continuum result from the lattice expression in the limit a/L → 0. This can be easily checked by noting that the sums (2.49) and (2.50) are dominated by terms with small lattice momenta ap µ → 0.
Finally we have performed some simulations with the openQCD code [33] at small values of the bare coupling in a pure SU (3) gauge theory. Using a 8 3 × 7 lattice and varying the value of the bare coupling (β = 60, 120, 240, 360, 480, 600, 720, 840, 960, 1080, 1200) we compare the analytical lattice prediction and the numerical results after collecting 10000 measurements of the gradient flow coupling for each value of β. We use the clover definition for G µν to compute the value of The lattice computation of can be checked in the following way: plotting All intercepts are of the order 10 −4 and compatible with zero within errors. The difference inÊ 0 (t, x 0 ) for different values of c, x 0 or between the continuum and the lattice result varies between 5% and 10%. We note that this last test is highly non-trivial since it is done for arbitrary x 0 at ρ = 1 on a small lattice where cutoff effects tend to be larger.

Definition of the flow coupling
Using our continuum result we define the gradient flow coupling for non-abelian gauge theories in the SF by means of This definition of the coupling is valid if the gauge field is coupled to fermions in arbitrary representations. As the reader may have noticed the scheme that defines the coupling depends not only on the quantities c, ρ, x 0 , but also on the value of the fermionic phase angle and the background field. In the simulations of the Schrödinger functional it is customary to include a phase angle θ in the fermionic spatial boundary conditions. In principle different values of θ are different schemes, although we have observed in some practical situations that the difference of the gradient flow coupling between θ = 0 and θ = 0.5 is below the 2%.
Up to now we have worked exclusively with zero background fields, but the generalization to other values is straightforward. It only requires the modification of the heat kernels to preserve the value of the boundary fields and a modified form of the propagator [32]. Nevertheless common wisdom suggests that cutoff effects are reduced for zero background field, therefore we prefer to work in this scheme. In this case the definition of the coupling is also symmetric about x 0 = T /2 and we choose that value to minimize boundary effects. Also choosing ρ = T /L = 1 seems reasonable and leaves us with a one-parameter family of couplings, parametrized by the smoothing ratio c.
By comparing the lattice and continuum behavior of the energy density as a function of c we can compute the leading order size of cutoff effects in the gradient flow coupling. As the reader can see in figure 2, the cutoff effects are large for small values of c, reach a minimum around c ∼ 0.5 and then grow again. We recall that with c = 0.5 the smoothing radius is equal to L/2, and therefore one is effectively smoothing over all the lattice. For c = 0.3 cutoff effects are smaller than 10% for a lattice of size L/a ≥ 8, while for c = 0.4 even the L/a = 6 lattice has cutoff effects of about 10%.
This figure suggests using c = 0.5 as a preferred scheme, but later, when lattice simulations enter into the game, we will see that the statistical errors of the coupling also grows with c, and therefore in practice it is better to stay with c ∈ [0.3, 0.5], probably depending on the particular case, but this is the subject of the next section.
We would also like to comment that if one is performing numerical simulations with the Wilson gauge action, one can benefit from smaller cutoff effects by using the lattice the coupling is given by Obviously both definitions of the coupling differ only by cutoff effects. We finally want to mention that it is possible to define analogous couplings by using only the spatial components G ik (t)G ik (t) . In a lattice simulation one stays further away from the boundaries by not including plaquettes with links in the time direction. This may result in smaller cutoff effects, although this point needs further investigations.

Non-perturbative tests
In this section we would like to analyze the gradient flow coupling numerically. We want to estimate both the size of cutoff effects and the numerical cost of evaluating the new gradient flow coupling. The main result of this section is that both quantities depend on the particular scheme via the parameter c. When c is increased cutoff effects decrease, but the numerical cost increases. We find that the window of values c ∈ [0.3, 0.5] allows a very precise determination with a mild continuum extrapolation.

Line of constant physics
As framework for our tests we choose a set of N f = 2 Schrödinger functional simulations at a line of constant physics as given through where g 2 SF is the traditional SF coupling and m the renormalized PCAC mass. From reference [34] we know that the physical volume is roughly L 1 ∼ 0.4 fm. The available five different ensembles are lattices with L/a ∈ {6, 8, 10, 12, 16} at T = L with vanishing boundary gauge fields and a fermionic phase angle θ = 0.5. Each ensemble consists of at least 8000 configurations separated by τ meas = 10 molecular dynamic units (MDU). We refer the reader to the appendices in [35] for any unexplained detail concerning the physics and run parameters.
We would like to measure the value of the gradient flow coupling in these ensembles. Since they have been tuned to have constant SF coupling, equivalent to constant volume in the continuum, we define the function Ω(u; c, a/L) = N −1 (c, 1, 1/2, a/L) · t 2 E(t, T /2) u=4.484, msea=0, θ=0.5 that at fixed c has the gradient flow coupling as continuum limit  We will also use some ensembles with larger lattices (L 1 /a = 24, 32, 40) but lower statistics. These have been defined by a slightly different line of constant physics: and thus differ only by cutoff effects. The statistics for this second set of ensembles is smaller (∼ 800 measurements). For these additional lattices we define a function similar to Ω according to Ω(u; c, a/L) = N −1 (c, 1, 1/2, a/L) · t 2 E(t, T /2) u=σ(ũ), msea=0, θ=0.5 that also has the same continuum limit. All ensembles have been tuned to have constant SF coupling only with some statistical accuracy. This propagates into an uncertainty in the determination of the functions Ω(u; c, a/L) and Ω(u; c, a/L). This error can be estimated by simple propagation of errors, for example applied to Ω(u; c, a/L) we have δΩ(u; c, a/L) = ∂Ω ∂u δu . (4.7) To evaluate this uncertainty we use another ensemble (labeled 12 * in table 2) with a slightly different value of β but also tuned to have vanishing quark mass. By evaluating both the SF coupling and Ω on this ensemble we can numerically estimate the derivative in equation (4.7) 4 . This source of error in fact dominates the error budget of Ω(u; c, a/L), which anticipates that the new coupling is numerically more precise. Figure 3 shows the gradient flow coupling as a function of c for the different ensembles. For c ∈ [0.3, 0.5] we observe a monotonic behavior of g 2 GF with a/L. As figure 4 shows this seems to be the scaling region of the gradient flow coupling for lattices with L/a > 8, and therefore this is the region on which we will focus from now on.

Numerical results and computing cost
In figure 5 we present the noise-to-signal ratio as obtained in our analysis as function of c for the individual lattices. We observe that the noise-to-signal ratio increases with increasing c (see table 2). Although our statistics does not allow us to draw definite conclusions, we observe a behavior compatible with a powerlike scaling of R NS with c. The behavior seems to be universal and independent of a/L in contrast to the traditional SF coupling, that has a divergent variance when approaching the continuum [37]. The product √ N meas R NS can directly be translated in the cost of obtaining the new gradient flow coupling with some precision. 8000 measurements are enough to achieve a precision of 0.1% for c = 0.3, while for c = 0.5 the precision decreases to 0.35%.
In figure 6 we plot the integrated auto-correlation time τ int for the different ensembles. The lattices L/a = 6, 8 have τ int ∼ 0.5 which means that the available configurations are too far separated in Monte Carlo simulation time to detect any auto-correlations in the chain. The L/a = 10, 12, 16 lattices show a clear increase of auto-correlations with increasing L/a. As can be inferred from figure 6 this increase is compatible with a scaling L/a  (9) 3.6(2) - Table 2: Lattice run parameters (see [35] for more details) and results for the gradient flow coupling at smoothing ratio c ∈ {0.3, 0.4, 0.5}. Errors are computed using the Γmethod [36]. We show the values of the gradient flow coupling g 2 GF (L 1 ) and of the function Ω(u; c, a/L). Furthermore, we quote the integrated auto-correlation time τ int of g 2 GF (L 1 ), estimated in units of the measurement frequency τ meas = 10 MDU (which is the same for each lattice), and the noise-to-signal ratio R NS . The lattice labeled 12* is used to estimate (∂Ω/∂u).  ∼ 1/a 2 towards the continuum. This is much in accordance with the conjecture of [38,39] for the scaling behaviour of the HMC algorithm in an interacting theory, since we do not expect non-zero topological charge sectors to contribute significantly at this small physical   volume.

Cutoff effects
For the continuum approach of the function Ω(u; c, a/L) (and also of Ω(u; c, a/L)) we observe a behaviour dominated by a linear scaling in (a/L) 2 . Hence, we choose Ω(u; c, a/L) = ω(u; c) 1 + A(u; c) · (a/L) 2 (4.9) as fit ansatz to extract the continuum limit ω(u; c) = g 2 GF (L) and the leading cutoff effects A(u; c). Figure 7 shows examples of these extrapolations for three representative cases c = {0.3, 0.4, 0.5}. We observe that cutoff effects decrease with increasing c. Quantitatively that decreases by a factor 2 when c is increased from 0.3 to 0.5 (see figure 7). Since in general cutoff effects of the SF step scaling function are very small, and given the fact that cutoff effects of Ω(u; c, a/L) change with c, we think that the cutoff effects in Ω(u; c, a/L) and Ω(u; c, a/L) are dominated by the lattice spacing dependence of the new gradient flow coupling. Therefore we expect both functions to show roughly the same scaling behavior. Although the points corresponding to Ω(u; c, a/L) have not been used for the previous continuum extrapolations, we have added the points to the plots in figure 7. The data fits well into the expected scaling behavior.
Summarizing figure 7 and table 2 we can say that when c is increased from c = 0.3 to c = 0.5 the relative cutoff effects decrease by about a factor of 2. We remember from previous sections that this decrease of cutoff effects comes at the expense of an increased relative statistical error (about three times larger when going from c = 0.3 to c = 0.5).

Mass dependence
Last but not least we have studied the mass dependence of the gradient flow coupling g 2 GF as it was done for g 2 SF in [35]. For this purpose we have generated an ensemble with L/a = 8 at the same value of the bare coupling as the available one, but with a non-zero quark mass. Actually the bare parameters of the simulation correspond to the lattice labeled as 8 * in [35], and the interested reader is encouraged to consult the original work for more details. Defining the dimensionless PCAC quark mass z = Lm, we obtain to be compared to the corresponding value of 1.4(4) for the Schrödinger functional coupling. The mass dependence of the gradient flow coupling as defined in the present paper is smaller by an order of magnitude.

Conclusions
The gradient flow can be used to define a renormalized coupling at a scale µ = 1/ √ 8t. In this work we have studied the perturbative behavior of the gradient flow in the Schrödinger functional. By setting the renormalization scale proportional to the linear size of the SF box, µ = 1/ √ 8t = 1/cL, we have defined a family of running coupling constants valid for an arbitrary SU (N ) gauge field coupled to arbitrary fermions. Since this coupling definition does not depend on any scale but the finite volume, it can be used for finite-size scaling purposes.
The coupling constant can be defined for different values of the background field in the SF. Since one expects cutoff effects to be smaller for the case of vanishing background field this is the case that we have studied in more detail. From our perturbative analysis we have been able to study the size of cutoff effects to leading order. Cutoff effects tend to be relatively large for either small or large values of c, but very mild for c ∈ [0.3, 0.5].
As an example for c = 0.3 the difference between the coupling in a L/a = 8 lattice and the continuum is around 10%, while for c = 0.5 the difference between the continuum and the value in an L/a = 6 lattice is below 4%.
We have analyzed a total of five ensembles tuned to have a constant SF coupling in a physical volume of L ∼ 0.4 fm. The cutoff effects observed in the gradient flow coupling on these ensembles shows a similar overall behavior as expected from earlier perturbative considerations. Provided that the smoothing ratio c is choosen wisely we can conclude that the gradient flow coupling has mild cutoff effects even for small lattices.
We have analyzed the numerical cost of evaluating the gradient flow coupling, and see that it increases with c. Nevertheless it can be computed very precisely with a modest numerical effort. For the case studied in detail, a box of size L ∼ 0.4 fm, we find that roughly 8000 measurements are enough to obtain a precision of 0.1% for c = 0.3, even on our larger lattice L/a = 16. The worst analyzed case at, c = 0.5, this precision drops to the still very good figure of about 0.35%.
What is the most convenient scheme? From a practical point of view there is no need to settle this discussion -and thus the unique definition of the scheme -immediately since one measures the gradient flow coupling at different values of c in any case while integrating the Wilson flow.
In summary we conclude that the gradient flow coupling in the SF has several practical advantages. It is valid for arbitrary SU (N ) gauge fields coupled to fermions in any representation, having the universal two loop beta function. It has mild cutoff effects and can be cheaply evaluated numerically with high precision. It has a very small dependence on quark masses and is naturally defined with vanishing background field, avoiding the generation of new ensembles just to compute the running coupling.
We think that the computation of the next order in the perturbative behavior of the gradient flow coupling in the SF is very interesting, in particular in connection with the determination of the Lambda parameter. This computation can also shed some light on an optimal value of c for the matching between different schemes. There are many applications for this new coupling, but we are particularly interested in using it to compute the step scaling function of QCD. We believe that this can be achieved with a very high accuracy with a modest computational effort.

Acknowledgments
In particular, the authors would like to thank R. Sommer for suggesting us (with immense patience) to look into this very interesting problem and also for guiding us in the steps of this work.
We have profited from the generosity of many individuals who have provided us with code for several purposes. We thank M. Lüscher for making available a code to integrate the flow equations numerically, to H. Simma for the invaluable help in adapting the code to our needs, and finally to M. Marinkovic (SF-MP-HMC) and the openQCD team (M. Lüscher, S. Schaefer and J. Bulava for implementing SF boundary conditions) for their very useful codes that we have used at different stages of this project to generate new ensembles.
A. R. also wants to thank A. Gonzalez-Arroyo, S. Sint, U. Wolff and F. Virotta for the many illuminating discussions concerning the role of the boundary conditions, the SF and the Wilson flow. Thanks, guys! This work is supported by the Deutsche Forschungsgemeinschaft in the SFB/TR 09 and the John von Neumann institute for computing. Our numerical computations were performed on the PAX cluster at DESY, Zeuthen.

A Notation
Momenta are always defined to be for the periodic, spatial directions and in the time direction (x 0 ). In the continuum sum over momenta are abbreviated by while in the lattice we have finite sums If we impose Dirichlet/Neumann boundary conditions on the interval the Laplacian has eigenfunctions and eigenvalues given by respectively. The corresponding eigenvalues are given through and these functions obey the completeness relations On the lattice derivatives are substituted by finite differenceŝ and defining the family of lattice momentâ Since the heat equation is a linear equation one can construct heat kernels with different boundary conditions with a linear combinations of the heat kernel in Euclidean space. In particular is, by construction, both a solution to the heat equation and is periodic in x with period L. Therefore it corresponds to the heat kernel on S 1 . Following a similar reasoning it is easy to see that (B.5) are heat kernels with Dirichlet and Neumann boundary conditions on the interval [0, T ] respectively. It is worth mentioning that, via the Jacobi imaginary transformation, the periodic heat kernel is nothing more than the third Jacobi theta function 5 This last expression, that can also be obtained using Poisson summation formula, turns out to be convenient for our computations.

B.2 Discrete heat kernels
When performing computations of the Wilson flow on the lattice we find a type of heat equation in which the time variable is continuous but the Laplacian is substituted by a discrete version Taking appropriate boundary conditions into account, we call fundamental solutions of this equation discrete heat kernels. The most easy way to construct them is by noting that one can formally writeK = exp t∂ x∂ * x (B.8) and now the task of finding the heat kernels is reduced to finding eigenvalues/eigenfunctions of the discrete Laplacian with the correct boundary conditions. By recalling the notions of 5 A standard reference with the same conventions used here is [41] appendix A one can immediately writê

B.3 Properties
The following properties are straightforward, but fundamental to easily reproduce our results: Finally we note that the construction of heat kernels in more than one dimension is done simply by multiplying one-dimensional heat kernels. For examplê is a 4-dimensional heat kernel with periodic b.c. in the three space dimensions and Dirichlet b.c. in the time dimension.
C Gauge fixing and gluon propagator in the SF

C.1 Gauge fixing and boundary conditions
Here we will deal with gauge fixing in the lattice formulation of the SF. Basically we adapt the contents of section 6 of [2] to our specific problem of zero background field. In the SF admissible gauge transformations must leave the boundary conditions of the gauge field invariant. In our particular case of zero background field this means that The lattice forward derivative∂ µ maps any infinitesimal gauge transformation w(x) to the gauge field∂ µ w(x). As the reader can check the boundary conditions for w(x) imply the correct boundary conditions for the gauge field∂ µ w(x). The corresponding operator is denoted by d. Its adjoint maps any gauge field into an infinitesimal gauge transformation and is defined by (d * A, w) = −(A, dw) . (C.7) The explicit transformation done by d * is This is precisely the gauge fixing function that we have to add to the action. In explicit form it is given by (C.9) To understand how the boundary conditions of A 0 (x) arise in the continuum, one can extend the domain of definition of A 0 (x) to −a ≤ x 0 ≤ T and choose to fix the additional variables with the condition Now the gauge fixing function is simply given by Note that the additional variables are not dynamical, but completely determined by the previous condition. The trick is purely a matter of convention, but now equation (C.10) can be interpreted as a boundary condition for A 0 (x). Our gauge fields have therefore Neumann boundary conditions at x 0 = 0, T , except for the zero momentum mode that has mixed boundary conditions.

(C.22)
Note that in the continuum limit D 00 (p, x 0 , y 0 ) obeys Neumann boundary conditions for p = 0, while the spatial zero momentum mode has mixed boundary conditions.

D Adaptive size integrators for the Wilson flow
On the lattice the flow equation has the form Following the advice in [27], we use the third order Runge-Kutta scheme given by We simply want to point out that with one extra computation, one can get a second estimate of V µ (x, t + a 2 ) This distance can be used to estimate the error made by the lower order integrator d = max x,µ dist(V µ (x, t + a 2 ), V µ (x, t + a 2 )) . (D.9) With this information one can adjust the step size so that the error in the integration does not exceed certain tolerance δ. After each integration step is updated according to −→ 0.95 3 δ d (D. 10) and if d > δ the integration step is repeated.
The reason why this scheme is efficient for the particular case of the Wilson flow is related to its smoothing properties. Close to t = 0 the configuration is rough and therefore a very fine integration is needed. But as t increases the configuration is more and more smooth, and one can have a very precise integration with a large . In practical cases it has saved us around a factor 4 in the number of integration steps. It also has the advantage that is tunned automatically, and one does not need to worry about the step size, but only plug in the desired tolerance.