A non-perturbative exploration of the high energy regime in $N_\text{f}=3$ QCD

Using continuum extrapolated lattice data we trace a family of running couplings in three-flavour QCD over a large range of scales from about 4 to 128 GeV. The scale is set by the finite space time volume so that recursive finite size techniques can be applied, and Schr\"odinger functional (SF) boundary conditions enable direct simulations in the chiral limit. Compared to earlier studies we have improved on both statistical and systematic errors. Using the SF coupling to implicitly define a reference scale $1/L_0\approx 4$ GeV through $\bar{g}^2(L_0) =2.012$, we quote $L_0 \Lambda^{N_{\rm f}=3}_{\overline{\rm MS}} =0.0791(21)$. This error is dominated by statistics; in particular, the remnant perturbative uncertainty is negligible and very well controlled, by connecting to infinite renormalization scale from different scales $2^n/L_0$ for $n=0,1,\ldots,5$. An intermediate step in this connection may involve any member of a one-parameter family of SF couplings. This provides an excellent opportunity for tests of perturbation theory some of which have been published in a letter [1]. The results indicate that for our target precision of 3 per cent in $L_0 \Lambda^{N_{\rm f}=3}_{\overline{\rm MS}}$, a reliable estimate of the truncation error requires non-perturbative data for a sufficiently large range of values of $\alpha_s=\bar{g}^2/(4\pi)$. In the present work we reach this precision by studying scales that vary by a factor $2^5= 32$, reaching down to $\alpha_s\approx 0.1$. We here provide the details of our analysis and an extended discussion.


Introduction
The Standard Model seems to describe all high energy physics experiments carried out to date, in some cases with extraordinary accuracy (cf.[2] for the most recent PDG review).For processes involving the strong interactions the precision is usually less impressive, due to our limited ability to extract quantitative information from QCD.One of the main tools is perturbation theory (PT) in the strong coupling, α s , and there has been significant progress in high order perturbative QCD calculations, with renormalization group functions now available up to 5-loop order in the MS-scheme [3][4][5][6][7].However, before a perturbative result can be confronted with experimental observables, the transition from quarks and gluons to hadronic degrees of freedom needs to be modelled in some way.Such models come in various shapes and forms, from "hadronisation Monte Carlo" in jet physics to "quark hadron duality" in QCD sum rules.A common problem then consists in assigning systematic errors to the model assumptions.A further issue is the reliability of PT itself, given that the series is only asymptotic.To some extent, the reliability can be assessed within PT itself, by comparing different orders of the expansion, or by increasing the energy scale, µ, such that α s (µ) becomes small, due to asymptotic freedom.Unfortunately, the rapidly increasing complexity of higher order calculations means that typically only a few terms in the perturbative series are available.In addition, the energy scale is often defined by the kinematics of the physical process under consideration.The variation of the scale is then rather limited and to assign an error to the perturbative result is difficult 1 .
In this work we carry out a systematic investigation into the reliability of PT.We do this by directly comparing non-perturbative QCD observables to their perturbative expansions, over a wide range of scales.Lattice QCD, together with a careful treatment of the continuum limit, is currently the only way to obtain such non-perturbative results, subject only to standard assumptions such as locality and universality.The main reason why this is rarely done is the usual limitation of any numerical approach: on a finite system it is very expensive to simultaneously resolve very different length scales.Most lattice QCD projects aim at hadronic low energy physics, and the space-time volume, L 4 , must then measure several femto metres across in order not to distort the hadronic states of interest.At least for massive single particle states, the finite volume effects are exponentially suppressed [9] and one may then pretend to be in infinite space time volume, up to a systematic error which is often below the percent level.On the other hand, with current lattice resolutions of, say, L/a < 100 this means that the cutoff scale set by the inverse lattice spacing, 1/a, reaches a few GeV at most, and the deep perturbative high energy regime seems out of reach.It is important to realize that this limitation is only due to the requirement that the lattice covers a physically large space-time volume.If this constraint is lifted, there is nothing that prevents simulations at very high energies, albeit in physically tiny space-time volumes.The observables2 we consider in this situation are all normalized as effective couplings, which run with L, the scale set by the finite space-time volume.In order to achieve this we set all quark masses to zero and scale all other dimensionful parameters proportionally to L, thereby obtaining a mass-independent scheme.In the high-energy regime, PT can be used to relate to more commonly used  The combination shown here yields directly the lowest order cofficient, b 0 of the β-function as (σ(u) − u)/u 2 = 2b 0 ln 2 + O(u).The dashed lines show the perturbative 2-loop behavior.The purple 1-sigma band shows our result (fit C in table 5).Data points for N f = 0, 2, 3, 4 are taken from the literature [18][19][20][21].
schemes such as the MS scheme of dimensional regularization.Moreover, by combining the idea of a finite volume scheme for the coupling with recursive step-scaling techniques [10,11], one may both determine the scale L in units of some hadronic scale, and reach the perturbative high energy regime without ever requiring enormous lattice resolutions, L/a.Obviously, the finite space-time volume constitutes an integral part in the definition of these observables.PT must then be adapted to this situation.While the Euclidean spacetime signature used in lattice QCD is advantageous in PT, all the sophisticated tools of standard PT in (infinitely extended) momentum space are of limited use.
As part of the project to determine α s (m Z ) from low energy hadronic input in 3-flavour QCD [1,12,13], our collaboration has applied these techniques to a 1-parameter family of finite volume couplings in Schrödinger functional (SF) schemes, for which the 3-loop βfunction is known [14][15][16][17].We have measured these couplings in numerical simulations and for a range of lattice sizes with unprecedented precision.Extrapolation to the continuum limit of this data allows us to carry out stringent tests of renormalized perturbation theory for energy scales ranging from about 4 GeV to 128 GeV.A first account of our results has appeared in a letter [1] and we here provide the details of this work and a more extended analysis.
The technique, used earlier for between N f = 0 and N f = 4 quark flavours [18][19][20][21], allows one to non-perturbatively verify the close-to perturbative running of the coupling and observe the small effects of dynamical quarks, as illustrated in figure 1.A preview of our final result is included in the figure, demonstrating our advanced precision.
The paper is organized as follows: section 2 uses a continuum language to explain how our QCD observables are defined and collects the relevant perturbative results from the literature.We also comment on "non-perturbative effects" which are associated with secondary minima of the action.Section 3 then presents the lattice set-up, the numerical simulations and statistics produced, and discusses the perturbative improvement of the data.The impatient reader might skip this section and directly pass to section 4. There, after the discussion of the continuum extrapolated results and associated systematic errors, the comparison to renormalized perturbation theory is performed before we conclude in Section 5. Finally, a technical appendix presents the models we used for the sensitivity of the data to a variation of the two O(a) boundary counterterm coefficients c t and ct .

