The gradient flow coupling at high-energy and the scale of SU(3) Yang-Mills theory

Using finite size scaling techniques and a renormalization scheme based on the Gradient Flow, we determine non-perturbatively the $\beta$-function of the $SU(3)$ Yang-Mills theory for a range of renormalized couplings $\bar g^2\sim 1-12$. We perform a detailed study of the matching with the asymptotic NNLO perturbative behavior at high-energy, with our non-perturbative data showing a significant deviation from the perturbative prediction down to $\bar{g}^2\sim1$. We conclude that schemes based on the Gradient Flow are not competitive to match with the asymptotic perturbative behavior, even when the NNLO expansion of the $\beta$-function is known. On the other hand, we show that matching non-perturbatively the Gradient Flow to the Schr\"odinger Functional scheme allows us to make safe contact with perturbation theory with full control on truncation errors. This strategy allows us to obtain a precise determination of the $\Lambda$-parameter of the $SU(3)$ Yang-Mills theory in units of a reference hadronic scale ($\sqrt{8t_0}\,\Lambda_{\overline{\rm MS}} = 0.6227(98)$), showing that a precision on the QCD coupling below 0.5\% per-cent can be achieved using these techniques.


Introduction
Yang-Mills (YM) theories play a central role in our understanding of natural laws. They lie at the heart of the unification between the electromagnetic and weak interactions and at the foundations of Quantum Chromo Dynamics (QCD). The (pure) SU (3) YM theory shares many of the interesting features of QCD. It is a strongly coupled non-abelian gauge theory where the fundamental degrees of freedom, the gluons, are not part of the spectrum of the theory. The theory has an intrinsic energy scale given by the Λ-parameter, and at energy scales much larger than Λ it is well-approximated by perturbation theory (i.e. the theory is asymptotically free [1,2]). Connecting this high energy perturbative regime of the theory with its low energy spectrum is a difficult multi-scale problem, very similar to the one faced when one wants to determine the strong coupling or quark masses in the Standard Model (SM). On the other hand, the pure gauge theory is much more tractable from a computational point of view when lattice field theory methods are employed. The currently known simulation algorithms are far more efficient in the case of the pure gauge theory than in the case of QCD. In summary, the SU (3) YM theory is very interesting in its own, and a perfect laboratory for testing new lattice techniques and ideas before applying them to QCD.
In this work we are interested in the determination of the intrinsic energy scale of the SU (3) YM theory, i.e., its Λ-parameter. In the case of QCD this is equivalent to determining the value of the strong coupling. We make use of lattice field theory methods [3,4], which allow us to determine non-perturbatively the running of the gauge coupling. In particular, thanks to the techniques of finite size scaling [5,6], this running can be computed over a wide range of renormalization scales. These two ingredients of our analysis are what allows us to bridge the 3 orders of magnitude that separate the typical energy scales of the hadronic, strongly coupled regime of the theory, and the high energy perturbative regime, without making any assumption on the validity of perturbation theory (see the discussion in [7,8]).
A similar study using the Schrödinger Fucntional (SF) coupling was one of the earliest applications of step scaling techniques in this context [9]. More than 20 years later, we are now capable of a much more precise computation. This not only because of the increase in computational power over the years, but also thanks to the development of new theoretical tools. In particular, the Gradient Flow (GF) [10,11] allowed us to introduce new coupling definitions [11,12,13,14] which are very compelling for step scaling studies. These couplings are constructed from observables with very small variance and very special properties under renormalization (see [15] for a comparison among different coupling definitions). As a result, GF couplings permit to achieve an excellent statistical precision, especially in the low energy regime of the theory, where the traditional SF coupling struggles to produce precise results.
GF-based couplings, however, have their issues, too. One of the main issues are the relatively large cutoff effects that have been observed in several studies (see [15] for a discussion). Despite a solid theoretical understanding on the anatomy of these discretization effects [16], these remain the main source of concern in most applications. Thus, using GF couplings one can achieve very high statistical precision, but accurate continuum results require the simulation of large lattices. There are also other challenges if one wants to use these couplings for a precise determination of the strong coupling. First, the known perturbative coefficients of the β-function show bad convergence [17,18]. This implies that even simple estimates of the theoretical perturbative uncertainties in the extraction of α s based on GF couplings, even at energies as high as the electroweak scale, are about half a percent. Second, the relative statistical precision on the GF couplings is typically δα GF /α GF ∝ const., while the traditional SF coupling has δα SF /α SF ∝ α SF . This means that eventually the numerical cost of the GF couplings will be larger than for the SF couplings when α 1. These issues explain the strategy followed by the ALPHA collaboration in the determination of α s [8,19,20,7], where the SF coupling is still the coupling of choice at high energies (see also refs. [21,22] for recent reviews).
In this work we will explore in detail the high energy sector of the SU (3) gauge theory using the Gradient Flow. One of our aims is to assess whether precise results for α s can be achieved using GF couplings. As mentioned earlier we will have to face the large truncation errors in these schemes due to the bad behaviour of their perturbative series. We will show that, in fact, due to the large truncation effects, it is very challenging to obtain precise results for the Λ-parameter using only these schemes. We will also pay special attention to the linear O(a) effects that typically arise in finite volume renormalization schemes that break translational invariance. They are the main source of systematic effects in computations based on the SF couplings [7], and we will see that GF schemes allow us to substantially reduce these effects.
The paper is organized as follows. Section 2 introduces the GF and SF schemes in our preferred finite volume setup. Section 3 introduces our numerical setup and our determination of the running of the GF coupling at high energies. In this section we shall discuss in detail the matching with the asymptotic perturbative behaviour. In section 4 we complete our determination of the running coupling in the non-perturbative domain and match our finite volume renormalization schemes with the infinite volume reference scales t 0 , and r 0 . Section 5 presents our final results and offers a comparison with the data available in the literature, while section 6 presents the conclusions of our work. A few appendices are moreover included to address some more technical details. In appendix A we provide evidence that our results for the GF coupling, that span a factor 3 in lattice spacing, are suitable to produce accurate continuum limit extrapolations. Appendix B details our procedure to estimate the boundary O(a) effects. Appendix C presents a more traditional analysis of our high energy data, and the related extraction of the Λ-parameter. Appendix D discusses the perturbative improvement we applied to the SF coupling. Finally, appendix E contains all raw data measurements of our simulations.

Running couplings and the Λ-parameter
The SU (3) Yang-Mills theory is a fairly simple theory. Formally, it is defined in terms of a gauge potential A µ which lives in the Lie algebra su(3) of SU (3), and by the action: where F µν is the field strength tensor, and g is a dimensionless free parameter: the gauge coupling. The renormalization of the theory requires us to introduce a renormalized coupling,ḡ ≡ḡ(µ), through some suitable renormalization condition; different conditions define what we refer to as different renormalization schemes for the coupling. Renormalized couplings depend explicitly on the energy scale, µ, at which their defining conditions are imposed. The dependence on this scale is encoded in their β-functions, β(ḡ), which are defined by the renormalization group (RG) equations: 3) The β-functions have an asymptotic perturbative expansion: where the first two coefficients, b 0 = 11 (4π) 2 and b 1 = 102 (4π) 4 , (2.5) are universal, i.e., they are independent of the specific renormalization scheme chosen for the coupling.
The scheme dependence only enters through the higher-order coefficients b k , with k > 2. 1 This implies that although at "low" energy different coupling definitions can behave very differently as a function of µ, at high-energy these differences must eventually disappear, as all definitions share the property of asymptotic freedom:ḡ(µ) µ→∞ → 0. The first-order RG equation (2.3) has as an implicit solution given by, 2 (µ) exp{−I g (ḡ(µ), 0)}, I g (g 2 , g 1 ) = where Λ is a constant of mass dimension one. Its value depends on the exact renormalization scheme chosen for the coupling. From eq. (2.6) it is clear that, given the knowledge of the β-function, the Λ-parameter is all that is needed to infer the value of the couplingḡ(µ) at any renormalization scale µ. In particular, the coupling must be a function of µ/Λ, which means that Λ defines what is to be considered "low" or "high" energy. The Λ-parameter is a compelling quantity to determine. First of all, its definition does not rely on perturbation theory: it is non-perturbatively defined once the corresponding coupling and β-function are. Second, it is a renormalization group invariant. As such it does not depend on any renormalization scale, i.e. dΛ/dµ = 0. In addition, even though its value depends on the scheme, this dependence can be computed analytically. Given two renormalized couplingsḡ X andḡ Y and the one-loop perturbative relationḡ with c 1 a pure number, one can easily show using eqs. (2.4)-(2.6) that the corresponding Λ-parameters are exactly related by: These properties make any Λ-parameter a natural reference scale for both the low and high energy regimes of the theory. In particular, any dimensionfull renormalization group invariant quantity, like for instance any "hadronic" quantity, must be proportional to Λ (or some power of it). 2 The proportionality constants that relate these quantities to Λ are fundamental numbers that characterize the pure Yang-Mills theory, as they are given solely by the dynamics. Lattice field theory is the only known framework that allows us to extract low-energy physics from first principles. In particular, the value of Λ in units of a typical hadronic scale µ had , for whichḡ(µ had ) = g had , can in principle be obtained by employing these techniques and eq. (2.6). The determination requires of course the knowledge of the corresponding β-function for all energies larger than µ had . The basic strategy for this computation is to divide this energy range into two parts. While the nonperturbative β-function is computed from µ had up to some energy µ PT , for whichḡ PT =ḡ(µ PT ) ḡ had , perturbation theory is used for energies larger than µ PT . For large enough µ PT , the relevant integral entering the definition of Λ/µ had can be approximated as (cf. eq. (2.6)): 9) where I N g is analogously defined as I g by replacing the β-function with its perturbative expression up to some order N , i.e., The ratio of interest, Λ/µ ref , is thus approximated by: (2.11) There are two important points we must stress here. First, in the extraction of the Λ-parameter the truncation errors are formally of O(ḡ 2N −2 PT ) (cf. eq. (2.11)). Since the running of the coupling at high energies is logarithmic in µ/Λ, reducing the size of these truncation errors of a given factor, requires a significantly larger change in the energy scale. Second, one must always remember that eq. (2.11) is only an asymptotic statement. In principle, non-perturbative corrections (e.g. power corrections) are also present when approximating the integral as in eq. (2.9). As we shall clearly see in the following, these issues imply that in order to accurately estimate the systematic uncertainties coming from the use of perturbation theory at high-energy, as well as to keep these uncertainties small, we must study the applicability of perturbation theory over a wide range of energies, reaching up to very large scales. Therefore, a precise determination of the Λ-parameter in terms of low-energy scales is no easy task.

