Determination of the QCD $\Lambda$-parameter and the accuracy of perturbation theory at high energies

We discuss the determination of the strong coupling $\alpha_\mathrm{\overline{MS}}^{}(m_\mathrm{Z})$ or equivalently the QCD $\Lambda$-parameter. Its determination requires the use of perturbation theory in $\alpha_s(\mu)$ in some scheme, $s$, and at some energy scale $\mu$. The higher the scale $\mu$ the more accurate perturbation theory becomes, owing to asymptotic freedom. As one step in our computation of the $\Lambda$-parameter in three-flavor QCD, we perform lattice computations in a scheme which allows us to non-perturbatively reach very high energies, corresponding to $\alpha_s = 0.1$ and below. We find that (continuum) perturbation theory is very accurate there, yielding a three percent error in the $\Lambda$-parameter, while data around $\alpha_s \approx 0.2$ is clearly insufficient to quote such a precision. It is important to realize that these findings are expected to be generic, as our scheme has advantageous properties regarding the applicability of perturbation theory.


INTRODUCTION
The fundamental parameter of the strong interactions, the coupling α MS (µ) =ḡ 2 MS (µ)/(4π), is an essential input parameter for theory predictions of high energy processes in particular the physics at the LHC [1][2][3]. Conventionally the running α MS (µ) is quoted at the electroweak scale, µ = m Z . There the coupling is weak, α = O(1/10), perturbation theory (PT) is usually accurate. In particular α MS (m Z ) is essentially equivalent to the renormalization group invariant Λ-parameter because the function is known precisely by replacing the renormalization group β-function by its perturbative expansion β pert s (g) = −g 3 l b −1 n=0 b n,s g 2n ; in the MS-scheme β pert MS (g) is known up to l b = 4 loops [4,5].
At lower energies, µ m Z , the perturbative uncertainty in approximating β s ≈ β pert s in eq. (2) is generally not negligible. It is ∆Λ s /Λ s = ∆ϕ s /ϕ s = c l b α l b −1 + . . . with coefficients c l b , which are, for l b ≤ 4, of order one in the MS scheme and expected to be so in "good" schemes in general.
While the MS scheme makes sense only perturbatively, physical schemes defined beyond the perturbative expansion are easily derived from short-distance QCD observables O s (µ) = c s It is clear that high energies µ (small α s ) and at least l b = 3 are needed if one aims for a precision determination of α MS (m Z ). Replacing high energy by just a larger l b is dangerous because the perturbative expansion is only asymptotic, not convergent, and non-perturbative "corrections" can be large. In particular, whether one has lost control is difficult to detect because our knowledge of non-perturbative physics is very incomplete. Thus it is a challenge to reach an accuracy of a few percent in Λ MS equivalent to sub-percent accuracy in α MS (m Z ). Unfortunately, the determinations which quote the smallest uncertainties do typically not come from observables at large µ and uncertainties are dominated by systematics such as unknown higher order perturbative and non-perturbative terms. Both the Particle Data Group [6] and the Flavour Lattice Averaging Group [7] are therefore not just taking weighted averages of the individual determinations to arrive at their world averages.
Here we consider a family of observables (schemes) where lattice simulations allow one simultaneously to reach high precision and high energy before using PT. Then PT at µ = O(m Z ) can be employed with confidence. In addition one can check its applicability at lower scales.

arXiv:1604.06193v2 [hep-ph] 27 Sep 2016
The crucial feature enabling the study of PT at high energy with continuum extrapolated non-perturbative lattice results is that we use a finite volume renormalization scheme [8,9]. QCD is considered inside a small volume of linear extent L with boundary conditions and observables which do not contain any other scale. Details will be presented below. The renormalization scale then is and the continuum limit of lattice simulation results is easily taken for L/a 1, with modestly sized lattices. This is the strategy of the ALPHA collaboration but so far it was mostly restricted to unphysical models with an insufficient number of quark flavors [9][10][11]. For the interesting case of N f = 3 QCD, the strategy was applied by the CP-PACS collaboration [12]. We now have very precise results for N f = 3 which allow us to see important details previously hidden by uncertainties (see also [13]).
In this letter we discuss the most essential step: the accuracy of PT for couplings α 0.2 and our resulting precision for Λ. We will see that it is crucial to nonperturbatively reach α ≈ 0.1 to have confidence in PT at the 3-4 percent level in Λ. On the other hand at α ≥ 0.15 and using the three-loop beta-function, one of our schemes (ν = −1/2) shows a 10% systematic error in Λ. This is not a statistical fluctuation as we will demonstrate by eq. (21).
Given that a priori our scheme has favorable properties for PT and that other tests of perturbation theory with similar precision and similarly small α are presently not available, our result gives reason for concern in determinations of α MS (m Z ) from µ values of a few GeV. This kind of lack of accuracy of PT may be one of the sources of the spread of results reviewed in [6].

