The gradient flow running coupling with twisted boundary conditions

We study the gradient flow for Yang-Mills theories with twisted boundary conditions. The perturbative behavior of the energy density $\langle E(t)\rangle$ is used to define a running coupling at a scale given by the linear size of the finite volume box. We compute the non-perturbative running of the pure gauge $SU(2)$ coupling constant and conclude that the technique is well suited for further applications due to the relatively mild cutoff effects of the step scaling function and the high numerical precision that can be achieved in lattice simulations. We also comment on the inclusion of matter fields.


Introduction
Non-abelian Yang-Mills theories play a central role in our understanding of high energy physics. Despite their apparent simplicity (they depend on one free parameter, the coupling constant), they have proven to produce very rich phenomena at the quantum level.
The invariance under conformal transformations present in the classical theory is broken by quantum corrections and a typical scale of the interactions (usually characterized by the Λ-parameter) comes into play. The coupling constant becomes scale dependent and asymptotic freedom [1,2] guarantees that at high energies its value is small. A perturbative analysis has proven to be very successful in making predictions in this regime, but at low energies the coupling becomes large and non-perturbative effects set in. Confinement characterizes the interaction in this regime and perturbation theory provides little information.
Finite volume renormalization schemes together with non-perturbative numerical simulations [3][4][5] play a major role in understanding how and when this transition from the perturbative to the non-perturbative regimes of YM theories happens. By identifying the linear size of a finite volume box with the renormalization scale a controlled and nonperturbative running of the coupling constant can be carried over from the very perturbative regime to the non-perturbative one.
There are many practical issues that have made the Schrödinger Functional [4] (SF) the preferred scheme to perform the previously mentioned program (see [6][7][8][9] for some applications). In the SF one imposes Dirichlet boundary conditions on the spatial components of the gauge fields at x 4 = 0, T , while they remain periodic in the spatial directions. By choosing appropriately the values of the fields at the time boundaries one can define a running coupling, usually denoted g 2 SF , with many good properties from a practical point of view. It is easy to evaluate numerically and precise, especially in the perturbative regime.
Other schemes have been proposed: In the twisted Polyakov loop scheme (TPL) the gauge fields are embedded in a torus with twisted boundary conditions and the ratio of Poliakov loops between the twisted and periodic directions is used to define a running coupling [5]. This scheme has some good properties, related to the manifest invariance under translations. In particular O(a) improvement in the pure gauge case is guaranteed, and therefore one expects smaller cutoff effects. On the other hand the observable used to define the coupling tends to be more noisy than the SF coupling.
When one wants to study QCD-like theories, or other models with fermions coupled to the gauge field, similar considerations apply. In the SF scheme [10,11] one can couple an arbitrary number of fermions in any representation to the gauge field. The coupling constant is defined in a similar way, and maintains its good properties. Fermions induce additional boundary counterterms that one needs to compute to have O(a 2 ) scaling. Moreover the boundary conditions for the fermion fields typically break chiral symmetry, and therefore one also needs bulk improvement, although this last issue can be addressed by modifying the boundary conditions of the fermion fields [12]. On the other hand, the TPL scheme can not be used with an arbitrary number of fermions in the fundamental representation. The twisted boundary conditions put a constraint on the number of fermions that can be coupled to the gauge field. But when this scheme can be used, it guarantees a better scaling towards the continuum. O(a) improvement is automatic, even with Wilson fermions, provided that one works with massless quarks [12,13]. Nevertheless the observable used to define the coupling (a ratio of polyakov loops), tends to be more noisy than the SF coupling, especially in the non-perturbative domain.
Recently the nice properties of the gradient flow regarding the renormalization of composite operators have introduced new observables as candidates for a running coupling definition [14,15]. By introducing an extra parameter (the flow time t) one defines a one parameter family of gauge fields. The evolution of the gauge field in flow time is given by a diffusion-like equation that drives the gauge field towards a classical solution of the YM equations, and therefore is a smoothing process. The key point is that the flow field does not require any renormalization at positive flow time (t > 0), and correlation functions of the smooth field can be considered as observables for a running coupling definition.
In particular the energy density of the flow field where G µν (t) is the field strength of the gauge field at flow time t, can be used to give a non-perturbative definition of the coupling at a scale µ = 1/ √ 8t. Moreover one can identify this renormalization scale with the linear size of the box and define a running coupling.
This idea was first applied to the case of a periodic box [16]. The dynamics of Yang-Mills theories in a periodic box contains contributions from zero momentum modes that are not Gaussian and therefore have to be treated exactly [17] (see also the review [18] and references therein). Making a long story short, this leads to a definition of the running coupling that is non-analytic in g 2 MS . The author wants to stress that there is nothing wrong with such a coupling definition, but often one wants to make contact with perturbation theory (i.e. when determinning the Λ parameter). A non-analytic definition of the coupling may make contact with perturbation theory at a larger energy scale and looses many nice properties, like a universal 2-loop beta function. Moreover perturbative computations in this setup are usually more involved and difficult because one has to treat the zero mode (toron) contribution non perturbatively.
These difficulties can be avoided with a different choice of boundary conditions. For example in the SF, if the boundary values are chosen wisely, zero momentum modes become incompatible with the boundary conditions, and therefore these non Gaussian modes are absent of the dynamics. A definition of a running coupling using the gradient flow in this setup has been proposed in [19]. It leads to a coupling definition analytic in g 2 MS and with a universal two loop beta function.
The infamous problem of topology freezing that affects large volume simulations [20,21], has recently been shown to also affect step scaling studies [22,23]. To overcome this problem it has been proposed [23] to use a setup with mixed boundary conditions. Since in this scheme the fields satisfy open boundary conditions at x 0 = T , topological freezing is avoided [23,24]. On the other hand the fields satisfy Dirichlet boundary conditions at x 0 = 0, and the absence of zero momentum modes and the analyticity of the observable used to define the coupling is guaranteed.
In this paper we propose another alternative based on using twisted boundary conditions for the gauge fields as in the TPL scheme. These twisted boundary conditions guarantee that the action has a unique minimum up to gauge transformations, and therefore zero momentum modes are not present in the dynamics. Moreover the formulation is manifestly translation invariant, since the scheme is defined in a torus, guaranteeing the absence of O(a) effects, even when one works with massless Wilson fermions. This paper is organized as follows. Section 2 presents a small review of the twisted boundary conditions. Readers interested in more details are encouraged to look at the review [25] and the recent paper [26]. Section 3 studies the perturbative behavior of the gradient flow with twisted boundary conditions. Again we encourage the reader to consult the original works [14,15] for a better understanding of the properties of the gradient flow. Section 4 uses the previous results to define a non-perturbative coupling that runs with the size of the box. Finally in section 5 we compute the running coupling for the case of an SU (2) pure gauge theory, and conclude that the modest size of cutoff effects and small variance of the observable when using Monte Carlo techniques to compute it, make it an interesting choice for further studies.