Finite size scaling
It is certainly not obvious how in the strategy presented above one can reach large energy scales via lattice field theory simulations. Indeed, the ultra-violet cutoff of the lattice theory set by 1/a, where a is the lattice spacing, has to be much larger than the largest energy scale one wants to reach; large discretization errors will otherwise affect the results. At the same time, in order to have finite-volume effects both in the coupling and in the hadronic quantities well under control, the infra-red cutoff set by the finite extent of the lattice, L, has to be at least a few femto-meters long. Given the fact that current computational resources allow us to simulate lattices with lattice sizes L/a ∼ O(10 2 ), this significantly limits the range of energy scales that can actually be covered in a single lattice simulation.
Finite-size scaling techniques overcome these problems by integrating the RG equations nonperturbatively in a finite-volume renormalization scheme [5]. Within this strategy the coupling is defined through a finite volume observable, and the renormalization scale is identified with the infrared cutoff, i.e., µ = 1/L. Said it differently, we define a renormalized coupling through a finite volume effect. In this way, high renormalization scales can be reached by splitting the computation into several lattices of smaller and smaller physical size. A quantity of primary role in finite-size scaling studies is the step scaling function: It is a discrete version of the β-function as it measures the change in the coupling when the renormalization scale is varied by a finite factor s. Note that once the step scaling function is known, the β-function can also be determined, by simply noticing that: (2.13)

Renormalization schemes
Analytic calculations in gauge theories are simplified by considering for the coupling the so called modified minimal subtraction (MS) scheme of dimensional regularization. Although this is a convenient choice for perturbative calculations, the MS scheme is only defined within perturbation theory and therefore only applicable at high energies. In this work we are going to employ different renormalization schemes, which are non-perturbatively defined in a finite Euclidean space-time volume. In order to apply the idea of finite-size scaling, the renormalization scale of these couplings must be linked to the physical size of the system, i.e. µ ∝ 1/L. At high energies, perturbation theory can then be used to relate any of these schemes to the more conventional MS scheme (cf. eq. (2.7)). Note that thanks to eq. (2.8), Λ MS can be implicitly defined non-perturbatively through any non-perturbative scheme, even though the MS scheme itself is intrinsically perturbative.

Boundary conditions and coupling definitions
When studying the pure Yang-Mills theory in a finite space-time volume the choice of boundary conditions for the fields matters. In this work we consider the SU (3) Yang-Mills theory with Schrödinger functional (SF) boundary conditions [6]. In this set-up, the gauge field is periodic in the three spatial directions with period L, while its spatial components satisfy Dirichlet boundary conditions at Euclidean times x 0 = 0, L, i.e., where C k (x), C k (x) ∈ su(3) are given external fields. A first compelling feature of this type of boundary conditions is that, for a proper choice of fields C k , C k , the action has a unique global minimum (up to gauge transformations). This avoids several complications when doing perturbation theory in a finite volume with respect to the general case [29]. Another advantage of using these boundary conditions is that the system can be probed by considering different boundary fields C k , C k . In the following we focus on a particularly convenient family of fields, given by the Abelian, and spatially constant fields of the form [9]: where η, ν are some (dimensionless) real parameters. Derivatives of the effective action of the pure Yang-Mills theory with SF boundary conditions with respect to the parameter η are renormalized quantities [6]. They can hence be used to define renormalized couplings. In particular, we can define a family of SF couplings as [6,9,30,31,8,7]: where different values for the parameter ν define different renormalization schemes. In fact, this family of couplings can be obtained from a linear combination of two distinct observables both defined for ν = 0, specifically, where, in terms of expectation values: The perturbative relation between the SF couplings defined above and the coupling in the MS scheme, is known to two-loop order, and it is given by [9,32,33,34] (α X ≡ḡ 2 X /4π): where r > 0 and a ν 1 (r) = a 1 (r) + 4πv 1 ν , a ν 2 (r) − (a ν 1 (r)) 2 = a 2 (r) − (a 1 (r)) 2 + (4π) 2 v 2 ν , and v 1 = 0.0694603(1), v 2 = −0.001364 (14).
The knowledge of the two-loop relation (2.20) allows us to determine from the 3-loop β-function in the MS scheme [35,36], the 3-loop β-function of the SF couplings, which means to determine the first non-universal coefficient (cf. eq. (2.4)): Using eq. (2.8) in conjunction with eq. (2.20), we can also obtain the ratio of the corresponding Λparameters: The result of eq. (2.24) shows that for values of |ν| = O(1), the b SF,ν 2 coefficients are "naturally" small compared to the lowest-order results of eq. (2.5). Within perturbation theory, one is hence keen to expect that for these SF ν schemes truncation errors in the β-function are small at small values of the coupling α ≡ḡ 2 /(4π) 1. Other compelling coupling definitions are possible within the SF framework. Of particular interest for our study are couplings defined through the Yang-Mills gradient flow (GF). The Gradient Flow evolves the gauge field according to a diffusion-like equation: where D µ = ∂ µ + [B µ , ·] denotes the gauge covariant derivative of the field B µ , and is the corresponding field strength tensor. The flow time t has units of length squared, and B µ (t, x) can be seen as a smoothed version of the original gauge field A µ (x) over a length scale ∼ √ 8t. The remarkable property of the flow fields is that gauge invariant operators made out of these fields are renormalized quantities for positive flow times, t > 0 [37]. This suggests that, for instance, the dimensionless quantity 3 can be used to define renormalized couplings at a scale given by the (inverse) flow time, e.g. µ = 1/ √ 8t. Clearly, this coupling definition does not rely on having specific SF boundary conditions, and for what matters in having a finite-volume either. Hence, one can take for eq. (2.15) the convenient choice: C k = C k = 0. Due to the explicit breaking of rotational symmetry by the boundary conditions, one can in fact obtain two independent coupling definitions by considering either the magnetic or the electric components of the energy density, eq. (2.28) [13], i.e., x)}, (2.30) 3 We use the convention that tr{T a T b } = − 1 2 δ ab , where T a , a = 1, . . . , 8, are generators of su (3).
where N m , N e are some constants that guarantee the correct normalization of the coupling [13]. In order to properly define a coupling suitable for finite-size scaling, we must relate the renormalization scale at which the coupling is defined with the size of the finite-volume, we thus set: µ = 1/ √ 8t = 1/(cL), where the constant c is part of the scheme definition. In this work we will exclusively take c = 0.3; the merits of this choice have been discussed in ref. [13]. We note while passing that, at fixed flow time, the limit c → 0 corresponds to the analogous GF coupling definition in infinite space-time volume [11]. In addition, we choose to measure the flow energy densities for x 0 = L/2 in order to maximize the distance of the observables from the space-time boundaries. When looking at the couplings on the lattice, this minimizes the O(a) contaminations coming from the lattice action (cf. Sect. 3.1).
Perturbative computations for flow quantities are usually involved. Already the 1-loop relation of the GF coupling with the MS coupling in infinite volume [11] is a challenging computation. On a finite volume the computation is even more involved (see [38]). The two-loop relation requires substantial effort [17,18], and for our particular choice of boundary conditions the result relies on novel methods [39,40,41,18] within the framework of Numerical Stochastic Perturbation Theory (NSPT). The results are [18]: where the coefficients k m/e 1,2 are collected in table 1. As for the case of the SF couplings, the 2-loop relations (2.31) allow us to infer the 3-loop coefficients of the β-functions of the GF couplings using the known results in the MS scheme, this yields: It is interesting to compare these results with those of the analogous GF coupling definition in infinite space-time volume (c = 0) [11,17]. For this scheme, the electric and magnetic definitions coincide, and we have [17]: Given the above results, it is clear that while the magnetic results for b GF,m 2 are significantly larger than for the infinite volume case, those for the electric components are similar. 4 In all these cases, however, the 3-loop coefficient is "unnaturally" large and of opposite sign if compared with the lower-order ones (cf. eq. (2.5)). It is also significantly larger that the SF ν schemes previously discussed (cf. eq. (2.24)). Already from a purely perturbative point of view, one is thus worried that higher-order corrections to the β-function may be large, even if the coupling is relatively small. These concerns will be in fact confirmed by our non-perturbative investigation.
To conclude, for the following it is also useful to work out the perturbative two-loop relation between the SF and the GF couplings. Combining eq. (2.20) with (2.31), we have:  3 The gradient flow coupling at high-energy