SF couplings
In order to apply the recursive step-scaling techniques to lattice QCD, it is desirable to define renormalized QCD couplings in a finite space-time volume, L 4 , and in the chiral limit.Such finite volume renormalization schemes are quark mass independent by construction [22], and the renormalization scale is set by µ = 1/L.It is then possible to apply recursive finite size scaling methods and trace the scale evolution over a wide range without the need for very large lattice sizes, L/a [10].Still, these requirements leave many options, such as the boundary conditions for the fields and the exact choice of observable.We here choose Schrödinger functional boundary conditions [23,24]: these introduce a gap in the spectrum of the Dirac operator, so that numerical simulations can be performed directly at zero quark masses, without the need for any chiral extrapolation.Moreover, perturbation theory remains tractable in this framework, as the absolute minimum of the action is unique up to gauge equivalence.For the observable we choose the traditional SF coupling [25,26] and a 1-parameter family of close relatives [27].The most important reason for this choice is the existence of a 2-loop calculation in this case [14,15], which, in combination with [16,17] allows to infer the 3-loop β-function for these schemes.Furthermore, the values of the 3-loop β-function coefficients are reasonable and enable us to make contact with the asymptotic perturbative regime at energy scales in the range O(10-100) GeV.
In the future one might also consider the more recent coupling definitions based on the gradient flow [28,29].The QCD 3-loop β-function is currently known in the case of infinite space-time volume [30], and there is progress for the case of a finite volume with SF boundary conditions [29] using numerical stochastic perturbation theory [31][32][33].These results seem to point to a 3-loop β-function coefficient which is significantly larger than in the MSand SF-schemes.This indicates that gradient flow couplings may not be ideal for matching with the asymptotic perturbative regime.Furthermore, cutoff effects are typically larger with the GF couplings than with the traditional SF coupling [13], so that larger lattice sizes are required.This partially offsets other computational advantages.Obviously, further studies are required and one should re-assess the situation once more perturbative information becomes available.

SF ν schemes
In the continuum the Schrödinger functional is defined as the Euclidean path integral, with the Euclidean continuum action S = S g + S f , Here, g 0 denotes the bare coupling constant, F µν is the field tensor associated with the gauge field A µ , and is the covariant derivative acting on the quark fields.It includes a constant U(1) background field which we set to θ µ = (1 − δ µ0 )θ, with the choice θ = π/5.
In the spatial directions L-periodic boundary conditions are imposed on all fields.At the time boundaries the fermionic fields satisfy [24] with the projectors For the gauge field one has with the boundary values C k and C k .The boundary condition at x 0 = 0 refers to the gauge transformed field, The integration over the SU(3)-valued and spatially periodic gauge functions Λ(x) in Eq. (2.1) ensures gauge invariance of the Schrödinger functional.The spatially periodic Λ(x) fall into different topological sectors labelled by an integer n, which is related to the topological charge of the gauge field, ) provided the Chern-Simons action of the boundary gauge fields C k , C k vanishes (which is the case for the choice below).The value of the gauge action in each sector n is then subject to the usual instanton bound [23] Using the gauge invariance of the Schrödinger functional under the transformations, Λ(x) → Ω(0, x)Λ(x) , (2.12) one may convert the integral over gauge functions Λ to a sum over n, with Λ in Eq. (2.6) replaced by fixed representatives Λ n for each topological sector.In particular one often sets Λ 0 = 1.
We now focus on Abelian and spatially constant boundary gauge fields, with traceless and diagonal 3 × 3-matrices φ and φ .Their diagonal elements ) still depend on 2 real parameters, η and ν.In the temporal gauge and the topological charge zero sector the field equations with these boundary conditions are solved by, which corresponds to a constant chromo-electric field, Inserting the field tensor into the gauge action, S g , one obtains which, for given η (and independently of ν) constitutes the absolute minimum of the action [23].One may thus define the effective action as a function of this background field, and its perturbative expansion, . The SF couplings ḡ2 ν (L) can be defined through In fact the ν-dependence is explicit, since both 1/ḡ 2 (L) and v(L) are ν-independent.In terms of the effective action, v(L) reads Note that the ν-independence of Γ 0 [B], implies that v(L) has a perturbative expansion starting at O(1).This ensures the correct normalization of the whole 1-parameter family of couplings, namely ḡ2 ν = g 2 0 to lowest order.Finally we remark that the entire 1-parameter family is determined by the expectation values, defined in terms of the functional integral, Eq. (2.1), at ν = 0. Once the lattice regularization is in place both quantities will thus become observables in numerical simulations.

β-functions and perturbative relations to the MS-coupling
The SF couplings are defined independently of perturbation theory and thus the same is true for their β-functions, where the asymptotic expansion on the r.h.s.starts out with the standard universal coefficients b 0,1 for and the 3-loop coefficient is given by (2.28) The 3-loop coefficient has been obtained by matching the coupling at the 2-loop level to the MS-scheme, where the β-function is now even known to 5-loop order (b 3 and b 4 in our notation) [3][4][5][6][7].For later use we collect the numerical values for N f = 3 QCD in table 1, together with the SF ν scheme results for the 3 choices of the parameter, ν = −0.5, 0, 0.3, which we selected for more detailed analysis in Section 4. Comparing the MS to the SF ν scheme we note that the respective 3-loop β-functions coincide for ν ≈ −0.3.In general, ν-values of O(1) are reasonable from a perturbative point of view.
Closely related to the β-functions are the step-scaling functions which connect couplings at scales which differ by a factor 2. Defining the precise relationship is, √ and the perturbative expansion of σ(u), is thus determined in terms of the coefficients of the β-function, with the first 3 given by Finally, we quote the relation between the SF and the MS couplings, in terms of α = ḡ2 /(4π) at the scales µ = 1/L and sµ, respectively, with s > 0. One finds In order to connect to the SF ν couplings for ν = 0 we need the expansion of v in the coupling ḡ.Defining the expansion is known to second order, where the coefficients for N f = 3 evaluate to 4πv 1 = 1.797887(5) , (4π) 2 v 2 = −0.741(14) . (2.38) we obtain the 2-loop relation, Inverting perturbatively and combining with the previous equations we have where In the perturbative matching of couplings one occasionally applies the principle of "fastest apparent convergence", which implies that s = s is chosen such as to make the one-loop coefficient, c ν 1 (s ), vanish.This is the case for and with this choice one obtains the relation, (2.45)