THE SF SCHEME
Our scheme is based on the so-called Schrödinger functional (SF) [14]. There are several introductory texts on the topic with emphasis on different aspects, from the general field theoretic concept [15] to detailed descriptions [16,17] and a review of concepts, history and recent results [18]. Here we just summarise those aspects which are needed to judge our findings below. Dirichlet boundary conditions are imposed in Euclidean time, for k = 1, 2, 3. The gauge potentials A µ are taken periodic in space with period L. 1 The six dimensionless just depend on the two real parameters η, ν, which multiply the Abelian generators of SU(3). With these boundary conditions the field which minimizes the action is unique up to gauge equivalence [9] and denoted by A µ = B class µ . In the temporal gauge, x 0 /L and corresponds to a constant color electric field. A family of couplings [21],ḡ ν , is then obtained by taking 1/O ν in eq. (3) to be the η-derivative of the effective action. This yields a simple path integral expectation value, which is well suited for a Monte Carlo evaluation in the latticised theory. Small fluctuations around the background field generate the non-trivial orders in PT. It is worth pointing out that the whole one-parameter family of couplings can be obtained from numerical simulations at ν = 0, as the ν-dependence is analytically known, in terms of the ν = 0 observablesḡ 2 ≡ḡ 2 ν=0 andv. Advantageous properties of these couplings are: ν +O(ḡ 6 ν ), for ∆ stat the statistical error at a given length of the Monte Carlo sample. This property makes it possible to maintain high precision at high energy. 2. The typical ∼ µ −1 , µ −2 renormalon contributions [22] are absent since the finite volume provides an infrared momentum cutoff. Instead, the leading known non-perturbative contribution is due to a secondary stationary point of the action [23] at g 2 0 [S(B sec ) − S(B class )] = 32.9. It generates corrections to PT of order which evaluates to O(10 −6 ) for α = 0.2. At such values of α, fields with non-zero topology are even further suppressed given that 3. The β-function is known including its three-loop term, and for reasonable values of ν the three-loop term is of order one as it is in the MS scheme.
4. As we will see discretisation effects are very small; at tree-level of perturbation theory they are O((a/L) 4 ). They are known to two-loop order in PT [24] and we can subtract those pieces [25].
The downside of the SF scheme is that the coefficient s(a/L) diverges like (L/a) 1/2+z for large L/a and is not that small in general. Here z is the dynamical critical exponent of the algorithm while the 1/2 in the exponent is due to the variance of the observable [25]. High statistics is needed and our computation is limited to L/a ≤ 24. A second issue is the acceleration of the approach to the continuum limit through Symanzik improvement. With our Dirichlet boundary conditions the Symanzik effective Lagrangian contains terms located at the time-boundaries. These are responsible for O(a) effects. We cancel them by corresponding improvement terms with coefficients c t andc t known only in PT, see below.

STEP SCALING FUNCTIONS AND Λ-PARAMETER
The non-perturbative energy dependence of finite volume couplings is constructed from the step scaling function [8] where m = 0 ensures the quark mass independence of the scheme [26]. The step scaling function corresponds to a discrete version of the β-function and is computed as the continuum limit a/L → 0 of its lattice approximants Σ ν (u, a/L). The conditionsḡ 2 ν (1/L) = u and m = 0 then refer to a (L/a) 4 lattice and fix the bare coupling and bare quark mass of the theory.ḡ 2 ν (1/(2L)) is to be evaluated for the same bare parameters on a (2L/a) 4 lattice.
We will use the ν = 0 scheme for a reference, dropping the index ν. The scale L 0 is defined by a value u 0 and the conditionḡ The solution of the implicit equation for u k+1 , k = 0, 1, . . . gives a series of couplings u k = g 2 (2 k /L 0 ). With a few steps one reaches µ = 1/L n = 2 n /L 0 = O(m Z ) and the perturbative ϕ at this high scale will give a good approximation to L 0 Λ Note that thanks to eq. (7) and the exact relation between Λ-parameters [9,20] r ν = Λ/Λ ν = e −ν×1.25516 , the same combination L 0 Λ can be obtained in any scheme with ν = 0. Whether different values of ν, number of steps (n) and perturbative orders (l b ) give consistent results is an excellent way to test the reliability of perturbation theory.