Lattice set-up
We regularize the SU (3) Yang-Mills theory in terms of the link variables U µ (x) ∈ SU (3), on a lattice of size L/a in all four space-time dimensions, where L is the physical size of the lattice and a is its spacing. For the lattice action we take the standard Wilson (plaquette) gauge action: where the sum is over all the plaquettes of the lattice, and U p denotes the product of the gauge links around the plaquette p. β = 6/g 2 0 , with g 0 the bare gauge coupling. Due to our choice of SF boundary conditions we included in eq. (3.1) the weight factor w(p): The coefficient c t can in principle be tuned to cancel the O(a) discretization errors stemming from the boundary of the lattice [6]. Unfortunately, however, only perturbative estimates are currently available [9,32,33]. In the following, we employ the two-loop result [33]: and estimate through dedicated simulations the effect of having the coefficient truncated at this order (cf. Appendix B). On the lattice, the SF boundary conditions are imposed by setting [6]: where C k (x), C k (x) are either equal to (2.15) with η = ν = 0, or to C k = C k = 0, depending on whether we are interested in measuring the SF or the GF couplings. Starting from these definitions a lattice regularization of the SF couplings naturally follows [6,42,9]. We refer the reader to the original references for the details. For the case of the GF couplings, instead, there is quite more freedom in their lattice definition. First of all, we need to specify a discretization for the flow equations (2.26). A popular choice is the Wilson flow (no summation over µ) [11]: where V µ is the lattice flow field and ∂ x,µ S W [V ] is the force deriving from the Wilson action, eq. (3.1). The Wilson flow describes the continuum flow equations up to O(a 2 ) errors. 5 This can be improved to O(a 4 ) by considering the Zeuthen flow (again no summation over µ) [16]: where ∂ x,µ S LW [V ] is now the force deriving from the Symanzik tree-level O(a 2 ) improved (Lüscher-Weisz) gauge action, S LW [44]. 6 Further details can be found in ref. [16]. Here we just want to comment that the term proportional to ∆ µ = ∇ * µ ∇ µ is included for all links, except for those temporal links touching the SF boundaries at x 0 = 0, L. In this case we set ∆ 0 = 0.
In order to define the GF couplings on the lattice we also need to specify a discretization for the E e/m -fields entering the definitions, eqs. (2.29). We consider the following two options (cf. ref. [16]). In the case where the lattice flow is given by Wilson flow (3.5), we choose to discretize E e/m in terms of the clover definition of G µν (see also ref. [11]). On the other hand, in the case where the Zeuthen flow (3.6) is used, we take the O(a 2 ) improved combination: E e/m = 4 3 E pl e/m − 1 3 E cl e/m , where E pl and E cl , are the energy densities discretized in terms of the plaquette action density and the clover definition of the flow field strength tensor, respectively [16]. Based on an analysis in terms of Symanzik effective theory the Zeuthen flow/improved observable combination is preferable, as this choice does not introduce O(a 2 ) effects when integrating the flow equations or when evaluating the operators at positive flow times. In this case, indeed, O(a 2 ) discretization effects come only from the action, eq. (3.1) (which also introduces O(a) effects), and from the incomplete knowledge of a flow improvement coefficient, c b (g 0 ), at t = 0 (see [16] for a complete discussion). The numerical experience gained so far seems to indicate that this results in a better O(a 2 )-scaling for the Zeuthen/improved observable combination than the Wilson flow/clover one (see e.g. ref. [19]).
Before giving the final expression for our lattice definition of the GF coupling we must address one last important point. It is well-known that numerical simulations of the SF at lattice spacings a 0.05 fm and with L 0.5 fm, tend to show large autocorrelation times due to the infamous problem of topology freezing [45,46,43]. As suggested in [46], this problem can be circumvent by defining the coupling within the sector of topologically trivial gauge fields. The continuum definitions (2.29)-(2.30) are replaced in this case by: where δ Q is a Dirac δ-function that enforces the topological charge Q of the gauge fields integrated over in the functional integral to be zero. We note that this modification actually defines different renormalization schemes than eqs. (2.29). On the other hand, these two set of couplings are indistinguishable from a perturbative point of view and thus share the very same perturbative results given in Sect. 2.3. On the lattice, we can finally define the new GF couplings through the expression: To define the topological charge on the lattice we use the clover discretization of the flow strength tensor [11]: measured at flow time √ 8t = cL. For the discretization of the flow equations used for Q we will always use the one employed for the discretization of E e/m entering the coupling definition. In addition, since on the lattice Q is not integer-valued, we replace the Dirac δ-function with:  Table 2: Lattice norms, eq. (3.8). In all cases we use c = 0.3. The labels W/Z refers to the Wilson/Zeuthen discretization and the labels m/e to the magnetic/electric GF coupling. See text for more details.
We conclude by noticing that, as proposed in [13], the normalization factorsN e,m (c, a/L) are better computed in lattice rather than continuum perturbation theory, using the very same lattice discretization employed in the simulations (which includes the definition of the lattice action, flow, and observable). 7 This guarantees that the exact relation:ḡ 2 GF, e/m = g 2 0 + O(g 4 0 ), holds. All discretization effects are hence removed at tree-level in perturbation theory. For completeness, we collect in table 2 the relevant values of the coupling norms used in this study.

Lattice step-scaling function
As discussed in Sect. 2.2, the strategy to determine the non-perturbative running of a renormalized coupling using finite-size scaling techniques relies on the computation of the step-scaling function where for a finite-volume renormalization scheme, µ ∝ 1/L. The corresponding β-function can then be determined from σ s (u) using relation (2.13). On the lattice, it is actually straightforward to measure a lattice approximation of the step scaling function. Indeed, we define the latter as: It is computed by measuring the renormalized coupling on lattices of size L/a and sL/a, at the same value of the bare coupling g 0 . The continuum step scaling function eq. (2.12) is then obtained by taking the continuum limit at a fixed value of the renormalized couplingḡ(µ), i.e., In the following section will apply this strategy to the GF couplings and discuss the determination of their β-functions at relatively high-energy scales. As we shall see, the approach to the perturbative asymptotic regime is dramatically slow for these schemes. This poses some severe limitations if one aims at extracting the Λ-parameter to high-precision using these coupling definitions.

Datasets
We measured the GF couplings on lattices with L/a = 8, 10, 12, 16, 20, 24, 32, 48, and for values of the bare coupling, β ∈ [6,11]. In total we collected between 1, 000 − 40, 000 measurements depending on the exact ensemble. The complete list of simulation parameters and corresponding results is given in Appendix E. This choice of parameters covers a range of renormalized couplings: and it allows us to determine the lattice step-scaling function Σ 2 (u, a/L) for L/a = 8, 10, 12, 16, 24, and Σ 3/2 (u, a/L) for L/a = 8, 16, 32. 8 Concerning the simulation algorithm, we used a combination of heatbath [4,48,49] and overrelaxation [50] as suggested in ref. [51]. In particular, we chose to alternate 1 heat-bath sweep with L/a over-relaxation sweeps. Since measuring the coupling (i.e. integrating the flow equations) is numerically more expensive than performing a Monte Carlo update, we repeated this process L/a times between subsequent measurements. In this way we took into account the expected O(a 2 ) scaling of the integrated auto-correlations of our observables. As a result, for basically all values of the simulation parameters, we obtained coupling measurements which are completely uncorrelated. The only exceptions are a few simulations for our largest lattices, L/a = 24 − 48, withḡ 2 GF ∼ 10; here the measured integrated autocorrelation times are τ int ∼ 1. In any case, we always take autocorrelations into account in our analysis through the Γ-method [52,53,54], implemented along the lines described in ref. [55]. We moreover note that in order to integrate the flow equations we use the adaptive step size integrator described in ref. [13]. This results in an improvement in computer time close to a factor 10 on our largest lattices compared to a fixed step-size integration of the flow equations.