Perturbation theory and the Λ-parameter
There are various ways to define a target precision for α s .Instead of referring to the coupling in some scheme at some scale it is attractive to instead refer to the Λ-parameter.
Given the coupling ḡx (L) in a scheme x the Λ-parameter in this scheme is a special solution of the Callan-Symanzik equation of the form with Note that this definition is independent of perturbation theory provided the coupling and its β-function are defined non-perturbatively.In practice, however, one would like to evaluate the Λ-parameter at a large energy scale µ = 1/L such that the integral in the exponent can be safely evaluated in perturbation theory.The exact scheme-dependence of the Λparameter is obtained by the one-loop matching of the respective couplings.Labelling the schemes by x and y, ḡ2 one obtains the exact relationship Note that this allows one to indirectly define Λ MS non-perturbatively, thereby justifying its use as a reference definition.With the perturbative matching coefficients of the previous subsection we obtain the relationships (for where Λ and Λ ν are the parameters for the SF and SF ν scheme, respectively.In particular, the ratio s of scales used in Eq. (2.44) is given by the ratio of the respective Λ-parameters.

On exponentially suppressed corrections to perturbation theory
The perturbative expansion of the path integral generates an asymptotic series, with zero radius of convergence.In applications one then hopes that, for the accessible range of couplings, the perturbative series provides a good quantitative description of the observable.The observables we consider here, the couplings in the SF ν schemes, are defined non-perturbatively in Euclidean space-time, with an infrared cutoff provided by the finite space-time volume.These properties are advantageous for perturbation theory, in particular, there should be no renormalon problem [34].Lattice QCD provides very good non-perturbative control of these observables, for couplings α in the range 0.1−0.2(cf.Section 3).Before testing perturbation theory, we would like to identify exponentially suppressed terms in the coupling which might preclude a good quantitative description of the non-perturbative data.Such terms are associated with local minima of the action, e.g.those corresponding to the classical solutions of the field equations.Given the instanton bound, Eq. (2.10), and the absolute minimum S g [B] = 2π 2 /g 2 0 of the action (Eq.(2.19), with η = 0), contributions from the |Q| = 1 instanton sector to our observables are accompanied by a suppression factor exp(−6π 2 /g 2 0 ) = exp(−3π/(2α)) and are therefore numerically irrelevant for our range of couplings.We may then ask the question whether there are further secondary minima of the action which are less strongly suppressed.Hence we are looking for a secondary minimum B * µ of the gauge action in the Q = 0 sector, which satisfies (2.51) In fact there are "large" gauge transformations at x 0 = 0 corresponding to gauge functions ω(x) which are topologically trivial but are not subject to the gauge fixing procedure around B µ [23].In order to find potential secondary minima we have resorted to a numerical experiment in the lattice discretized theory.More precisely, we have first performed numerical simulations of the pure SU(3) Yang-Mills theory on a lattice with linear extent L/a = 8, (at β ≡ 6/g 2 0 = 5.7), and generated a long Monte Carlo history of about 64.000 configurations, corresponding to 128.000 MDU, using the same simulation code as for our N f = 3 QCD simulations (cf.Section 3).Every 5 th gauge configuration has been taken as initial condition for the gradient flow equation [28], which we then integrated up to very large flow times t, corresponding to c = √ 8t/L = 10; the gradient flow is a smoothing operation and drives the gauge field towards a local minimum of the action.At large flow times we selected the gauge field configurations in the Q = 0 sector3 .Apart from the background field, Eq. (2.17), we have indeed found a single further local minimum.In order to check for its stability and to obtain its continuum limit, we have performed similar simulations on finer lattices with L/a = 12, 16, 24, and bare couplings such as to keep ḡ2 (L/2) = 2.77 approximately fixed.After extrapolation to the continuum and in the temporal gauge we find that this secondary minimum corresponds to the spatially constant Abelian field, (2.54) The boundary conditions at x 0 = 0 thus are given as and (2.56) The gauge function ω(x) is thus non-constant in the x 3 -direction, which induces the shift by ±2π in 2 of the angles of C ω 3 , in addition to the permutation of the colour 2-and 3-components of φ, Eq. (2.16).Obviously the spatial directions can be permuted, so this minimum has a 3-fold degeneracy.Hence, the classical field B * µ is Abelian and spatially constant, but with boundary values, transformed by the gauge function To find the gap in the gauge action we insert the non-zero components of the field tensor into the gauge action Eq.(2.2), with the result Hence the gap, ∆S, is found to be 10/3 in units of π 2 /g 2 0 which is 2/3 below the Q = 1 instanton threshold.This leads to a suppression factor exp(−g 2 0 ∆S/(4πα)) = exp(−5π/(6α)).For the range of couplings in our study, this factor varies from a few times 10 −6 to below 10 −10 , which renders such a non-perturbative contribution completely negligible.

Lattice action
We choose the standard Wilson plaquette action for the gauge fields and three flavours of non-perturbatively O(a) improved Wilson fermions.The lattice action is then given by S = S g + S f , with The gauge field action S g is a sum over all oriented plaquettes p on the lattice, with the weights w(p), and the parallel transporters U (p) around p. With the gauge field boundary conditions given in terms of the Abelian fields, Eq. (2.15), the gauge part of the action is completely specified by setting w(p) = 1 except for timelike plaquettes touching one of the boundaries for which w(p) = c t .The Dirichlet boundary conditions for the quark fields look exactly the same as in the continuum, cf.Section 2.
Like in the continuum we take the fermionic fields to be spatially periodic and implement the phase θ = π/5 via a constant U(1) background field the Wilson-Dirac operator in the fermionic action (3.2) takes the form, which includes the Sheikholeslami-Wohlert term [35].For the clover leaf definition of the field strength tensor, F µν , we refer to [36] and the improvement coefficient c sw (g 0 ) is set non-perturbatively using the result from [37].Finally, the fermionic O(a) boundary counterterm action is specified by [36] The 2 boundary counterterm coefficients, c t (g 0 ) and ct (g 0 ) are set to their perturbative two-and one-loop expressions, respectively [15,25], with the known perturbative coefficients for N = 3 colours given by c (1) t We notice a significant cancellation in the one-loop term c (1) t between the gluon and fermion contributions.We interpret the resulting relative size of one-and two-loop terms for N f = 3 as an accident and not a sign for a poor behaviour of the series in general.

Lattice observables
Like in the continuum, the basic observables 1/ḡ 2 and v are given as expectation values, Eq. (2.25), of gauge invariant fields, which are now obtained as ηand ν-derivatives of the lattice action 4 (3.1,3.2).The lattice version of the Abelian background field takes the form, with B µ (x) the continuum expression, Eq. (2.17).Cutoff effects with such Abelian gauge fields are known to be small [23].Indeed, the η-derivative of S g [V ] yields the lattice renormalization constant which converges to 12π with O(a 4 ) corrections.We will use this lattice definition of k in order to ensure ḡ2 = g 2 0 exactly at lowest order.Note that this also holds for ḡ2 ν , since v vanishes identically at tree level.
On the lattice with Wilson quarks, the chiral limit is not sharply defined, and one also needs to specify the exact definition used.For given bare coupling g 0 , we require the PCAC quark mass, to vanish on an (L/a) 4 lattice with the Abelian boundary conditions, Eq. (3.3).Here f A (x 0 ) and f P (x 0 ) are Schrödinger functional correlation functions defined e.g. in Eqs.(2.1) and (2.2) of [38], and ∂ 0 , ∂ * 0 are the forward and backward lattice time derivatives, respectively.Finally, the improvement coefficient, c A , is set to its perturbative 1-loop value [39,40], since a non-perturbative estimate is not available for N f = 3 and our choice of gauge action.Given that we do not attempt to reach the low energy, hadronic regime, we expect one-loop perturbation theory to work reasonably well for c A .The chiral limit is now defined by m(L) = 0, and, for given bare coupling g 2 0 ≡ 6/β, the bare mass m 0 for which this equation holds, defines the critical mass parameter or, equivalently, the critical κ, With these conventions we may now define the lattice observables.Specifying the value u of the coupling ḡ2 (L) at vanishing quark mass m(L) defines our approach to the continuum limit, and other lattice observables are then well-defined functions of u.In particular v gives rise to 2 lattice observables which differ by the chiral limit definition.The appearance of 2 lattice versions for ω(u) is a consequence of the definition of the lattice step-scaling functions through which requires simulations on lattices with resolutions L/a and 2L/a, at the same bare parameters.In particular, the simulations on the 2L/a-lattices are performed at the bare mass parameters for which the PCAC mass vanishes on the L/a lattice.Finally, we also consider the lattice step-scaling functions for ḡ2 ν , at non-zero values of ν.

Perturbatively improved lattice observables
In order to accelerate the approach to the continuum limit one may use perturbation theory to subtract the lattice artefacts order by order in the coupling from the non-perturbative data [41].The 2-loop calculation in [15] has been carried out in the very same lattice regularized theory, and the two-loop lattice artefacts in the ν = 0 step-scaling functions, are indeed available to this order.With the coefficients for N f = 3 from table 2, one may thus define the improved step-scaling functions, up to loop order i = 2.By construction, the leading cutoff effects for i = 0, 1, 2 are then given by5 and are thus suppressed by additional powers of the coupling.The term linear in a/L is due to the incomplete cancellation of the O(a) boundary effects and could be eliminated by a non-perturbative determination of c t and ct .We will come back to the question of remnant O(a) effects in Subsect.3.7.
For the observables Ω and Ω one parametrizes the cutoff effects by 2 functions, and ˜ .For Ω we have with perturbative expansion and analogous equations hold for Ω and ˜ .Unfortunately, the published results of the 2-loop calculation do not allow for the extraction of the cutoff effects for this case, so that the perturbatively improved observables, ) are only available to 1-loop order, i = 1, with the coefficients 1 and ˜ 1 given in table 2.
The same remark applies to the step-scaling function Σ ν for ν = 0. Using the notation, the one-loop coefficient is given by where v 1 is the expansion coefficient of the continuum function ω(u), Eq. (2.38).Values for δ ν 1 can be inferred from table 2, for N f = 3 and the lattice sizes relevant for this study.Values of the coefficients for N f = 3 and the relevant lattice sizes, as required for perturbative cancellation of lattice artefacts up to 2-loop order in Σ, and to one-loop order in Σ ν , Ω and Ω, cf.text.