Twisted boundary conditions
The twisted boundary conditions were first introduced by 't Hooft [27] to characterize confinement. But for us, they will be used simply as a tool to study the renormalization of Yang-Mills theories. The main observation is that the requirement for physical quantities to be periodic can be accomplished by fields that change by a gauge transformation under translations over a period.

Gauge fields
We will consider SU (N ) gauge fields in a four dimensional torus of size L 4 . The twisted boundary conditions are implemented by requiring the field to gauge-transform under the displacement of a period where Ω µ (x) are known as the twist matrices. The uniqueness of A µ (x + Lμ + Lν) requires that the twist matrices have to obey the relation where n µν is an anti-symmetric tensor of integers modulo N called the twist tensor. It is easy to check that under a gauge transformation, Λ(x), the twist matrices change according to but the twist tensor n µν remains unchanged. Therefore all the physics of the twisted boundary conditions is contained in the twist tensor, and the particular choice of twist matrices is irrelevant. One can restrict the gauge transformations to those that leave the twist matrices unchanged. It is easy to check that the necessary and sufficient condition for the gauge transformations is to obey the periodicity condition The reader interested in knowing more about the twisted boundary conditions is invited to consult the review [25]. Here we will use a particular setup: we choose to twist only one plane (the x 1 − x 2 plane) by choosing n 12 = −n 21 = 1, while the rest of the components of the twist tensor will be zero. This means that our gauge connections will still be periodic in the x 3 and x 4 directions. As we will see, this choice guarantees that the action has a unique minimum (modulo gauge transformations), and therefore it turns out to be a convenient choice for perturbative studies. This is the reason why the very same choice has been made before to define the Twisted Polyakov Loop running coupling scheme [5], or for other perturbative studies [28]. We will closely follow the notation and steps presented in [26], a reference that the reader interested in more details should consult.
A convenient implementation of twisted boundary conditions consists in using spacetime independent twist matrices. In particular for the periodic directions we set the twist matrices to one We will use latin indexes (i, j, · · · = 1, 2) to run over the directions in the twisted plane, while greek indexes (µ, ν, · · · = 0, . . . , 3) will run over the four space time directions. The consistency relation Eq. (2.2) implies the following condition for the twist matrices Notice that the boundary conditions for the gauge field with this choice of the twist matrices are and A µ = 0 is a valid connection. In fact we will show that it is the only connection compatible with the boundary conditions that does not depend on x, and therefore it is the unique minimum of the action modulo gauge transformations. Eq. (2.7) defines a generalization of the Dirac algebra. It can be shown [25] that there is a unique solution for the matrices Ω i modulo similarity transformations. Introducing the color momentum,p i = 2πñ i N L with n i = 0, . . . , N − 1 it is easy to check that the N 2 matrices where α(p) are arbitrary phases, are linearly independent and obey the relation Moreover all but Γ(p = 0) are traceless, and therefore they can be used as a basis of the Lie algebra of the gauge group. This means that any gauge connection can be expanded as The prime over the sum means that the termp i = 0 is absent in the sum, as required for a SU (N ) gauge group. Notice that the coefficientsÂ µ (x,p) are functions (not matrices) periodic in x. Therefore one can do an usual Fourier expansion and obtain with the usual spatial momentum Finally we define the total momentum as the sum of the color and space momentum P i = p i +p i , P 3,4 = p 3,4 . Noting that any P µ can be uniquely decomposed in the space momentum and color momentum degrees of freedom we can safely write Γ(P ) instead of Γ(p). Our main conclusion is that any gauge connection compatible with our choice of boundary conditions can be written as In particular the only connection that does not depend on x is given byÃ µ (P ) = 0. In general the matrices Γ(P ) are not anti-hermitian, but one can choose the phases α(P ) of equation (2.8) so that this condition is enforced (2.14) In this case, the Fourier coefficientsÃ µ (P ) satisfy the usual relatioñ 15) and the Γ matrices are normalized according to We finally note that a simlar expansion is possible on the lattice, with the only difference that the space momentum will be restricted to the Brillouin zone.