SIMULATIONS
We used the standard Wilson plaquette action and three massless O(a)-improved [27,28] quarks simulated by a variant of the openQCD code [29,30]. At eight couplingsḡ 2 (1/L) in the range 1.11−2.02, we simulated pairs of lattices L/a, 2L/a with L/a = 4, 6, 8 and at three couplings we also included L/a = 12.
Between 80k and 300k independent Monte Carlo measurements were taken on each lattice. As we have already noted, non-trivial topology is very suppressed in these small volumes [31]. Therefore topology freezing [32,33] is irrelevant here.
A critical issue for any lattice computation is the removal of discretization effects. In preparation of our continuum extrapolations we apply both Symanzik improvement of the action and perturbative improvement of the step scaling function [25]. In comparison to earlier work, we here propagate the estimated uncertainty of those O(a) improvement coefficients which are only known perturbatively into the errors of the step scaling functions. They can then be assumed to be free of linear a-effects within their errors. Details are found in the supplementary material attached at the end of the paper.

CONTINUUM EXTRAPOLATIONS AND RESULTS
As the residual linear a effects are treated as an uncertainty, we can proceed with continuum extrapolations linear in a 2 . First we look at the data in figure 1. They are statistically compatible with having no a-effects for L/a ≥ 6; for N f = 0 this was found with similar precision for L/a ≥ 5 (see Fig. 3 of [34]).
Both the continuum limit of the step scaling function and its cutoff effects are smooth functions of the coupling. This motivates global fits of the form where i is the order of PT to which cutoff effects are removed in eq. (23). We performed various such fits in order to assess the systematic errors which result from the assumptions made in the fit functions. We parameterize the cutoff effects by a polynomial in u, with the correct asymptotics for small u, where the case of neglecting cutoff effects is covered by n (i) ρ = 0. The continuum step scaling function is naturally parameterized by a polynomial in u, Lower order coefficients are fixed to their known perturbative values while s 3 ("n c = 1") or s 2 , s 3 ("n c = 2") are fit parameters. A selection of such fits are illustrated in table I. Instead of the parameters of the continuum step scaling function the table shows directly the extracted L 0 Λ, where L 0 is defined through eq. (11) and the value u 0 = 2.012. Recalling eq. (7) and usingv = 0.1199(10) (see next section) we havē Apart from the form of the fit, L 0 Λ depends on the value of n where eq. (13) with β ν = β pert ν is used. Since we insert β pert ν at three-loop, the residual dependence on the coupling is O(α 2 (1/L n )).
The observed behavior, figure 2, is consistent with a dominatingly linear dependence of L 0 Λ on α 2 (1/L n ). For ν = 0 the slope is not very significant and for ν = 0.3 it about disappears, but for ν = −0.5 it is quite large and outside errors.
This suggests to perform alternative fits, where the continuum step scaling function is parameterized by an effective four-loop β-function, adding a term b eff 3 g 9 to the known ones. The determined L 0 Λ are then automatically independent of n and we include b eff 3 instead of u n=4 in the table. For ν = −0.5 the effective fit value is larger fit  than it should be in a well-behaved perturbative expansion. We will come back to this issue shortly, but first we give our result for L 0 Λ. We take the standard polynomial fit to σ (for ν = 0) with α n ≈ 0.1 (u n ≈ 1.2). A typical perturbative error of size ∆(ΛL n ) = α 2 n ΛL n is then a factor 3 or more below our statistical errors. We quote (withḡ 2 (1/L 0 ) = 2.012) (21) , (19) with the known Λ MS /Λ [9,20]. This is the result of fit C. It is in perfect agreement with all variations of the global fit, even with fit G, which neglects all cutoff effects but uses only data with L/a ≥ 8. It has a rather conservative error. If an even more conservative result is desired, one may take the one of fit D, L 0 Λ = 0.0303 (13).

ACCURACY OF PERTURBATION THEORY
While b eff 3,ν is large for ν = −0.5, it does have an error of around 50%. A much better precision can be achieved by considering directly the observable with coefficients v 1 = 0.14307, v 2 = −0.004693 [11]. In contrast to the step scaling function ω(u) does not require pairs of lattices, so that the continuum extrapolation can be performed using data for the entire range of lattice sizes L/a = 6, 8, 10, 12, 16, 24. Improvement and fits for obtaining the continuum limit are carried out in analogy to those of Σ ν . Figure 3 shows the result of two different fits with fit parameters d k in The overall band of the two fits may be taken as a safe estimate of the continuum limit. As an example we find ω(2.012) = 0.1199(10) for both fits, leading to eq. (18). In the above analysis we did not use data with L/a = 6. Including them yields only tiny changes and excellent χ 2 values.
A good measure of the deviation from two-loop perturbation theory is at α = 0.19. It is quite large and statistically significant beyond any doubt. If one attempts to describe this by perturbation theory, the three-loop coefficient v 3 has to be too large for perturbation theory to be trustworthy at α = 0.2. Again, we come to the conclusion that α ≈ 0.1 needs to be reached non-perturbatively before perturbation theory becomes accurate.