Simulation parameters and statistics
Using the openQCD code [42,43] we have simulated lattice sizes L/a = 4, 6, 8, 10, 12 around 9 values of the coupling ḡ2 (L) = u in the range 1.1 − 2.0, cf.table 3.At the same bare coupling g 2 0 = 6/β and bare quark mass am 0 = 1/(2κ) − 4 we then doubled the lattice sizes and simulated for 2L/a = 8, 12, 16 and, in 3 cases also for 2L/a = 24, cf.table 4. Starting from the L/a = 12 lattices we have tried to approximately match the values of the coupling for ν = 0 at L/a = 4, 6, 8, so as to be able to do continuum extrapolations of the step-scaling function at individual values of the coupling, without the necessity for large interpolations of the data.
As a target precision we chose the criterion, which is reached for most of our data except for some L/a = 10 lattices.These lattices were however not used for the step scaling procedure as we did not generate corresponding  configurations on 2L/a = 20 lattices.Except for some checks we also refrained from using lattices as small as L/a = 4 and thus do not list the results here.However, the L/a = 10 data and the 2L/a = 8 data are used for the continuum extrapolation of Ω and Ω, respectively, and are therefore included in the tables.
Note that the choice of the reference value ν 0 = 0.3 is rather arbitrary.In fact, the data in the table for ḡ2 , ḡ2 ν 0 =0.3 and v, with their statistical errors enables the reconstruction of the coupling at any value of ν, using Eq.(2.39) and straightforward error propagation, We have checked that this reconstruction does indeed reproduce the result of a direct data analysis at a given ν-value, provided that the treatment of autocorrelations is done consistently for the couplings at all ν-values and v.We find that the precision for the ν = 0 coupling, Eq. (3.30), translates to higher values for other choices of ν, for instance we find an increase of 20 percent for ν = 0.3 (from tables 3 and 4), and ca.50 percent for ν = −0.5 from Eq. (3.31).All statistical errors were determined using the Γ-method [44].For our observables, one even has to be careful that one sums up the autocorrelation function sufficiently far.Still the final autocorrelation times range from values somewhat below 2 MDU for weak coupling and small L/a, to about 8 MDU at larger coupling and L/a = 24.Further details on the performance of our algorithms will be reported in [45].

Treatment of statistical errors
When forming the step scaling function Σ(u, a/L) there are statistical uncertainties both for ḡ2 (L), table 3, and for ḡ2 (2L), table 4.These are propagated to the error of Σ(u, a/L) with u the central value of the estimate of ḡ2 (L), via To estimate the required derivative ∂Σ/∂u we differentiate the 3-loop truncation of the continuum function, σ(u), Eq. (2.31), corrected for the known lattice artefacts at one-and two-loop order for ν = 0 and ν = 0, respectively, cf.Subsect.3.3.For ν = 0 this leads to and similarly for ν = 0 with δ ν 1 from Eq. (3.29), the unknown δ ν 2 set to zero and with the scheme dependence of s 2 (via b 2 , Eq. (2.32)), taken into account.As a cross check, we also estimated the derivative directly from the data and found the differences to be negligible.
For the study of the observables Ω and Ω we proceed similarly: to obtain the derivative with respect to u we first perform a rough continuum extrapolation neglecting the errors on u.The resulting polynomial fit function ω(u) ≈ 0.14307 − 0.004693 × u + 0.0077906 × u 2 − 0.0105266 × u 3 + 0.0023996 × u 4 , (3.34) is then differentiated to provide an estimate for ∂Ω/∂u and ∂ Ω/∂u, neglecting any L/adependence of the derivative.