The non-perturbative β-function at high-energy
In this section we study the viability of using the GF couplings to extract the Λ-parameter at highenergy. To this end, we first introduce a convenient high-energy scale, µ ref , by specifying a relatively small value for the GF couplings. Specifically, we define this scale in terms of the magnetic component of the coupling, and set:ḡ In the following we also need the corresponding value of the coupling in the electric scheme: . This is given in eq. (3.22), and we assume it known for the time being; we shall come back shortly to its determination. With these definitions at hand, the quantity we are interested to compute is: where GF may stand for either the magnetic or electric coupling scheme; clearly, I GF g is defined in terms of the proper β-function (cf. eq. (2.6)). Note that the results in the GF schemes are expressed in the MS scheme using the known relations between Λ-parameters, eq. (2.33).
To evaluate eq. (3.15) the necessary ingredient is the β-function in the range:ḡ 2 GF ∈ [0,ḡ 2 GF, ref ]. Our preferred strategy to obtain this is to consider a parametrization of the β-function of the form: where the coefficients b 0 , b 1 , b 2 are fixed to their perturbative values of eqs. (2.5), (2.32). This enforces the correct asymptotic behaviour of β(g) for g → 0. The coefficients p k are then determined by fitting our set of non-perturbative data. More precisely, we introduce the function: and define the χ 2 -function: is computed from the measured value of the GF coupling, u i , at a given i = (L/a, g 0 ), and from the corresponding result for the lattice step-scaling function Σ i s = Σ s (u i , a/L). In the above expression: parametrizes the cutoff effects in the data. Note that the data corresponding to different scale factors s can be combined into a single fit as long as we consider different parametrizations for the cutoff effects i.e., the coefficients ρ (s) k are independent parameters for different values of s. Concerning the nature of the discretization errors, as discussed in Sect. 3.1, our data is in principle affected by O(a) errors. Rather than including an explicit O(a) term in eq. (3.18), however, we decided to proceed in the following way. The data strongly supports the conclusion that within our statistical precision, the dominant discretization errors for our lattice resolutions and coupling values are of O(a 2 ) rather than O(a) (cf. Appendix A). We thus estimated through dedicated simulations the systematic effect on the value of the coupling induced by the deviation of the 2-loop result for c t (g 0 ), used in the simulations, from an educated guess for its actual, non-perturbative, value. We then added this systematic uncertainty to the measured values of the coupling, and performed the fits including this systematic effect. Appendix B contains a detailed discussion about this point.
Restricting to values of u ≤ḡ 2 GF, ref the dataset described in Sect. 3.3, the fits defined by eqs. (3.16)- (3.20) give in general excellent χ 2 's. More precisely, we consider fits with n b = 4, 5 in the parametrization of the β-function, and take n c = 2, 3 (n c = 2) to describe the cutoff effects of the data with s = 2 (s = 3/2). 9 The resulting fits all have a χ 2 /dof ∼ 0.5 − 0.9, where the fits to the electric components of the GF coupling always have smaller χ 2 's than those to the magnetic ones. The main reason for this is that the estimated O(a) uncertainties are significantly larger in the electric case, resulting in larger errors for the couplings and SSFs (cf. Appendix B).
In the case of the electric components of the GF coupling, in addition to the determination of the β-function, the evaluation of eq. (3.15) requires also the determination ofḡ 2 GF, e, ref =ḡ 2 GF, e (µ ref ). We can obtain the latter by performing a linear fit to our 19 data points forḡ 2 GF ∈ [2, 3], on lattices with L/a = 8, 12, 16, 24, 32, 48. Specifically, we fit the quantity, where a 0 , a 1 , ρ 0 , ρ 1 are fit parameters. The fit has a very good χ 2 ∼ 0.9, and allows us to obtain the precise result: This result is very stable under a change of the number of fit parameters, or of the data included in the fit. (Note that the data spans a factor 6 in the lattice spacing.) In practice, the tiny uncertainty that we obtain forḡ 2 GF, e, ref could be completely neglected in the following analysis, but we include it anyway.
The results for Λ MS /µ ref obtained according to the above strategy are summarized in table 3, and also presented in figure 1. Different parametrizations of the β-function and different data sets corresponding to different lattice discretizations of the GF couplings all produce consistent results.   There is also perfect agreement between the determinations from the electric and magnetic components. This is a highly non-trivial test, since the ratio of Λ-parameters in these schemes is ∼ 1.15. Only when converted to the common scheme Λ MS our results agree. As our final result for Λ MS /µ ref we quote the analysis based on the Zeuthen flow data for the electric components of the GF coupling with L/a > 10, and with n b = 5, n (s=2) c = 3 in the parametrization of the β-function. This gives: The central value of any other fit that uses the Zeuthen flow discretization with L/a > 8 is well included in this error band. The quoted uncertainty is rather conservative as we discard all data with L/a = 8, 10. For completeness, we also give the corresponding fit parameters for the β-function: . (3.25) Before moving to the next section, we want to mention that in order to gain further confidence in our determination of Λ MS /µ ref we also considered some alternative fitting strategy. One example is to determine the β-function entering eq. (3.15) by fitting our dataset according to the χ 2 -function (cf. eq. (3.18)): where the function G s (u i ) is the solution of the equation, with F defined in eq. (3.17). As before, in the expression for F we can simply consider a parametrization of the β-function as eq. (3.16). Alternatively, we can take: It is clear that in the limit x → 0 this parametrization reproduces the correct perturbative expression for the β-function. The nice feature of this choice is that we can now evaluate explicitly the function F (a, b) entering the χ 2 functions. Indeed, Once the function F is determined, it is straightforward to compute (3.15). Without entering into more details, very similar conclusions apply for these alternative fits as for our preferred strategy. The fits describe our dataset well, and different discretizations and schemes all give results for Λ MS /µ ref which are perfectly compatible with eq. (3.23).

Comparison with perturbation theory
It is instructive to plot our results for the β-functions and compare them with their 3-loop perturbative predictions. This is illustrated in figure 2. It is clear from the figure that the GF based schemes show a very poor convergence to their expected perturbative behaviour. For the electric scheme, the nonperturbative data is barely compatible within errors with the 3-loop β-function at the most perturbative point, where α ∼ 0.08. The magnetic scheme shows even poorer convergence, with a clear deviation at α ∼ 0.08 between our data and the 3-loop prediction. In this case, looking even beyond the range covered by the data, our parametrization for the β-function appears to deviate from its 3-loop approximation down to couplings as small as α ∼ 0.05. With hindsight, these conclusions are not too surprising, given the large value of the b 2 -coefficients in these schemes (cf. eq. (2.32)). The slightly better convergence of the electric scheme with respect to the magnetic one made us favour as our final result, eq. (3.23), a determination based on the electric rather than the magnetic coupling. Either way,  however, a determination of the Λ-parameter based on the GF schemes is very much limited in the attainable precision due to these issues. The effect of the poor convergence to perturbation theory of these schemes on the determination of Λ MS /µ ref can be quantitatively studied by examining the dependence of Λ MS /µ ref on the scale µ PT at which perturbation theory is used to extract it. More precisely, we can look at: where, as usual,ḡ GF, X ≡ḡ GF (µ X ), and GF may refer to either the electric or the magnetic components of the GF coupling. The function ϕ GF is defined by eq. (2.11) in terms of the GF β-function and its 3-loop perturbative approximation. As anticipated in the above equation, in the limit α PT → 0 the function φ GF approaches Λ MS /µ ref with O(α 2 PT ) corrections. The latter are expected to be small if and only if 3-loop perturbation theory is a good approximation for the β-function at scales of O(µ PT ). As stated above, however, this is not the case for the GF schemes, even at the largest energy scales we explored. At values of the couplings α PT = 0.2, for instance, where 3-loop perturbation theory is typically expected to be accurate, we find: The approximations φ GF, e , φ GF, m thus deviate from our final result eq. (3.23) by ∼ 6% and 13%, respectively. These numbers are between 3 to 6 times larger than our uncertainty on eq. (3.23).
Even at the largest energy scale reached by our simulations, where α PT ∼ 0.08, both schemes show a significant deviation from eq. (3.23) (cf. figure 3). As anticipated, the issue of poor convergence has a dramatic effect on the precision that can be attained for Λ MS /µ ref using GF based schemes. Due to the large corrections in eq.  .30)) as a function of the matching scale with perturbation theory (represented by α PT =ḡ 2 PT /4π). At α PT ∼ 0.08 (our most non-perturbative data point) the NNLO predictions using the electric and magnetic schemes disagree by ∼ 2.5%.
for instance, has an error 50% smaller than eq. (3.23). Unfortunately, however, at this value of α the results in different schemes differ by a few standard deviations. These values are also significantly different from φ GF (0), indicating that a safe contact with perturbation theory has not been reached. In conclusion, using GF based schemes, 3-loop perturbation theory is not accurate enough to extract the Λ-parameter with ∼ 1% error at α ∼ 0.1. The conclusions of this section are further corroborated by a more traditional step-scaling strategy to extract Λ MS /µ ref . We refer the interested reader to Appendix C for more details.