Matter fields
The inclusion of matter fields interacting with gauge fields with twisted boundary conditions is not completely straightforward. To understand why it is better first to consider how to include fermion fields in the fundamental representation. Since the twist matrices tell us how fields change under translations, one naively expects but one can easily see that this choice is not consistent, since the value of the field ψ(x + Lî + Lĵ) depends on the order in which we perform the translations due to the noncommutativity of the twist matrices. This difficulty can be avoided by introducing more fermions, or what usually is called a "smell" degree of freedom [29]. If α, β = 1, . . . , N s are indices that run over the N s smells of fermions, and a, b = 1, . . . , N run over the color degrees of freedom, the boundary conditions of the fermions read This means that a fermion smell becomes a linear combination of the gauge transformed fermion smells under a translation. θ i are in principle arbitrary, but introduced for convenience to remove the zero-momentum modes of the Dirac operator. These phases have to be chosen such that they are not elements of the gauge group (i.e. e ıθ ∈ SU (N )). This choice of boundary conditions for the fermion fields is consistent, but they require the number of smells to be equal to the number of colors. One can easily extend the construction to the case when the ratio N s /N is an integer, but in general one can not have an arbitrary number of fermions in the fundamental representation.
On the other hand fermions in the adjoint representation transform in the same way as the gauge fields, and therefore any number of fermions would be compatible with the twisted boundary conditions.
Regardless of the representation but assuming that the matter fields are compatible with the twisted boundary conditions, O(a) improvement for massless Wilson quarks is automatically satisfied since fields live on a torus, and the boundary conditions do not break chiral symmetry (see [12,13]).