Quality of tuning to the chiral limit
An important aspect of Wilson fermions is the need to tune the bare quark mass parameter (parameterized by κ) to a critical value, such that chiral symmetry is restored up to cutoff effects.For our choice of condition m(L) = 0, with the PCAC mass of Eq. (3.15), we have performed extensive tuning runs which enable a precision such that, at all stages of the calculation [45].The corresponding values for κ are given in table 3. What is the tolerance of a slight mistuning of the mass?Using 1-loop perturbative results from ref. [26] for the mass dependence of ḡ2 and v we obtain, in the continuum limit, This should be compared with the target statistical precision, which is, for ν = 0, given in Eq. (3.30).We follow ref.[19] and allow for an uncertainty of about 1/3 of the statistical error.Neglecting small cutoff effects in the mass derivative and for N f = 3 this yields the bounds, for the ν-values that we chose for more detailed analysis in Sect. 4. We note that the achieved precision of the mass tuning, Eq. (3.35), stays well within these bounds, by at least a factor 10.Even if these perturbative estimates turned out to be significantly off the mark, e.g. by a factor 2, the systematic error associated with imperfect quark mass tuning would still be negligibly small and can thus be safely ignored.

Lattice artefacts linear in a/L
Despite the use of a non-perturbatively O(a) improved bulk action the very presence of the time boundaries in the Schrödinger functional creates lattice artefacts linear in a.
In principle these could be cancelled by an appropriate non-perturbative tuning of the improvement coefficients c t and ct , Eqs. (3.1, 3.7).In practice, however, we are currently limited to the use of perturbative estimates, Eqs.(3.8, 3.9).Hence some remnant linear aeffects in our data cannot be excluded.Instead of including a corresponding term in the fit ansatz for the continuum extrapolations we try to estimate the size of these uncertainties and include them as an additional systematic error.Using a combination of simulations and perturbation theory we have produced a model for the sensitivity of our data to a variation of c t and ct .The details are deferred to appendix A, where we obtain linearized shifts of the data, for instance, and analogously for a shift c t = ct + ∆c t .Hence, the model yields an estimate of the data that would have been obtained if the simulations had been performed at slightly different values c t and c t .To complete the model we thus need an educated guess for ∆c t (g 0 ) and ∆c t (g 0 ) such that the difference between a fully non-perturbative definition of c t and ct and the perturbative estimates (3.8, 3.9) is likely to be covered.We here choose i.e. a term of the neglected order with an effective coefficient.In the case of c t which is known to 2-loop order, cf.Subsect.3.1, we use a geometric progression and define (2) We note that particularly the choice for ∆c t is likely an overestimate, due to the accidental cancellation of the gluonic and fermionic terms observed in Subsect.3.1.
There are several options for the inclusion of this systematic error.We chose to proceed as follows: we first perform continuum extrapolations ignoring potential O(a) errors in both the original and the shifted data.We then take the spread of a given observable as an additional systematic error and add it in quadrature.Obviously this assumes that this systematic error is subdominant.We have therefore dismissed all continuum extrapolations where this turned out not to be the case.We will discuss the impact of these variations on the continuum extrapolations in the next section.

Continuum extrapolation of the step-scaling function
We now proceed with the continuum extrapolation of the data for the step-scaling function, for our default scheme with ν = 0.The 19 available data points for lattice resolutions L/a = 6, 8, 12 are shown in figure 2. Simulation parameters have been chosen such as to have approximately matched u-values between different L/a, and this is seen in the vertical line-up of the data.The fact that the data are so close together at given u-value illustrates that cutoff effects in the SF scheme with the chosen lattice regularization are generally small, even without perturbative improvement.While our data enables a more traditional  5).The data points are the approximations at finite L/a = 6, 8, 12 taken from table 4 with errors from Eq. (3.32).
continuum extrapolation, u-value by u-value, we have done this only as a cross-check.Our preferred strategy is to simultaneously fit all data to a global ansatz of the form (4.1) Here i = 1, 2 denotes the order of perturbative improvement of Σ and i = 0 refers to unimproved data.In general, such global fits have the advantage that an interpolation of the data to common u-values is not required.More importantly, however, the expected smooth u-dependence of the step-scaling function both on the lattice and in the continuum limit, is automatically built into this ansatz.As anticipated in the last section, we assume leading cutoff effects to start at O(a 2 ), with the linear a-effects being treated as systematic errors.Our fit ansätze for the cutoff effects thus are of the form, and the assumption of no lattice artefacts, ρ (i) = 0, is referred to by n ρ = 0.For the continuum step scaling function we consider polynomial fits with n c = 2 parameters, or 1-parameter fits (n c = 1), where s 0,1,2 are fixed to their perturbative values Eqs.(2.32).As the lattice artefacts are generally small at most n ρ = 2 parameters are required to obtain excellent fits to the data.
A selection of our fits is given in table 5.As an example we consider a 4-parameter fit (fit D) with n c = n ρ = 2 to the 2-loop improved data at ν = 0, Note that the fit coefficient c 1 is not far from the perturbative value s 2 = 0.001151; it is therefore reasonable to fix this parameter to the perturbative one and only fit a next order coefficient.Hence the majority of fits in

16
Table 5: Overview of the continuum fit functions and results.The naming convention is the same as in ref. [1].The two errors in the fit parameters are the statistical and the total error respectively, where the total error includes the systematic uncertainty from a variation of c t and ct , added in quadrature.
Given the smallness of the cutoff effects, even fit G with n ρ = 0 parameters seems reasonable, if one restricts to data with L/a ≥ 8.For the 2 continuum fit parameters of fit G the results are, c 1 = 0.0006 (12), c 2 = 0.0011 (7), Cov(c 1 , c 2 ) = −0.86 × 10 −6 .(4.7) While the χ 2 /d.o.f = 13/9 = 1.44 does not look too good, a comparison with fits B and F (with n ρ = 1) indicates that this may be an accident.In fact the χ 2 -values are not a sharp criterion in our case, as these strictly refer only to the statistical errors of the data and the given fit functions used, and thus do not account for the systematic uncertainties from cutoff effects linear in a.
In order to quantify these systematic uncertainties we repeat the fits with the data shifted by varying either c t or ct , as explained in Subsect.3.7.For fits with a single continuum parameter, n c = 1, we then take the spread in central values for this parameter as a systematic uncertainties due to either c t or ct variations and combine them in quadrature with the statistical error to obtain a total error of the fit parameter.Thus, in table 5, the fits with n c = 1 show 2 errors, the first being the statistical and the second the total error.In all fits we find that the c t -uncertainty dominates the effect of the ct -uncertainty; for instance, for fit B we obtain (25) stat.(15) ∆ct (6) where the r.h.s.takes the form given in table 5.For fits with n c = 2 continuum parameters we proceed in the same way.However, rather than quoting a total error on the continuum fit parameters, we propagate these uncertainties to the observables in table 6, where the results from the n c = 2 fits D and G are given with both a statistical and total error.While the total errors for most fits are dominated by the statistical error, this is not the case of fit G, where the total errors are predominantly systematic, cf.table 6.This indicates that fits with n ρ = 0 are too rigid to account for the O(a) variation of the data.While n ρ = 1 fits B and F are acceptable, we settled for fit ansätze with n ρ = 2 and n c = 1 to data with L/a ≥ 6 as our preferred choice (fits A,B,C,E,H).Then, using the 2-loop improved data leaves us with fits C and E, which are essentially equivalent, and figure 2 shows σ(u) from fit C with its error band.