Non-perturbative matching to the SF scheme
Over the years the SF couplings have been widely employed in step-scaling studies [42,9,34,56,8], and recently their perturbative behaviour has been carefully investigated in the context of determining the Λ-parameter of QCD to a few per-cent accuracy. From these studies, we expect the matching of the SF schemes with perturbation theory to be quantitatively much better than what we observed above for the GF schemes. In particular, 3-loop perturbation theory may be precise enough for our purposes at scales where α ∼ 0.1. This in fact will be confirmed below.
Given these observations, in this section we consider a different strategy for determining Λ MS /µ ref . We will first match non-perturbatively the SF and the GF schemes, and then extract the Λ-parameter using perturbation theory in the SF schemes. To achieve this, we have measured the SF couplings on lattices of size L/a = 6, 8, 10, 12, 16, and at values of the bare coupling g 0 that match our measurements of the GF couplings on 2L/a lattices, i.e. L/a = 12, 16, 20, 24, 32, respectively. We thus combine a change of scheme with a change in the renormalization scale by a factor s = 2c, where c = 0.3 in our study. We collected from about 5 × 10 5 measurements on the smaller lattices with L/a = 6, up to 2 × 10 6 measurements for L/a = 16 (see table 17 for the exact figures). We measured both the SF couplingḡ 2 SF and the observablev, which gives us access toḡ 2 SF,ν for any value of ν (cf. eq. (2.18)).
In this section we focus exclusively on ν = 0. The consistency between determinations for different ν-values is discussed in Appendix C.
With these results at hand, we fit the data for the GF and SF couplings as: The functionρ(u) parametrizes the cutoff effects in this relation by a simple polynomial: We also use a polynomial for f (u): Using the perturbative relation, eq. (2.35), we could in principle fix the coefficients f 0 , f 1 , to their perturbative values. Here, however, we refrain to do so. Clearly, this implies that the non-perturbative matching between the SF and GF schemes is only trustworthy within the available range of couplings i.e., forḡ 2 GF, SF ∼ 1 − 2. We obtain good fits by taking, e.g., n f = 2, 3 and nρ = 2. We note however that our final results for Λ MS /µ ref do not depend significantly on this particular choice. We only see some dependence if we include or not our coarser lattices (i.e. L/a = 6 for the SF coupling and L/a = 12 for the GF couplings). Any disagreement disappears once we perturbatively improve the SF data to 2-loop order in lattice perturbation theory (see appendix D). Having observed this, in order to be conservative, we prefer to discard the coarser lattice spacing for computing the matching, and use the 2-loop improved data. All our fits are then good with a χ 2 /dof ∼ 0.5 − 1.
Once the function f (u) has been determined, the β-function in the SF scheme can be inferred from the one in the GF scheme using the relation: Note that since we are not interested in matching the GF schemes with their perturbative expressions, we can consider also fits for β GF (g) where the b 2 coefficient is treated as a fit parameter, rather than being fixed to its perturbative value (cf. eq. (3.16)). In this case, the results for β SF (g) only use as perturbative information the universal coefficients of the β-function, eq. (2.5). In particular the known value for b SF 2 is not used. We have determined β SF using the data for β GF both in the electric and magnetic GF scheme. The SF β-function determined this way is shown in figure 4b. As one can see from the plot, as expected there is agreement between the determination from the electric and magnetic GF schemes. The non-perturbative data then match the perturbative prediction in all the range of couplings we covered.
Similarly to what we did in the previous section for the GF coupling, we can now explore in the SF scheme the effect on Λ MS /µ ref of matching with perturbation theory at different energy scales, µ PT . To this end, we introduce the function (ḡ X, Y ≡ḡ X (µ Y )):  (cf. eq. (2.10)). This quantity is entirely analogous to eq. (3.30), the only difference being that once we arrive at the scale µ PT with the running in the GF scheme, we change to the SF scheme at the scale 2cµ PT , and we then use perturbation theory in the SF scheme (see figure 4a for a cartoon). Note that the factor s = 2c appearing on the r.h.s. of eq. (3.37) compensates the scale difference in matching the SF and GF schemes. The function φ SF has, of course, an asymptotic expansion of the form: This number has one of the largest uncertainties of all the analysis we considered, and its error covers all central values of the determinations which use only lattices with L/a > 8. In conclusion, thanks to the better perturbative behaviour of the SF schemes, by non-perturbatively matching the GF and SF schemes we can reliably quote an uncertainty significantly smaller than in eq. (3.23). Similar conclusions are obtained using a more conventional step-scaling strategy based on the GF couplings, and considering different values of ν for the SF scheme. We refer the interested reader to Appendix C.

Connection to an hadronic scale
To compute the Λ-parameter in units of an hadronic quantity, like t 0 [11] or r 0 [28], we must relate this quantity with the technical scale µ ref . In this section we consider two different strategies to compute this relation. Both rely on the determination of the β-function of the GF coupling at relatively lowenergy scales. In the first strategy, we first use the β-function to relate µ ref with a convenient low-energy finite-volume scale, µ had . We then relate, in a second step, this scale with the infinite volume scales t 0 and r 0 . We shall refer to this strategy as "fixed scale determination". In the second strategy, instead, we use the knowledge of the β-function to match directly µ ref with t 0 and r 0 ; we shall refer to this as "global analysis".

The β-function at low-energy
We begin by defining two convenient low-energy finite-volume scales, one for each GF scheme. We do so through the conditions:ḡ 2 GF, e (µ had,e ) = 10.95 =ḡ 2 GF, m (µ had,m ) ; (4.1) note that µ had,e = µ had,m . The desired ratios of scales µ ref /µ had can now be inferred once the β-function in the two schemes is known in the range of couplings: In this range, we find convenient to employ a parametrization of the β-function of the form: This is completely analogous to eq. (3.28), the difference being that here we do not fix the coefficients p 0 , p 1 , p 2 , to their perturbative values. As we saw for eq. (3.29), this parametrization allows us to express in a straightforward way the function F of eq. (3.17) in terms of the parameters of the βfunction, i.e.: As done previously for the β-function at high-energy, the fit coefficients p k can then be determined by minimizing the χ 2 -function, eq. (3.18) (or eq. (3.26)); we again consider the parametrization (3.20) for the cutoff effects. Once F is determined, the desired ratios of scales are given by: We consider fits to data withḡ 2 GF ∈ [2−4.7], and their corresponding SSFs. We obtain good fits choosing n b = 2, 3, and we typically take n c = 2, 3, although the final results are pretty much independent on the number of parameters used to parametrize the cutoff effects. Figure 6 collects a landscape of results. As one can see from the figure, we have excellent agreement between different analysis as long as we discard the coarser lattices with L/a = 8. As our final results we quote: To obtain the sought-after relation √ 8t 0 µ ref we now need to compute √ 8t 0 /a and aµ had for several values of the lattice spacing and extrapolate their product to the continuum limit, i.e.,     [58]. For r 0 /a, instead, the results are from † † [59] and [58]. The data labeled § is obtained from the values of r c /a of [60] together with r c /r 0 = 0.5133 (24) which gives the quoted values of r 0 /a. This is done by determining the finite volume and hadronic scale in units of the lattice spacing at matching values of the bare coupling, β = 6/g 2 0 , which guarantees that the value of a is the same up to scaling violations. The quantity t 0 /a 2 is known from the literature over a wide range of β. In table 4 we collected the results used in our study and the corresponding references. Results for t 0 /a 2 are known for values of the lattice spacing as fine as 0.025 fm. Simulations at very fine lattice spacings deal with the problem of topology freezing [45] either by simulating very large physical volumes [57], or by using open boundary conditions in the Euclidean time direction [58].
In order to determine aµ had , instead, for a given set of lattice sizes L/a, we must find the values of the bare coupling, β = β had (L/a), for which the conditions (4.1) are satisfied. At these β's we then have: aµ had = a/(cL) (cf. eqs. (2.29),(2.30)). Given our large set of data we can consider lattices with sizes L/a ∈ [12,32]. 10 The values of β had (L/a) are then easily found by performing interpolations of the data in β at fixed L/a. The values of β had (L/a) determined this way are collected in table 5, where the results for different lattice discretizations of the GF couplings are given.
From tables 4 and 5, it is clear that the values of β for which the large volume quantities t 0 /a 2 are available do not match exactly those which correspond toḡ 2 GF (µ had ) = 10.95 for our L/a's. To obtain the products (4.7) at matching β-values we can thus proceed in either of the following two ways: 1. Given the results of table 4, we fit the data for log(t 0 /a 2 ) as a function of g 2 0 . In practice, we obtain a good description of the data by using a simple polynomial fit of degree 5. In this way we can infer the values of log(t 0 /a 2 ) at the bare couplings of table 5, corresponding to the finite volume results for aµ had .
2. We fit the results for β had (L/a) as a function of L/a, or equivalently aµ had . In this case, we obtain good fits using polynomials of degree 3. Given this functional form, we can then find the value of aµ had at the bare couplings where the large volume quantities t 0 /a 2 have been computed (cf. table 4).