SUMMARY AND CONCLUSIONS
Our chosen definition of α s (µ) allows us to compute it with very good precision through lattice Monte Carlo simulations. In particular we have controlled the errors due to the discretisation of the theory also at large µ. Known non-perturbative corrections are parametrically very small: O(e −2.6/α ). In other words we have an excellent scheme to test the accuracy of PT in a given region of α.
In fact, we have a family of schemes, depending on ν. For small positive ν, the couplings follow perturbation theory very closely in the full investigated range 0.1 ≤ α ≤ 0.2 as illustrated by the flatness of Λ in figure 2 extracted from eq. (13) with the three-loop β-function.
However, for negative ν, e.g. ν = −0.5, values of α just below 0.2 are not small enough to confirm perturbative behaviour. The observablev, figure 3, shows that the α-dependence seen in figure 2 is not just a statistical fluctuation. We could take the continuum limit ofv with very high precision and eq. (21) shows a clear deviation from the known perturbative terms, corresponding to l b = 3, at α = 0.19.
We conclude that it is essential to reach α = 0.1 in order to be able to achieve a precision around 3% for the Λ-parameter. Fortunately we have access to that region and can quote such an accuracy in eq. (19). While of course schemes exist where three-loop running holds accurately down to smaller energies -for example ν = 0.3 produces flatness in figure 2 as far as we can tell -to know whether a chosen scheme possesses this property is difficult unless one has control also over the α ≈ 0.1 region. Once that is achieved larger α are not much needed any more.
What we reported in this letter is part of our determination of a precise value for Λ MS . As our next step, we will soon connect L 0 to the decay constants of pion and kaon, as explained above and in [35].
We thank our colleagues in the ALPHA collaboration, in particular M. Bruno, C. Pena, S. Schaefer, H. Simma and U. Wolff for many useful discussions. We thank U. Wolff for a critical reading of the manuscript. We would also like to show our gratitude to S. Schaefer and H. Simma for their invaluable contribution in adapting the openQCD code. We thank the computer centres at HLRN (bep00040) and NIC at DESY, Zeuthen for providing computing resources and support. S.S. acknowledges support by SFI under grant 11/RFP/PHY3218. P.F. acknowledges financial support from the Spanish MINECO's "Centro de Excelencia Severo Ochoa" Programme under grant SEV-2012-0249, as well as from the MCINN grant FPA2012-31686. This work is based on previous work [18] supported strongly by the Deutsche Forschungsgemeinschaft in the SFB/TR 09.

SUPPLEMENTARY MATERIAL
Here we add some details on the Symanzik improvement of the action and perturbative improvement of the step scaling function [25]. In particular we discuss how the uncertainties in improvement coefficients are handled.

Improvement of the action
Apart from the bulk O(a) improvement term, the complete removal of linear (in a) discretization errors requires a boundary improvement coefficient c t in the gluon action [14] and a coefficientc t in the fermion action [27]. Regarding the former, the known two-loop accuracy [24], c t = 1 + g 2 0 (−0.0890 + 0.019141N f ) + g 4 0 (−0.0294 + 0.002N f + 0.000(1)N 2 f ) , (22) is expected to be sufficient [34] since we are in the weak bare coupling region. Still, we propagate the small deficit of O(a)-improvement into our errors. As an uncertainty in c t we use the full two-loop term, adding ∆ ct Σ = 0.0234g 4 0 |∂ ct Σ| in quadrature to the statistical errors. The derivative ∂ ct Σ was obtained from a numerical estimate of δ b (ḡ 2 ) ≡ L Similarly, for the coefficientc t (g 0 ) = 1 − 0.01795g 2 0 + O(g 4 0 ) [36] we use the full one-loop term as an error estimate. Its effect is much smaller than the one of the uncertainty in c t , since it contributes only through quark loops.
These errors are responsible for around 30% of the uncertainties of some of our central results in TABLE I. Throughout error propagation is carried out as detailed in [13].

Perturbative improvement of the step scaling functions
In addition, we can improve the observables to a given order i in perturbation theory but to all orders in a via Σ (i) (u, a/L) = Σ(u, a/L) with δ 1 , δ 2 known [24].