4.2
The SF coupling for ν = 0 at scales L n = L 0 /2 n We now use the continuum fit functions for the step-scaling function at ν = 0 to evaluate the coupling at different scales L n = L 0 /2 n , separated by factors of 2. Our starting point is the reference scale L 0 , defined implicitly by The value 2.012 corresponds to the largest value of the coupling u for which the stepscaling function is known.In physical units the scale L 0 has been determined to be around 1/(4 GeV) [12].We note that σ(2.012) defines the coupling ḡ2 (2L 0 ), so that the lowest energy scale reached with the SF coupling is around 2 GeV.Recursive application of the continuum step scaling function, σ(u), allows us to obtain, in the continuum limit, the couplings at This defines the couplings u n as a set of observables, with our data enabling the recursion up to n = 5, thereby covering a total scale factor of L −1 /L 5 = 2 Table 6: Results for the couplings u n = ḡ2 ν (L n ), the Λ-parameter evaluated at u n , cf.Eq. (4.15), in units of the reference scale, L 0 (4.9), and the effective β-function coefficient, b eff 3 (4.14),for most fits of table 5. Results for L 0 Λ obtained with fits E, F and H are given in table 7.

Effective and fitted β-function
Given σ(u) in terms of 1 or 2 continuum parameters c k , one may translate this result into an effective 3-loop coefficient of the continuum β-function.For convenience we define b(g 2 ) = −gβ(g) so that Then Eq. (2.30) becomes, which can be solved for b eff 3 , with the result, Note that b eff 3 will depend on the value u where it is measured.Extracting this coefficient at different values of u should yield consistent results in the perturbative regime, and this is indeed the case for the ν = 0 data, cf.table 6.
This motivates a different parameterization of our fits with a single continuum parameter, namely via a 4-loop coefficient b fit 3 in the β-function as a fit parameter 7 .This is the purpose of fits E, F and H, cf.table 5, where we have taken σ(u) to be defined by Eq. (4.12) with b(u) = b 3loop (u) + b fit 3 u 5 and inserted σ(u) into Eq.(4.1).The resulting values for the fit parameter b fit 3 are given in table 5.This representation of our continuum results is very practical.While the fit function in Eq. (4.4) allows us to find the couplings at scales which are separated by a factor 2, the β-function readily yields the scale ratio separating two given couplings.