Determination of √ 8t 0 µ ref : global analysis
As can be seen from table 4, the flow scale in lattice units t 0 /a 2 has been determined precisely even at very fine lattice spacings. This data, together with our results for the β-function at low energy in Sect. 4.1, calls for an alternative strategy to determine directly the product √ 8t 0 µ ref . 10 In principle data with L/a = 48 is also available. However, fixing the value of the coupling toḡ 2 GF = 10.95 requires in this case an interpolation of data points which are significantly distant from the target value (cf. table 9). We thus prefer not to include this data in our analysis. Nonetheless, we have checked the effect of including it and found complete agreement for our final result, but with a slightly smaller error.      We start by performing once again a fit to the data for log(t 0 /a 2 ) of table 4 by a polynomial in β = 6/g 2 0 of degree 5. This gives us a parametrization: in the range β ∈ [5.96, 7.01]. We then consider our GF coupling data in the very same range of β, for lattice sizes L/a = 16, 20, 24, 32, 48. For the ease of presentation we focus on the results for the magnetic components of the coupling; the analysis using the electric scheme gives compatible results. We will denote these generically as u FV =ḡ 2 GF,m (µ FV ; a/L, g 0 ). Using the parametrization of the β-function, eq. (4.2), valid in this range of couplings, we can compute the scale factor between the reference scale µ ref and the scale µ FV corresponding to a given coupling u FV =ḡ 2 GF (µ FV ). This is simply given by (cf. eq. (4.3)): where u ref ≡ḡ 2 GF (µ ref ), of eq. (3.14). Note that the scale factor associated to a given u FV is only given up to discretization errors, the reason being that we are considering the value of the coupling at a certain L/a and g 0 . Given the results for R, we recall that aµ FV = a/(cL). The parametrization (4.10) then allows us to determine the second matching factor: which relates the given finite-volume scale aµ FV with our hadronic quantity in lattice units. Combining finally the factor R of eq. (4.11) with Q of (4.12) we obtain: which, up to discretization errors, is the product √ 8t 0 µ ref we were after. Note that by construction the dependence on u FV cancels in the leading term of the product, eq. (4.13): only discretization errors thus depend on u FV . The latter turn out to be relatively small. This can be appreciated in figure  9, which shows the different results for the factor S t 0 . As we can see from the figure, the product S t 0 is approximately constant as a function of u FV with variations of at most 2% when changing u FV ∼ [3.5, 11.5]. This observation suggests to fit our results for S t 0 (u FV , a/L, g 2 0 ) as: where a 0 and ρ k are fitting parameters. The continuum result √ 8t 0 µ ref , corresponding to the coefficient a 0 of the fit, shows very little dependence on the degree of the polynomial describing the cutoff effects, as long as n c > 1. As our final estimate we quote: which is based on the Zeuthen-flow data in the magnetic scheme for L/a ≥ 16 and uses n c = 2. This result is in agreement with our previous result, eq. (4.9). We note that including also lattices with L/a < 16 gives completely compatible results with eq. (4.15), but with smaller uncertainties.

Determination of r 0 µ had
The strategies we used to compute √ 8t 0 µ had can also be applied for the determination of r 0 µ had . In this case, however, the situation is a little complicated by a few technical issues. First of all, the different results for r 0 /a available in the literature which we could consider [58,60,59] use different discretizations for the relevant observables. Secondly, two of these determinations, refs. [60,59], are more than 20 years old. A note of concern in this case is thus the issue of topology freezing [45], which was not known at the time of these computations. This issue is potentially more significant at the very fine lattice spacings simulated in ref. [60], while the results of ref. [59] are confined to substantially larger spacings. The more recent determination of r 0 /a in [58], on the other hand, employs open boundary conditions. In addition, their values for r 0 /a cover a factor two in lattice spacings and reach down to small spacings comparable to those of [60] (cf. table 4). For these reasons we consider safer to restrict our attention in the following only to the results of [58,59].
Most of the results of [58] have been obtained at values of the bare coupling g 2 0 which lie outside the range of bare couplings of table 5 (cf. table 4 where the results of [58] are labelled with a ). This means that the fixed scale strategy considered in Sect. 4.2 cannot really be applied to this data set. A fixed scale determination can only be considered for the results of [59] (which are labelled by † † in table 4). For these we opt for strategy 2 of Section 4.2. We thus use the same fit for β had (L/a) as a function of L/a considered there, which is based on a 3rd degree polynomial. Given this functional form, we then find the value of aµ had = a/(cL) at the values of g 2 0 where r 0 /a is known, and compute r 0 µ had . The continuum extrapolation lim a→0 r 0 a × aµ had = r 0 µ had , (4.16) corresponding to the r 0 /a data of ref. [59] are shown in figure 10. Using the results of these extrapolations and those of eq.  The global approach described in Section 4.3 can on the other hand also be used for the r 0 /a data of ref. [58]. We thus begin by fitting the data for r 0 /a as a function of the bare coupling g 2 0 . We do this separately for the data set of [58] (labelled [B] in figure 11) and for the data set of [59] previously used (labelled [A] in figure 11); this because the two computations use different observable discretizations.  [58] in the range of β covered by the finite volume simulations (table 5) are also plotted. These however are not used for any continuum extrapolation (see main text for the details).
With these functional forms at hand, we then determine the product of the running factor µ ref /µ FV (cf. eq. (4.11)) and the matching factor r 0 µ FV (cf. eq. (4.12)). The result: is expected to be constant, up to scaling violations. Figure 11 shows the quantity S r 0 for the two data sets. It is clear that scaling violations within each data set are very small. Looking at the continuum extrapolated values, obtained after fitting separately the two data sets to a constant with cutoff effects parametrized by a 1st degree polynomial (cf. eq. (4.14)), we find:

The Λ-parameter
We are now ready to express Λ MS in terms of the hadronic scales t 0 and r 0 . Our first main result is: The first factor on the r.h.s. of this equation comes from the non-perturbative high-energy determination of the GF β-function in the electric scheme, combined with a non-perturbative matching of the GF and SF schemes. Thanks to this strategy, perturbation theory in the SF scheme could be safely used to run to infinite energy and obtain the result of eq. (3.39). The second factor derives instead from  Figure 11: The product S r 0 (u FV , a/L, g 2 0 ) of the running factor R and the matching factor Q r 0 (analogous to eq. (4.12) for r 0 ). This has to be constant and equal to r 0 µ ref up to discretization errors; the plot clearly shows that these are very small for both data sets. In different colours we have the two data sets: [A] is ref [59], [B] is ref [58] (cf. table 4), while different symbols correspond to different lattice sizes. The grey error band is the result of eq. (4.18) (see main text for more details).
our determination of the β-function at low-energies in the magnetic GF scheme (cf. eq. (4.5)). The third and last factor is the non-perturbative matching of the magnetic GF coupling with the gradient flow scale t 0 , eq. (4.8). The result in eq. (5.1) would not practically change, if we were to choose the scale µ had,e defined in the electric scheme in the second and third factor. We would also obtain a perfectly compatible result ( √ 8t 0 Λ MS = 0.6227(94)), if we replaced the last two factors with √ 8t 0 µ ref determined via the global approach, eq. (4.15).
Our final error estimate (∼ 1.6%) is very conservative. All three factors are determined by dropping the two coarsest lattice spacings at our disposal, and we typically chose as final result the analysis with the largest uncertainty. The only exception was when we favoured the result of eq. (3.39) over eq. (3.23). We found however compelling to match with the asymptotic perturbative regime in the SF rather than the GF scheme. Using the GF scheme would have indeed increased our final error by more than 50%, due to the bad convergence of its perturbation theory; this despite of the fact that we had access to the non-perturbative running of the coupling up to very large energy scales, whereḡ 2 GF ∼ 1. Looking at the different contributions to our final uncertainty, it is clear from eq. (3.39) that most of the uncertainty comes from the determination of the non-perturbative running from α = 0.2 to α ∼ 0.08. Given the high-precision we aimed for, however, we found compulsory to reach, nonperturbatively, these high-energy scales, and accurately test the applicability of perturbation theory. Figure 12 then illustrates the separate contribution to the total error squared from the different simulations and some other sources. Most of our error comes from the most expensive simulations, while for instance the contribution from the uncertainty on t 0 is completely negligible. Also the systematic uncertainty deriving from our ignorance of the boundary O(a) counterterm (labelled as "c t effect" in the figure) contributes only little to the final uncertainty. Moreover, in our strategy the SF coupling is only used at very high energies, where it is measured very precisely and with negligible cutoff effects. Consequently, its contribution to the final error is very small. In summary, our final uncertainty is completely dominated by the statistical uncertainty in the measurements of the GF coupling on our finest/largest lattices. The results could thus be improved significantly by just investing more computer time on these simulations.    [61]. These are: ALPHA 98 [62], QCDSF/UKQCD 05 [63], Brambilla [64], Flow 15 [65], Kitazawa 16 [66], Ishikawa 17 [67].
A completely analogous analysis leads to:  Figure 13 shows a comparison of our computations with some other determinations available in the literature. More precisely, we have included the results entering the last FLAG average [61]. We find a significant discrepancy with some of these determinations, in particular with those of QCDSF/UKQCD 05 [63], Flow 15 [65] and Kitazawa 16 [66]. These determinations extract the Λ-parameter from measurements of Wilson loops, and rely on 2-loop bare lattice perturbation theory at a scale of a few GeVs. Here we recall that our final results use continuum extrapolations performed with data that cover a factor two in lattice spacings, with two extra lattice spacings (reaching a factor 3 change in a) used to check the consistency of our results. The non-perturbative running is performed up to very high energy scales where α < 0.1, and our matching with perturbation theory has been performed in several renormalization schemes. Our determination satisfies the most stringent of the FLAG criteria.
Nevertheless there is a significant discrepancy with the FLAG average. Although we think that the FLAG criteria are conservative, this work shows that even computations that meet these criteria can differ by more than the quoted uncertainties. This discrepancies need further investigation, and once clarified the criteria used to rate different lattice determinations of the strong coupling might need to be readjusted. Also the authors believe that given the level of precision reached by current computations, one should probably consider √ 8t 0 Λ MS rather than r 0 Λ MS as the standard quantity of comparison.