The gradient flow in a twisted box
The gradient flow has recently proved to be an interesting tool to study several aspects of YM theories [14,15,[30][31][32]. By introducing an extra coordinate t, called flow time (not to be confused with Euclidean time x 0 ), gauge fields are smoothed along the flow according to the equation is the field strength tensor. The main reason why renormalization problems are highly simplified with the use of the gradient flow is that correlations functions made of the flow field B µ (x, t) do not need renormalization at positive flow time [15]. In particular the energy density is a renormalized quantity at a scale µ = 1/ √ 8t and can be used (at t > 0) to define a renormalized coupling [14]. Moreover one can use a finite box to define a finite volume renormalization scheme by running the renormalization scale with the linear size of a finite volume box [16,19,23] The constant c parametrizes the ratio between the renormalization scale and the linear size of the box L, and is part of the definition of the renormalization scheme. Being a finite volume renormalization scheme, the boundary conditions are relevant. The idea has been applied in a four dimensional torus with periodic boundary conditions [16], with Schrödinger functional boundary conditions [19], and also with mixed boundary conditions [23]. In this work we will give a definition of a running coupling in a four dimensional torus with twisted boundary conditions, but first we need to study the perturbative behavior of E(t) in a twisted box.
3.1 Perturbative behavior of the gradient flow in a twisted box: continuum

Generalities and gauge fixing
We are interested in the perturbative expression for E(t) , and in order to avoid some difficulties in the definition of propagators, it turns out to be convenient to fix the gauge of the flow field B µ (x, t). This can be achieved by studying the modified flow equation The superscript (α) recalls that covariant derivatives and field strength are made of the modified flow field B Therefore gauge invariant quantities are independent of α. Note that the previously defined gauge transformation Λ(x) belongs to the restricted set of gauge transformations that leave the twist matrices invariant (see equation (2.4)), and the boundary conditions of B (α) µ are also independent of α.

Flow field and energy density to leading order
The particular choice α = 1 simplify the computations, and we will use it for the rest of this section. The modified flow equation reads in this case In perturbation theory one re-scales the gauge potential with the bare coupling A µ → g 0 A µ , and the flow field has an asymptotic expansion in the bare coupling To leading order our flow equation (3.7) is just the heat equation expanding B µ,1 (x, t) in our preferred basis (2.13) one can easily solve (3.9) and obtain Finally our observable of interest also has an expansion in powers of g 0 One can easily obtain Finally using the expression for the gluon propagator one gets

Perturbative behavior of the gradient flow in a twisted box: lattice
When defining the gradient flow in the lattice one has to make several choices. These basically correspond to the particular discretizations of the action whose gradient is used to define the flow, as well as the discretization of the energy density and the choice of action that one simulates (i.e. Wilson/improved actions). First we will analyze the popular case where the Wilson action is simulated, and one uses the same action to define the flow (in this case is called the Wilson flow). The clover definition of the observable has been a typical choice [14] for a discretization of the energy density. Later we will comment on the general case.

Generalities and gauge fixing
On the lattice the gradient flow is substituted by a discretized version. There are several possibilities: one can use the Wilson action where the sum runs over the oriented plaquettes, and define the flow equation by equating the time derivative of the links with the gradient of the Wilson action In this case the gradient flow is usually referred as the Wilson flow. Some explanations of our notation are in order. If f (U µ (x)) is an arbitrary function of the link variable U µ (x), the components of its Lie-algebra valued derivative ∂ a x,µ are defined as . (3.18) In perturbation theory one is interested in a neighborhood of the classical vacuum configuration. In this neighborhood the lattice fields U µ (x) and V µ (x, t) are parametrized as follows: Again it is convenient to study a modified flow equation where the gauge degrees of freedom are damped. We will consider with V Λ µ (x, 0) = U µ (x) and the forward lattice covariant derivativeD Λ µ acting on Lie-algebra valued functions according tô The most natural choice for Λ(x, t) is the same functional used for gauge fixing where∂,∂ * denote the forward/backward finite differences. We again note that the boundary conditions of V Λ µ (x, t) are independent of α, since Λ(x, t) belongs to the restricted class of gauge transformations that leave the twist matrices unchanged.

Flow field and energy density to leading order
Again the choice α = 1 turns out to be convenient and we will stick to it from now on.
The modified flow equation reads The flow field can be expanded in powers of g 0 (equation (3.8)) and to first order in g 0 we have Expanding the flow field in our favorite Lie-algebra basis (equation (2.13)) one can write the solution to the previous equation is the usual lattice momentum. We can choose among different discretizations for the energy density. The most popular one consists in using the clover definition for G µν (x, t) [14]. To leading order we havê where∂ µ = 1 2 (∂ µ +∂ * µ ) is the symmetric finite difference. The energy density computed with the clover definition for the field strength tensor reads Using the definitionsP µ = 1 a sin (aP µ ) , (3.30a) and the lattice gluon propagator, one can easily obtain (3.31)

Some comments on different discretizations 1
In general the lattice computation of the leading order behavior of the energy density involves several choices of discretization: the action that one simulates (labelled (a)), the action whose gradient defines the flow evolution (labelled (f)), and finally the discretization used to compute the observable (labelled (O)). To leading order, these three choices can be expressed as choice of "actions" The matrices K (a) and K (f ) may (and should) contain a gauge fixing part, but not the one corresponding to the observable K (O) . In this way final results will be independent of the choices of gauge. The inverse of the K (3.34) Using this notation it is trivial to obtain the form of the flow field to leading order and noting that the reality of the action requires that H + (t, P ) = H(t, −P ), we can write the expression of the energy density to leading order as This formula allows an easy evaluation of the energy density, to leading order in perturbation theory, for any choice of discretizations. One general point that one can make is that if one uses the Wilson flow the matrix H(t, P ) can be chosen to be proportional to the identity (by an appropriate gauge choice), and therefore commutes with any other matrix. Moreover if the action that one simulates is the same as the one that we use to compute the observable, the product of matrices DK O together with the trace simply result in a factor 3, and therefore one obtains This means that without changing the flow, improving the action and the observable leads to exactly the same cutoff effects than if one does not improve anything (to leading order).

Tests
In order to check the previous computations one can perform several consistency checks. First it is obvious that the continuum result (equation (3.15)) is recovered from the lattice one (equation (3.31)) if one takes the limit a/L → 0. In the infinite volume limit boundary conditions are irrelevant, and therefore for L → ∞ one should recover the result of [14] that reads This result is reproduced from our expression equation (3.15) by simply noting that and therefore 1 (3.41) Finally recalling that there are N 2 −1 terms in the sum overp i (the termp i = 0 is explicitly excluded) one obtains To check the lattice computations we have performed some dedicated pure gauge lattice simulations. We use the plaquette action of an SU (2) gauge theory in two different volumes L/a = 4 4 and L/a = 6 4 . We collect 10, 000 measurements of E cl (t) for different values of t and β = 2/g 2 0 = 40, 80, 120, 200, 400, 560, 800, 960, 1120, 1280. In these large-β simulations the measured E cl (t) should reproduce the perturbative expression (equation (3.31)). Being more precise, we will study numerically the quantity We expect that R(g 0 , t) = O(g 2 0 ), and therefore by fitting the data from the simulations to a linear behavior R(g 0 , t) = m(t)g 2 0 + n(t)

Running coupling definition
The computations of the previous section guarantee that . On the other hand the properties of the gradient flow ensures that t 2 E(t) is a renormalized observable defined at a scale µ = 1/ √ 8t. This suggests that one can use t 2 E(t) for a nonperturbative coupling definition. Moreover if one keeps the product of the renormalization scale and the linear size of the box fixed (i.e. µL = 1/c = constant) the coupling will depend on no scale other than the linear size of the box, and therefore will be ideal for finite size scaling.
In full glory our coupling definition reads This coupling has a perturbative expansion and a universal two loop β-function where the universal coefficients, for the case of an SU (N ) YM theory are given by We point out that the same coupling definition is valid if one includes any number of fermions in any representation, as soon as they are allowed by the twisted boundary conditions (more details in section 2.2).

Cutoff effects in the twisted running coupling
The comparison of the lattice and the continuum computations of E(t) can give us an idea of the size of cutoff effects (to leading order in g 2 0 ) of the twisted gradient flow coupling. We are going to study in detail the case of lattice simulations using the Wilson action, the Wilson flow, and the clover definition for the observable. If we definê

Improved coupling definition
If one is computing t 2 E(t) non-perturbatively via lattice simulations, and one is using the Wilson action, the Wilson flow and the clover observable for the evaluation of the energy density observable, one can alternatively define the coupling via This last coupling definition has the same properties, but one expects an improved scaling towards the continuum limit, since the leading order cutoff effects ∝ g 2 0 have been removed thanks to the lattice factorN T (c, a/L).
In a similar way, any choice of discretizations that define a coupling can be normalized with a factor computed on the lattice (cf. section 3.2.3), leading to an improved scaling towards the continuum.

SU (2) running coupling
In this section we will compute the running coupling in SU (2) pure gauge theory to test if the twisted coupling definition is applicable for step scaling studies.
We will first recall the general strategy, introduced in [3], of the recursive procedure involved in the computation of the running coupling. The process starts with a nonperturbative coupling definition that depends on no scale other than the linear size of a finite-volume box. Of course in our case this role will be played by the twisted gradient flow coupling g 2 TGF (L). The step scaling function tells us how much the coupling changes when the renormalization scale is changed by a factor s σ s (g 2 (L)) = g 2 (sL) .
(5.1) Therefore it is a discrete version of the β−function. The value s = 1/2 is a typical choice, and is the one we will use for now on, although the basic idea is the same for any other value. If one knows how much is the coupling at a renormalization scale that corresponds to a large volume (lets call it L max ), and one knows the step scaling function, then one can obtain the value of the coupling at scales L max /2 k for k = 1, 2, . . . . Eventually, one will reach a very small box size (very large energy scale), where asymptotic freedom guarantees that one can safely make contact with perturbation theory.
To compute the step scaling function numerically one starts by measuring the coupling on some lattice of size L/a. The step scaling function is easily obtained by simply measuring the coupling on a lattice half as big (L/(2a)) while keeping the rest of the bare parameters constant. The step scaling function computed in this naive way will carry an implicit dependence of the lattice spacing (the cutoff), and therefore defines a lattice step scaling function Σ(g 2 (L), a/L) .

(5.2)
In order to obtain the continuum step scaling function σ(g 2 (L)), one simulates several pairs of lattices and takes a continuum limit: Σ(u, a/L) .

Simulation details
We will simulate SU (2) YM theory using the Wilson action where the sum runs over all oriented plaquettes. We simulate lattices of size L/a = 20, 24, 30, 36, and in order to compute the step scaling function also lattices of half this size (L/a = 10, 12, 15, 18). The range of values of β (between 2.75 and 12.0) translate to renormalized couplings g 2 TGF (L) between 7.5 and 0.6 (for c = 0.3), enough to cover both the non-perturbative and perturbative regions of the theory. Appendix A collects the values of the g 2 TGF (L) of our simulations. We will use a combination of heatbath [33][34][35] and overrelaxation [36] as suggested in [37]. In particular we choose to do one heatbath sweep followed by L/a overrelaxation sweeps. Since measuring the coupling (i.e. integrating the flow equations) is numerically more expensive than the Monte Carlo updates, we repeat this process 50 times between measurements.
In total we collect 2048 measurements of the coupling for each lattice size, each value of β, and several values of c ∈ [0.3, 0.5]. These measurements are collected in N r parallel runs (replica) of length N MC each so that N r × N MC = 2048. We check that there are no autocorrelation between measurements (i.e. τ int = 0.5 within errors), even for our larger lattices and larger values of c. We conclude that we can safely consider the measurements independent.
The Wilson flow equations are integrated using the adaptive step size integrator described in appendix D of [19]. With this scheme we make sure that the integration error in each step is not larger than 10 −6 .

Data analysis
For each L/a we have computed the value of the twisted gradient flow coupling at different values of β (we call it g 2 TGF (β; L/a)). These data are fitted to a Padé-like ansatz This fit imposes the one-loop constraint to the data (i.e. g 2 TGF (β; L/a) → 4/β at large β), and has a total of 2M free fit parameters.
Alternatively, and to estimate the dependence of our results on the choice of functional form used to fit the data, we use a different Padé inspired functional form that also ensures the correct one-loop behavior at large β.
We obtain good fits (χ 2 /ndof ∼ 0.6 − 1.9) with M = 2 when using the functional form of Eq. 5.5 to fit the lattice data (i.e. 4 fitting parameters). When using the functional form of Eq.5.6 we need M = 4 to accurately describe the data on the small lattices (L/a = 10, 12) and M = 6 for the larger ones (L/a = 15, 18, 20, 24, 30, 36). It is important to stress that the data are statistically uncorrelated, since they correspond to different simulations.
In the figures 3 we show a couple of these fits. Our worst fit corresponds to the L/a = 24 lattice and the Pade fit gives a χ 2 /ndof = 1.69, while the Taylor fit results in a fit quality of χ 2 /ndof = 1.9. We see how in this case the two different functional forms interpolate differently between the data, giving us confidence that if one estimates the error of the interpolation using both functional forms, one will be on the safe side 2 .
We use resampling methods to propagate errors by using 4000 bootstrap samples. All fitting parameters derived from our original data are computed for each bootstrap sample. Interpolation points are computed for each bootstrap sample and each functional form. The final error of the interpolated point is computed using both functional forms and all bootstrap samples, and therefore takes into account not only the statistical uncertainty, but also the systematic effect due to the dependence of the interpolating functional form.

Step scaling function
We will first show the continuum extrapolations of the step scaling function Σ(u, a/L) at some representative values of u = 7.5, 3.75, 1.5. Figure 4 shows that these extrapolations are mild. We have used the value c = 0.3 that gives a precision in the data for the renormalized coupling between 0.15% and 0.25%. We stress that this systematic effects is taken into account in our analysis by using both functional forms to estimate the error of the interpolations (see the text for more details). (b): Fits to the data of the L/a = 36 lattice. As we can see, in this case both interpolating functions agree within errors, although the polynomial fits tends to have larger errors.
One of the advantages of the use of the twisted boundary conditions is the absence of O(a) cutoff effects, that are present for example in the Schrödinger functional due to boundary effects. Here the invariance under translations guarantee that the continuum limit can be safely taken by a linear extrapolation in (a/L) 2 .

Running coupling
As a final application, we will compute the running coupling. We will fix the scheme by setting c = 0.3. We start our recursion in a volume L max defined by the condition The lattice step scaling function and its continuum limit is computed as described in the previous sections. As figure 4 shows, the extrapolations towards the continuum are rather flat. The continuum limit values are used to further compute the values of the step scaling function at larger renormalization scales (smaller volumes), up to L min = L max /2 26 , where g 2 TGF (L min )| c=0.3 = 0.5324(84).
Since the same functional form (fitting parameters) are used recursively to compute the values of the coupling at different scales, one has to propagate errors taking into account the correlations correctly. This is done in the spirit of the resampling methods in the most naive way: one uses as input for the coupling at a scale L all the bootstrap samples of the coupling from the scale 2L. We recall here that these bootstrap samples carry the information not only of the statistical uncertainties, but also of the dependence of our results on the functional form chosen to fit the data. Our results have carefully taken into account the two sources of systematic uncertainty: the continuum extrapolation and the choice of fitting function for our lattice data.   Figure 5 shows the running of the coupling from the low energies to the high energies, over a factor 2 26 change in scale, while table 2 contains the numerical values of the coupling at different renormalization scales. The fact that the absolute error in the renormalized coupling tends to be constant a large energies (small volumes), is a consequence of the error propagation, that dominates for large energies the error budget.
As a further consistency test, we have repeated the full running of the coupling using as scale factor to define the step scaling function s = 2 (i.e. we run from high to low energies), obtaining consistent results.
The Λ parameter can be extracted, in units of L max via using that µ = 1/cL. The previous formula is exact, but the last exponential is essentially unknown analytically. Nevertheless if one uses a value of g 2 TGF (L) where the difference between the two loop and the non-perturbative results are negligible, the effect of the last  exponential is also negligible. Of course this is more certain the smaller the coupling, but since the relative error of the coupling grows as the coupling decreases, this would translate in a larger error for the Λ parameter. Below we quote a couple of values as example.
We want to end this section with a small comment on the use of different values of c. The main point has already been raised in [19]: the larger the value of c, the larger the (relative) statistical error of the coupling, but the scaling towards the continuum seems better. This general behavior is consistent with the leading order in perturbation theory as we have seen. We will simply say that the relative error in the raw data increases with c, and roughly one can say that for c = 0.4 the relative error is two times larger than for c = 0.3, while for c = 0.5 the error is three times larger. This statement seem to hold true independently of the volume (i.e. of the value of g 2 TGF ).

Conclusions
In this work we have studied the YM gradient flow for fields obeying twisted boundary conditions on a 4D torus. The energy density of the flow field at positive flow time is used to define a non-perturbative coupling at a scale given by the flow time. Moreover one can make the flow time scale with the finite size of the box. In this way one obtains a coupling that depends only on a scale given by the size of the torus, and the usual techniques of step scaling can be applied. The perturbative behavior of the energy density of the flow field is computed to leading order both in the lattice and in the continuum. This allows us to estimate the size of cutoff effects of the coupling to leading order in perturbation theory. Results are consistent with other definitions of the running coupling using the gradient flow. Cutoff effects are mild.
In order to state the validity of this coupling definition beyond perturbation theory we have performed a numerical study in pure gauge SU (2). We have computed the running of the coupling from the very perturbative regime g 2 ∼ 0.5 to the non perturbative one g 2 ∼ 7.5, over a change by a factor 2 26 in scale. We have used lattices of sizes L/a = 10−36. The statistical precision that one can achieve with this coupling definition is very high: 2048 independent measurements are enough to achieve a sub-percent precision (0.15 − 0.25%). This precision is roughly independent of the physical volume, the value of the coupling or the lattice size. We have also shown that the continuum extrapolations of the step scaling function are rather flat. Moreover we have shown that this techniques can be used to reliably extract the Λ parameter.
The same coupling definition can be used if fermions are coupled to the gauge field. But the consistency of the twisted boundary conditions imposes a constraint on the number of fermions in the fundamental representation that we can use. In particular, what arguably is the most interesting application of these techiniques, namely the coupling constant of QCD (gauge group SU (3) with N f = 4 fermions in the fundamental representation) can not be attacked with this method. Nevertheless there are many applications of this work. On one hand N f = 3 flavour of quarks in the fundamental representation is already interesting for the strong interations. But it is probably in the context of conformal theories and physics beyond the standard model, with models like SU (3) with N f = 12 fermions in the fundamental representation, or with adjoint fermions, where the nice properties of this coupling definition (automatic O(a) improvement, analyticity and high statistical precision) can have a higher impact. In particular the setup with twisted boundary conditions have already been used for step scaling studies in the TPL scheme [38,39]. In this particular case the use of the gradient flow observable will lead to more precise results, as we have already seen in some preliminary results in the last lattice conference [40]. One can also consider the application of the same ideas to other choices of twisted boundary conditions. This is specially interesting in the context of reduced models, as we have also seen recently [41].   Table 3: Raw values of g 2 TGF (β; L/a).