Determination of the Λ-parameter
Once the coupling u n = ḡ2 (L n ) is small enough, it is justified to use three-loop perturbation theory for the β-function in the expression and determine the Λ-parameter in units of L n and thus in units of L 0 = 2 n L n .Note that the expansion of the integral in the exponent is unknown at order ḡ4 as this term requires the knowledge of the 4-loop coefficient b 3 which is not available in the SF scheme.Provided such higher order terms are small, the result for L 0 Λ should be independent of n and the way the integral is evaluated.For completeness we note that our default evaluation consists in the direct numerical integration, using the truncated 3-loop β-function without expansion of the integrand or the exponential function.
The results for Λ in units of L 0 are given in table 6, where Eq. (4.15) is evaluated for the coupling at scales L n , for n = 0, . . ., 5 and for the various fit functions.An alternative evaluation of the Λ-parameter is obtained with the fits E, F and H in terms of a fitted β-function.One simply inserts the β-function into Eq.(4.15) and evaluates the integral numerically between ḡ2 (L 0 ) = 2.012 and ḡ2 (0) = 0.The resulting Λ-parameters are given in table 7 and show a remarkable consistency.We will discuss the results further in Subsect. 4  The continuum extrapolation for Ω(u, a/L) and Ω(u, a/L) proceeds along the same line as for the step-scaling function.A difference is that both data sets can be constrained to the same continuum limit but require separate fit coefficients for the cutoff effects.Moreover, the lattice resolutions L/a cover the range 6 − 24, i.e. a factor of 4 in scale and thus allow for an excellent control of the continuum limit.
The global fit ansätze used here are and analogously for Ω(i) with ρ(i) .Here, i = 1, 0 refers to 1-loop improved data (cf.Subsect.3.3) or unimproved data, respectively.In the models for the cutoff effects we just include 2 quadratic terms in a/L for either data set, with coefficients ρ 1,2 and ρ1,2 , e.g.
and the powers of u are chosen according to the expectation from perturbation theory.As in the case of the step-scaling function, linear terms in a/L will be treated as systematic errors.
with fit parameters d k , k = 1, . . ., 4 and v 1 and v 2 set to the known perturbative coefficients, Eq. (2.38).We have also experimented with separate fits to Ω (i) and Ω(i) and find good overall consistency.Here, we restrict the discussion to combined fits of the Ω (i) and Ω(i) data, with a common continuum fit function, ω(u).We distinguish fits of type A and B with 3 and 4 continuum fit parameters, respectively.Hence, fits of type A have 3+2×2 = 7 parameters, while type B fits have 8 parameters.
With these fit ansätze one obtains decent χ 2 /d.o.f.values for the one-loop improved data, even when including all 52 data points with L/a ≥ 6 (cf.tables 3 and 4).Given this much data we may afford to exclude the L/a = 6 lattices, thereby reducing the number of data points to 44 .An example for the continuum function ω(u) thus obtained is ω(u)| fit A,i=1,L/a≥8 = 0.14307 − 0.004693u + 0.01284u 2 − 0.01480u 3 + 0.003349u 4 .(4.20) The fit has a χ 2 /d.o.f.= 33.5/37and the covariance matrix for the fit parameters is given by Note that the error encoded in the covariance matrix is only the statistical error.To account for the systematic effect estimated from the variation of the O(a) counterterm coefficients c t and ct (cf.Subsect.3.7), we here proceed in complete analogy with the analysis of the step-scaling function.In table 8 we quote 2 errors, the first statistical, the second including the effect of a c t and ct -variation.This only marginally increases the errors, as is evident from table 8.The fits to the unimproved data have higher χ 2 /d.o.f.values, emphasized in bold face in table 8, unless the L/a = 6 data are dropped.As mentioned above, χ 2 is not the full story, given that our fits assume the absence of a/L effects and this effect is taken into account afterwards by our c t , ct -variation.However, we do see that 1) these variations have a tiny effect on the continuum values and 2) still, for example, ω(1.11) of the large χ 2 fits is off significantly.These fits have to be discarded.The other ones, which cover a remarkable range of lattice spacings, are entirely consistent.These observations allow us to conclude that perturbative improvement works very well in our coupling range, our treatment of c t , ct -variations is safe (maybe overly conservative), and most importantly, resolutions a/L ≤ 1/6 are sufficient to apply our continuum one-loop two-loop Ω (1) (u, a/L) − ρ (1) (u, a/L) Ω(1) (u, a/L) − ρ(1) (u, a/L) Figure 3: The bands show the continuum fit functions for fits of type A and B to one-loop improved data for L/a ≥ 8 and the data points are one-loop improved data with the cutoff effects subtracted using the models ρ (1) (u, a/L) and ρ(1) (u, a/L) from the type A fit.
extrapolations which assume that O((a/L) 3 ) effects have a negligible effect.All this makes us very confident also in the continuum extrapolations of Σ, where the very small lattice spacings are not available, but where we have 2-loop perturbative improvement at our disposal.
We return to the specific discussion of ω.As our best value at the reference coupling u = 2.012 = ḡ2 (L 0 ) we choose the result of fit A to 1-loop improved data with L/a ≥ 8. ω(2.012) = 0.1199 (10), (4.22) which is required to define the starting point for the step-scaling procedure for SF ν schemes with non-zero values of ν (s.below).Another interesting value is ω(u) at the largest available coupling, u = 2.45, which correspond to α = 0.195, ω(2.45) = 0.1117( 13), (4.23) using the same fit function.As discussed further in [1], and as is evident from the large difference between 2-loop PT and the non-perturbative result in Figure 3, an unnaturally large next order perturbative coefficient would be required to perturbatively describe the function ω(u) at such values of the coupling.Finally, we comment on the different behaviour of the fits A and B, which is seen in figure 3 for small couplings, outside the range of the data.This illustrates the danger of using fit functions outside their range of validity.While fit A is constrained to produce the 2-loop perturbative result for ω(u), Eq. (2.37), fit B leaves the 2-loop coefficient v 2 as a fit parameter, d 1 (4.19).The result, (4π) 2 d 1 = −0.9(2.9), should be compared with Eq. (2.38).While the central value is not too far off, the large error illustrates the difficulty to estimate such asymptotic coefficients, even if precise data is available over a wide range of couplings.in α 2 is very pronounced for ν = −0.5.Therefore, it is a strong consistency check for our analysis that all values for L 0 Λ are compatible around α = 0.1, despite considerable deviations at larger couplings.This means we can confidently extract L 0 Λ in this regime.Our final value is obtained from fit C, taking the n = 4 estimate at ν = 0, viz which is slightly more precise than the result quoted in [1], due to a refined model for the O(a) boundary effects, cf.Subsect.3.7.For an even more conservative error estimate one could take fit D, again at n = 4 and ν = 0, which yields L 0 Λ = 0.0303 (11).
Using the fits E, F and H, in terms of the fitted β-function, the values in table 7 are obtained.The fact that these are all compatible, with very similar central values further boosts the confidence that our final result is very robust.Finally, coming back to the question raised in Section 2 about exponentially suppressed contributions, we emphasize that the consistency of our analysis with fits taking the same functional form as higher order perturbative terms provides indirect evidence for the absence of such non-standard terms within our numerical precision.

Alternative tests
So far, our strategy has been to first determine Λ in the SF scheme and then convert it to Λ MS .However, one might also match the SF to the MS-coupling at 2-loop order using Eq.(2.41) and then extract the Λ-parameter within the MS-scheme.While the perturbative precision is parametrically the same as before, we present this alternative view here, as it is closer to the strategy often used in phenomenological applications.
Within the MS-scheme we have, with µ = s/L, and where s is an additional scale parameter and p ν i (s) = c ν i (s)/(4π) i , cf.Eq. (2.41).The unknown 3-loop and higher order terms in the argument of ϕ MS will be neglected in the following.The function ϕ MS , Eq. (2.47) can be evaluated using up to 5-loop order for the β-function.For our range of α-values, the numerical difference between 4-or 5-loop order evaluation is found to be negligible.The dominant uncertainty is due to the 2-loop truncation of the perturbative conversion to the MS coupling, which multiplies the sensitivity to a change in the coupling, and thus induces an O(g 4 ) or O(α 2 ) uncertainty in the estimate of the Λ-parameter.As mentioned before, this is parametrically the same as previously, cf.Eq. (4.16).We now use the non-perturbative results for the SF ν -couplings from table 6 as input in Eq. (4.28) and study the dependence of the Λ-parameter estimates on the choice of scale L n , the scale factor s and the parameter ν. Figure 6 shows some typical results; at fixed ν and s we observe again an approximate linearity in α 2 , with the asymptotic values being compatible with our best estimate, Eq. (4.27).However, we note that the slope varies significantly as a function of s and ν.
We find that the choice of s = s , Eq. (2.44), which eliminates the one-loop term in the matching, Eq. (2.45), is often (but not always) a good one.For the cases ν = 0 and ν = −0.5, figure 5 shows the 1-and 2-loop matching coefficients to the MS-coupling, Eq. (2.41), as functions of the scale factor s. The values for s are roughly around 3, 5 and 2 for ν = 0, −0.5 and 0.3, respectively.While for ν = 0 (similarly for ν = 0.3) the two-loop coefficient is near minimal around s and stays positive (figure 5, left panel), a more complicated pattern is seen for ν = −0.5 (figure 5, right panel).
A common method to assign a systematic error to a perturbative uncertainty consists in a renormalization scale variation by a factor 2 in both directions, around some "optimal scale" (cf. the review of QCD in [2]).Taking the values s as our optimal values for the scale factor we can now assess how this method fares in our context.In figure 7 this systematic error is shown, together with the total errors, for the estimates of the Λparameter.As one might expect, this systematic error dominates the error at low energies, reduces proportional to α 2 and becomes negligible at higher energies.This is indeed the case for ν = 0 and ν = 0.3, where the systematic errors are seen to bracket the shaded area representing our reference result (4.27).However, this is not the case for ν = −0.5,where a significant underestimation of the systematic error is observed, cf.figure 7.
However, in all cases we note that for α ∼ 0.1, the systematic uncertainty of the matching with perturbation theory, obtained by varying s is well below the statistical uncertainties.Moreover the latter are in line with the errors obtained with our previous strategy.This further reinforces our previous conclusions: thanks to the high energies reached with the step scaling method, our uncertainties are fully dominated by statistical errors, systematic uncertainties being negligible.The spread of results obtained by the variation of the perturbative matching scale provides a way to assess the systematic uncertainties which works well with the SF schemes at ν = 0, 0.3, even at α ≈ 0.2 (although the Final result ν = −0.5 ν = 0 ν = 0.3 systematic uncertainty there is large).But the failure of this method for the SF scheme with ν = −0.5 indicates that this method may not always be reliable, particularly if the coupling is not small and cannot be varied.This illustrates that perturbative truncation errors are very difficult to estimate within perturbation theory, and that reaching high energies is crucial for a robust determination of the strong coupling.Indeed we see that for values α ≈ 0.1 there is nice agreement between all schemes and reasonable choices of the scale factor s, within errors which clearly allow us to meet the target accuracy of 3 per cent for the Λ-parameter.