Conclusions
Renormalization schemes based on the gradient flow have many attractive properties for applications in lattice QCD. Renormalized couplings are defined via observables with very small variance, which allows to attain great statistical precision. In addition, the coupling is given directly by an expectation value. Hence, there are no systematic errors associated with extracting properties from the large Euclidean time behaviour of a correlator. By using a renormalization scheme based on the gradient flow we have obtained a rather precise determination of the Λ-parameter in the SU (3) gauge theory. Along the way we have learned several important lessons that will come useful when applying this technology to the more relevant case of QCD.
Using finite size scaling techniques we have determined non-perturbatively the β-function in a range of couplings α ∼ 0.1 − 1. Our results indicate that at energy scales where α ∼ 0.1 contact with 3-loop perturbation theory is not safe for GF schemes if one aims at a precision ∼ 0.5% in α s .
One might argue that our particular setup (finite volume renormalization schemes with Schrödinger functional boundary conditions in the pure gauge theory), can cast some doubts on the general validity of our conclusions. In this respect one should first note that we have checked two different renormalization schemes (based on the electric and magnetic components of the energy density at positive flow time, respectively). The scheme based on the electric components, in particular, has a very similar perturbative behaviour to the corresponding infinite volume scheme (i.e. their non-universal 3-loop coefficients of the β-function are pretty close; see e.g. figure 2a). Secondly, if one considers the parametric uncertainty originating from the missing perturbative orders in extracting α s from the infinite-volume GF coupling in QCD one reaches similar conclusions to our non-perturbative study: the extraction of α s at the electroweak scale in the GF scheme carries a 0.5% theoretical uncertainty. If the extraction is performed at a few GeV (the energy scale typically accessible to large volume simulations), the theoretical uncertainty increases to ∼ 2 − 3%. Quark effects are absent in our study, but perturbatively their effect at high energy is small compared to the effect of the gluons [17]. The presence of quarks is thus expected not to change the picture very much. We conclude that the qualitative conclusions of our study are indeed general.
Our work allows us to precisely determine the pure gauge Λ-parameter in units of a typical hadronic scale (we considered both t 0 and r 0 ). Our result for r 0 Λ MS , with a precision ∼ 1.7%, shows some tension with other recent lattice computations, in particular with those where the MS-coupling is extracted from Wilson loops at an energy scale set by the lattice cutoff. One drawback of the GF couplings are the relatively large cutoff effects, which have been observed in many different applications (see [15] for an overview). Despite the fact that we have a solid theoretical understanding of the nature of these cutoff effects [16], they are still the main source of concern when considering GF-based observables. In order to have discretization effects under control we used 5 different lattice resolutions which cover a factor 3 in the spacing, and two different discretizations to integrate the flow equations and compute our observables. We see no deviation from pure O(a 2 ) scaling violations. Despite these observations, we quote results where the two coarsest lattice spacings are discarded, and we add a generous estimate for the O(a) boundary effects. We recall that we have performed the running non-perturbatively up to very large energy scales, we have matched with the perturbative behaviour in four different renormalization schemes (with their respective Λ-parameters varying by factors of two), and used at least two different methods to match with a large volume hadronic scale. All in all, the significant discrepancy with the results in the literature shows the difficulty in extracting the fundamental parameters of the Standard Model with high precision.
We conclude by pointing out that the results of this work represent a serious warning for any attempt of reducing the current uncertainty of the world average value of α s using lattice QCD and the GF schemes, especially if one aims at an infinite volume determination. Here the range of scales that can be explored is limited to α 0.25, completely insufficient to quote sub-percent precision in α s . On the positive side we have shown a viable strategy to reach a precision of 0.3% in α s . It combines the use of the GF schemes to determine the running non-perturbatively, and a non-perturbative matching at highenergy with the traditional SF schemes, that show small effects in the truncation of the perturbative series. Such a project would also require a precise and accurate determination of the hadronic scale in a theory with three or more active quarks [68].  . (A.2)). (Data from table 6.) As we can see, the slope is positive at weak couplings, while turns negative at strong couplings.

A Continuum limit of GF couplings
Several studies have shown that renormalized couplings defined through the gradient flow are affected by significant cutoff effects (see e.g. ref. [15] and references therein). These cutoff effects have been first carefully studied at tree-level of perturbation theory [69], and later more systematically in the context of Symanzik' effective theory [16]. Despite these efforts, however, they remain the main source of concern in applications of running couplings derived from the gradient flow.
In order to establish that the lattice resolutions employed in our study are fine enough to obtain accurate continuum extrapolations, we here analyse in detail the continuum extrapolation of the lattice step scaling function, Σ 2 (u, a/L), for the magnetic GF scheme, at 3 representative values of the coupling; these are: u 1 = 1.04784, u 2 = 3.5705, Note that we shall focus on our preferred choice of discretization for the observable, i.e., Zeuthen flow/improved definition (cf. Sect. 3.1). To perform the continuum limit of the lattice SSF at these coupling values, we first need to find the values of the bare coupling β for the chosen lattice sizes, for which the renormalized coupling, u =ḡ 2 GF, m (µ), is equal to the target values (A.1). We did this for lattice sizes: L/a = 8, 10, 12, 16, 24. The results of this tuning are given in the third and fourth columns of table 6. Once the values of the bare coupling are determined, at these β-values we compute the GF coupling on lattices of double size, i.e. with L/a = 16, 20, 24, 32, 48, respectively. The corresponding results are given in the last column of table 6. Finally, the lattice SSFs so obtained can be extrapolated to the continuum limit (see figure 14a).
We consider for Σ 2 continuum extrapolations linear in (a/L) 2 . 11 Moreover, in order to gain some insight on higher-order discretization errors we consider linear extrapolations in (a/L) 2 of 1/Σ 2 . We expect these two strategies to give compatible results up to O((a/L) 4 ) errors. We thus have: As a further test on the scaling properties of our data, we also study the effect of discarding our coarser lattice with L/a = 8 from the extrapolations.  Table 6: Data for the lattice step scaling function, Σ 2 (u, a/L), of the magnetic GF coupling, at the target couplings of eq. (A.1). The results of different continuum extrapolations are given.
Having this noticed, all our fits have good quality, and the agreement among different determinations of the continuum SSF, σ 2 (u), is in fact quite good (compare the rows marked by L/a = ∞ in table 6). 12 We thus conclude that for our choice of discretization, our data show no significant deviation from O(a 2 ) scaling. In addition, figure 14a shows that the slope of the continuum extrapolations is positive at weak couplings, while changes to negative at strong couplings. Somewhere around u ∼ 4, the data have no significant cutoff effects.
In summary, the detailed study presented in this appendix shows that once lattice sizes in the range L/a = 8 − 24 are considered, within the whole range of couplings u ∈ [1, 5] the continuum extrapolations of the lattice step scaling function of the GF coupling present no significant deviation from O(a 2 ) scaling within our precision.

B Boundary O(a) effects
It is well-known that due to the breaking of translational invariance in the time direction, even the pure Yang-Mills theory with Schrödinger functional boundary conditions suffers from O(a) discretizations effects [6]. These can in principle be entirely removed by a proper tuning of a single boundary counterterm coefficient, c t (g 0 ) (cf. eq. (3.1)). Unfortunately, however, in practice there is no compelling method to determine c t non-perturbatively. As a result, at present, given our choice of Wilson gauge action, c t is only known in perturbation theory to two-loop order [32,33]. The leading discretization errors in our data are thus parametrically of O(g 6 0 a/L). On the other hand, the results of the investigation in Appendix A, where we ignored any O(a) effects in the data, support the conclusion that these effects are in practice very small, and below our statistical precision. Nonetheless, in order to guarantee the high-precision of our results, we here want to address this potential source of systematic effects in detail. To this end, we measured the GF couplings for several values of c t which are shifted from the 2-loop result, c t , used in the simulations: We did this for L/a = 8, 10, 12, and at the β-values corresponding to the 3 couplings of eq. (A.1).
The results we obtained are collected in table 7. The deviation, ∆ḡ 2 GF, e/m , of these results from the couplings of eq. (A.1), is a clear measure of the sensitivity of the coupling on the boundary improvement coefficient c t . The data of table 7 is well represented by a fit: whereḡ 2 GF, ct is the GF coupling (either electric or magnetic) measured for a given c t , and ∆c t = c t −c t . For the fit coefficients a 0 , a 1 we obtain the results: As expected from general considerations, the electric components are the most affected by boundary O(a) effects [15]. The values of these fit coefficients, on the other hand, are basically the same for our two choices of discretization of the observable, i.e., Wilson flow/clover and Zeuthen flow/improved observable.
Having established the sensitivity of the GF couplings around the two-loop value of c t , eq. (B.1), in order to estimate the uncertainty to attribute to our data for the incomplete tuning of c t , we now need an estimate for the difference between c t and the non-perturbative value of c t . Having no information about the latter, a reasonable guess is to take for this deviation the full two-loop term of the series eq. (B.1). Given the fact that the coefficients of the series (B.1) appear to decrease with the order, our estimate can be considered a conservative one.
In conclusions, we add in quadrature to the statistical uncertainty of the GF couplings computed at c t = c t and for a given g 0 and L/a, the systematic uncertainty: with a 0 , a 1 given by eqs. (B.3). We stress once again that this is done in order to take into account possible O(a) effects in our data that might arise from the mistuning of the boundary conterterm coefficient c t . With the exception of our coarsest lattices, this effect turns out to be sub-dominant to the statistical errors.  (48) Table 7: Data to study the c t dependence of the GF couplings. The columns labelled by ∆ḡ 2 GF, e/m give the deviation of the GF couplings measured at c t = c t + ∆c t , from the values of eq. (A.1) measured at c t = c t .

C
Step-scaling determination of Λ MS /µ ref As an alternative determination of Λ MS /µ ref we also consider a more traditional approach based on the computation of the continuum step-scaling function of the GF coupling, combined with a nonperturbative matching of the GF and SF couplings. For the determination of the continuum SSF we employ the very same data for Σ 2 (u, a/L) entering the computation of the corresponding β-function (cf. Sect. 3.4). 13 More precisely, we consider data for the Zeuthen flow/improved observable discretization, we study both the electric and magnetic definition, and restrict our attention to couplings: u ≤ḡ 2 GF, ref , whereḡ 2 GF, ref ≡ḡ 2 GF, ref,e/m , depending on the chosen scheme. 14 The lattice SSF, Σ 2 (u, a/L), is then fitted according to the functional form: where σ 2 is a parametrization of the continuum SSF: with coefficients s 0 , s 1 , s 2 fixed to their perturbative values (cf. eqs. (2.5),(2.32)), while the function ρ (2) models the coupling dependence of the leading discretization errors in the data: (We recall that we assume the leading discretization errors to be O(a 2 ) as the O(a) boundary effects are taken into account as systematic uncertainties (cf. Sect. B).) The data with L/a > 8 is well described by the above functional form for any combination of n σ = 6, 7 and n c = 2, 3, giving a χ 2 /dof ∼ 0.5 − 1, depending on the exact fit. In particular fits to the electric GF coupling data always have smaller χ 2 than those involving the magnetic ones. Similarly to what we did for the β-function, in order to be conservative, we take as our preferred fits, fits to the data with L/a > 10. We then choose n σ = 6 and n c = 3. These fits yet have excellent χ 2 's, but have larger errors. The second step is the determination of the non-perturbative matching between the GF and SF couplings. More precisely, we here consider three different SF coupling definitions, corresponding to ν = −0.3, 0, 0.3. Similarly to what discussed in Sect. 3.5 for the ν = 0 case, fits of the form (3.33) give a very good description of the data; moreover also in this case the results for Λ MS /µ ref depend very little on the exact choice we make. We thus settle on fits where f 0 and f 1 are fixed to their perturbative values, n f = 3, and nρ = 2. Note that, in being once again conservative, we neglect the L/a = 6 SF and L/a = 12 GF coupling results in determining the matching, even though these are well described by the fit function too. We then opt for not considering the PT improvement of the SF data as this has anyway no significant effect once L/a ≥ 8.
Having all the basic ingredients, we can now proceed with the determination of Λ MS /µ ref . Starting from the value of the GF coupling, u 0 =ḡ 2 GF, ref , the knowledge of the continuum SSF allows us to infer the value of the coupling at the energies scales µ n = 2 n µ ref , n = 1, 2, . . ., by solving the recursion relation (ḡ 2 GF ≡ḡ 2 GF, e/m , σ 2 ≡ σ GF, e/m 2 ): Considering both the range of GF couplings we covered and the range where the non-perturbative matching with the SF couplings is available, our data allow us to perform n = 6 steps, i.e., starting from µ ref , we are able to increase the energy scale by a factor 64. For a given value of u n =ḡ 2 GF, e/m (µ n ) determined this way, we can now compute Λ MS /µ ref through: where either X = GF,ḡ X,n =ḡ GF, e/m (µ n ) and s n = 2 n , or X = SF,ḡ X,n =ḡ SF (2cµ n ) and s n = 2 n+1 c.
The value ofḡ SF (2cµ n ) is of course inferred from that ofḡ GF, e/m (µ n ), using the non-perturbative matching relation between the SF and GF schemes previously established. We then expect that the results for different values of n and/or schemes should all agree, up to O(ḡ 4 X,n ) corrections, asḡ X,n → 0. In Fig. 15 we present the results based on the SSF of the electric GF scheme. As one can see from the figure, the determination of Λ MS /µ ref obtained directly from the GF coupling suffers from large O(α 2 ) corrections. The situation is of course completely analogous to what we already discussed in Sect. 3.4 in terms of the β-function (cf. Fig. 3). In particular, if we read off the value of Λ MS /µ ref for say, n = 6, the result is: Although this is compatible with our final estimate, eq. (3.23), it is clear that given the trend of the results for different values of n, a reliable determination of both mean value and error, necessarily requires an extrapolation to α = 0. Quite different is the situation for the determinations obtained after switching non-perturbatively to the SF couplings. These show indeed much milder O(α 2 ) corrections, particularly so for the case ν = 0 (cf. Fig. 4c). At values of the coupling α ∼ 0.1, any discrepancy between different ν-values is negligible within our statistical errors, and we can thus safely quote: which corresponds to the value obtained for ν = 0 at n = 5. This is in perfect agreement with our final estimate eq. (3.39). For completeness we also present in Fig. 16 the corresponding results based on the magnetic GF scheme. It is no surprise that the results for Λ MS /µ ref coming from the magnetic coupling show larger O(α 2 ) corrections than what we have seen for the electric scheme (cf. Fig. 3). Considering for instance the results for n = 6, we have: Λ MS /µ ref = 0.0830 (10), which is clearly biased in both central value and error if compared to eq. (3.23) (note that α GF, m (µ n=6 ) ∼ 0.08!). On the other hand, once again switching non-perturbatively to the SF schemes solves the issue. Taking for instance the results for ν = 0 and n = 5 we obtain: Λ MS /µ ref = 0.0798 (9), well in agreement with the result of eqs. (3.39) and (C.8). In conclusion, the determination based on the SSF of the GF couplings rather than the β-function gives perfectly compatible results, reinforcing even further the robustness of our analysis and conclusions.

D Perturbative improvement of the SF coupling
In this appendix we give some details on the perturbative improvement of the SF coupling in the ν = 0 scheme, employed in Sect. 3.5.
In bare lattice perturbation theory, the SF coupling is given by, where the coefficients of this series can be written as: Referring to ref. [32] one can easily show that and which approximation should be good enough in practice. The coefficients m a 1 , m a 2 and m b 2 instead can be found in table 1 of ref. [32].
To work out the cutoff effects in the SF coupling to two-loop order we need to know the asymptotic behavior of m 1 and m 2 for L/a → ∞. This is given by (see table 2  where b 0 and b 1 are the universal one-and two-loop coefficients of the β-function, eq. (2.5). Given these results we can define, The two are of course equivalent up to O(g 8 0 ) terms. In Sect. 3.5 we considered the form eq. (D.12). Note that by using eq. (D.1) and the results of Sect. 2.3, we could in principle rexpress the perturbative improvement of the SF coupling in terms the GF couplings, rather than the bare one. Although this might appear more natural when determining the matching between the GF and SF couplings, we prefer not to do so and stick with the bare results. In any case, as already mentioned in the main text, the effect of the perturbative improvement is very small in practice (cf. table 8). Choosing a different option hence does not make any real difference. 15 Note that a slightly different value for m ∞ 2 was first given in ref. [32].           (47)