Conclusions and Outlook
Using numerical simulations and finite volume step-scaling techniques, we have studied a family of SF couplings, parameterized by ν, over a range of scales corresponding to energies of 4-128 GeV, thus differing by a scale factor 32. This, together with an unprecedented control of statistical and systematic errors represents a luxury which we have exploited to test the accuracy of perturbation theory.Choosing the Λ-parameter for ν = 0 in units of L 0 ≈ 1/(4 GeV) as a reference, its evaluation requires the knowledge of a coupling and its β-function between a finite and the infinite energy scale, where the coupling vanishes by asymptotic freedom.Perturbation theory to 3-loop order is available for the asymptotic scale dependence beyond an energy scale 1/L, which can be chosen anywhere in the range covered by our non-perturbative data, provided the ratio L/L 0 is known.By looking at the spread of values for L 0 Λ one therefore tests the accuracy of perturbation theory at the scale 1/L.Moreover, the exact relation between Λ-parameters of different schemes requires a one-loop matching of couplings which is known in all cases considered.It is therefore also possible to test the robustness of the Λ-parameter determination by using SF-schemes at various values of ν as an intermediate step.The result is neatly illustrated in Figure 4, where all data points should coincide up to a parametric uncertainty of order α 2 .We conclude that a target precision of better than 3% for L 0 Λ (which approximately corresponds to a 0.5% precision for α s (m Z )) requires non-perturbative data for a large enough range of couplings so that the perturbative truncation errors can be safely estimated.Our range of scales 4 − 128 GeV reaching down to α ≈ 0.1 allows us to reach such a precision.While some schemes may give compatible results even at α ≈ 0.2, it seems all but impossible to anticipate the quality of a given scheme if the coupling cannot be varied significantly.
With the hindsight of our 2.3% precision result for L 0 Λ, Eq. (4.27), we have also looked at an alternative test, which is close to procedures widely used in phenomenology.Namely, we have converted our non-perturbative observable, an SF ν -coupling with some choice for ν and L, to the MS-coupling where we allowed for a relative scale factor s in this perturbative conversion.Given the coupling in the MS-scheme the full machinery with up to 5-loops for the β-function [3][4][5][6][7] is available to extract the Λ-parameter.However, as in phenomenological applications, the limiting factor is the perturbative order in the conversion to the MS-scheme.We can perform this step at 2-loop order; for comparison we note that the 5-loop, O(α 4 ) result for the R-ratio [46] translates to 3-loop order when formulated as a conversion between couplings.Looking at the dependence of the scale factor, a common method consists in identifying an "optimal scale factor", s , and then vary this factor between s /2 and 2s to obtain a systematic error estimate (c.f. the review of QCD in [2]).It is a bit of an art to determine the "optimal scale factor", and some appeal to the kinematics or the physics of a given observable is often made in this context [2].We here applied such a procedure, choosing s close to the ratio of Λ-parameters, which means the one-loop coefficient in the conversion to the MS scheme is made very small.As illustrated in Figure 7, this procedure gives an error that shrinks proportionally to α 2 and often brackets the correct result.However, we have also found cases (e.g.ν = −0.5)where this procedure does not work and underestimates the systematic effect substantially, even at couplings around α = 0.15.We interpret this result as a warning: estimating errors within perturbation theory is notoriously difficult, and one may chance one's luck by being too aggressive in this step.
The work presented in this paper constitutes a major step in the α s -determination by the ALPHA-collaboration [12].Despite considerable improvements in the precision, this step currently still contributes the largest single error in this project.One may therefore hope for further progress, perhaps by combining the SF ν schemes with alternative schemes.Gradient flow couplings are obvious candidates, provided the problems with large cutoff effects can be solved [13,47], and the perturbative information is pushed at least to the same level as for the SF coupling.The latter step is possible based on numerical stochastic perturbation theory [31][32][33].Finally we note that, given the coupling results, similar nonperturbative tests of perturbation theory might also be performed using the quark mass parameters [48].

12 Figure 2 :
Figure 2: The step-scaling function for the ν = 0 SF-coupling.The band shows our result (fit C, cf.table5).The data points are the approximations at finite L/a = 6, 8, 12 taken from table 4 with errors from Eq. (3.32).

Figure 6 :
Figure 6: Determination of L 0 Λ MS at different physical scales (parametrized by the value of α in the x-axes), and using different renormalization scales (value of s) to match with the MS scheme.The left (right) panel uses the SF ν -scheme with ν = 0 (ν = −0.5),cf.text.

Figure 7 :
Figure 7: Statistical (interior error band) and total (exterior error band) uncertainties in the determination of L 0 Λ MS .The total uncertainty is computed by adding in quadratures the statistical and systematic uncertainties.The latter are computed varying the renormalization scale by a factor 2 above and below the value s .See text for more details.

Table 1 :
Coefficients in the asymptotic expansion of the β-function in different schemes.Note that the universal coefficients for N f = 3 are (4π)b 0 ≈ 0.716197, (4π) 2 b 1 ≈ 0.405285.

Table 3 :
Simulation parameters and results on the L-lattices.The hopping parameter κ was tuned such that the PCAC mass m(L), Eq. (3.15), vanishes.

Table 4 :
Simulation parameters and results on the doubled lattices.The hopping parameter κ was tuned such that the PCAC mass m(L/2) vanishes, cf.Eq.(3.15).
6= 64.The results for u n are collected in table 6, for the various fit functions representing σ(u). .6.

Table 7 :
L 0 Λ obtained with the fits to the coefficient b fit 4.5 Continuum extrapolation of Ω and Ω