{'fraction_non_alphanumeric': 0.03682532526085276, 'fraction_numerical': 0.01637575679505346, 'mean_word_length': 4.155653328905861, 'pattern_counts': {'":': 0, '<': 24, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Determination of alpha_s from the QCD static energy: an update

We present an update of our determination of the strong coupling alpha_s from the quantum chromodynamics static energy. This updated analysis includes new lattice data, at smaller lattice spacings and reaching shorter distances, the use of better suited perturbative expressions to compare with data in a wider distance range, and a comprehensive and detailed estimate of the error sources that contribute to the uncertainty of the final result. Our updated value for alpha_s at the Z-mass scale, M_Z, is alpha_s(M_Z)=0.1166^{+0.0012}_{-0.0008}, which supersedes our previous result.


I. INTRODUCTION
The strong coupling, α s , is the only free parameter of quantum chromodynamics (QCD) in the massless quark limit. The knowledge of its value with good precision and accuracy is of vital importance for the study of any process that involves the strong interactions [1]. It has, for instance, a big impact on cross-section calculations for the Large Hadron Collider, where its uncertainty is often one of the limiting factors on the precision reached by theoretical predictions. The present world average by the Particle Data Group (PDG) [2] quotes an uncertainty of about ±0.5% for α s at the scale of the Z-boson mass, M Z . But, when one looks at the individual values that enter in the average, one realizes that the uncertainties of many of the individual measurements are dominated by errors of theoretical origin. These errors are often difficult to precisely assess, and in fact some of the individual measurements are not in good agreement with each other. A flagrant example of that is given, for instance, by comparing the lattice result obtained in Ref. [3] with the value obtained from analyses of the event-shape variable thrust in e + − e − collisions of Ref. [4]. Disagreements like that question, to some extent, if the uncertainty quoted in the world average really reflects our current understanding of the value of α s . In that sense it is crucial to have an increasing corroboration of the value of α s , by extracting it from different independent quantities, and at different energy ranges, and to critically and exhaustively analyze all the theoretical errors that enter in each of the determinations. In Ref. [5] we presented a novel determination of α s , based on the comparison of lattice data for the energy between two static sources in the fundamental representation of QCD, i.e. the QCD static energy, at * Current address: Department of Physics and Astronomy, University of Iowa, 203 Van Allen Hall, Iowa City, Iowa 52242-1479 short distances with perturbative expressions. This extraction was possible due to recent progress both in the lattice evaluation and the perturbative computation of the static energy. It provided a competitive determination of α s , that stemmed from a perturbative calculation at three-loop order, including resummations of logarithmically enhanced terms. In this paper we present an update of the result of Ref. [5]. The main new ingredients are: (i) we include new lattice data, with smaller lattice spacings, and reaching shorter distances, (ii) we use perturbative expressions that are better suited to compare with lattice data in this wider distance range, and (iii) we include a more detailed assessment of the different errors sources that can contribute to the uncertainty of the final result; the discussion of the different errors sources comprises, in fact, the bulk of the present paper. The rest of the paper is organized as follows. In Sec. II we describe the lattice data we use, in particular explaining in detail the lattice systematic errors at short distances. Section III displays the perturbative expressions we employ. The bulk of the paper is contained in Sec. IV, where we explain in detail the analysis to extract α s , and discuss all the error sources that can affect our result. Section V contains some discussion of our results, as well as comparison with previous related works. Finally, in Sec. VI we present a short summary of the main results and conclude.

II. LATTICE DATA
In this paper we use 2+1-flavor lattice QCD data on the static quark-antiquark energy, E 0 (r), obtained using tree-level improved gauge action and Highly-Improved Staggered Quark (HISQ) action, procured by the HotQCD collaboration [6]. The strange-quark mass m s was fixed to its physical value, while the light-quark masses were chosen to be m l = m s /20. These correspond to a pion mass of about 160 MeV in the continuum limit, which is very close to the physical value. More precisely, we use lattice QCD data corresponding to the lattice gauge couplings β = 10/g 2 = 7.150, 7.280, 7.373, 7.596 and 7.825. The largest gauge coupling, β = 7.825, corresponds to lattice spacings of a = 0.041 fm. One may worry about the evolution of the topological charge on such fine lattices, but, as it was demonstrated in Ref. [6], the Monte-Carlo evolution of the topological charge is acceptable even for β = 7.825. The lattice spacing was fixed using the r 1 scale defined as The values of r 1 /a for each of the above gauge couplings have been determined in Ref. [6] and are summarized in Tab. I. We are mostly interested in the static energy at short distances, where lattice artifacts can be significant. Therefore, it is important to remove these artifacts and to estimate the corresponding systematic errors before comparing the lattice results with perturbation theory. Since we are interested in the cutoff effects at short distances, we could use lattice perturbation theory to estimate these effects, due to asymptotic freedom in QCD.
In the simplest case we could use the tree-level result. At tree level the static energy on the lattice is given by where D 00 (k 0 = 0, k) is the temporal gluon propagator on the lattice, which depends on the choice of the gauge action, and C F is the Casimir of the fundamental representation. The above equation defines the so-called improved distances r I (r) (for simplicity, from now on we omit the dependence on r and just write r I ), that renders the lattice Coulomb potential to be continuum-like and will be used below to reduce lattice artifacts. In Fig. 1 we show the ratio of the lattice static energy to the continuum static energy at tree level. As one can see, the discretization effects are small, at the level of few percent, even at distance r/a = 1. For the improved gauge action the discretization effects are below 1% level for r/a > 2, while they are still significant for the standard gauge action. To estimate the cutoff effects in the actual lattice calculations we need a continuum estimate of the static energy. Based on the free-theory result we assume that the cutoff effects are negligible for distances r/a > 2, and fit the lattice results at β = 7.825 in this region to Coulomb plus linear plus constant form, to obtain a continuum estimate. To reduce the cutoff effects we replace the distance r by r I in the fit (see below). The fit can be used as a continuum estimate for distances r > 0.3r 1  described in Ref. [5]. In Fig. 2 we show the ratio of the static energy at different lattice spacings (β) to the continuum estimate described above for r > 0.3r 1 . In terms of r in lattice units we cover the distance range r/a > √ 2. Obviously, the range in r/a extends to smaller values for smaller values of β. As one can see, the discretization effects are similar to the free theory, although slightly larger in some cases. Thus, the tree-level result gives us a fair estimate on the size of discretization effects, as expected. We can use the free-theory result to correct for the lattice artifacts at short distances, i.e. implement a tree-level improvement of the static energy. One way to correct for the lattice artifacts is to replace the distance r by the improved distance r I defined in Eq. (2). We will refer to this improvement scheme as improvement 1. Alternatively, we can fit the static energy to the form and correct the lattice data by subtracting a (1/r−1/r I ).
We call this improvement scheme improvement 2. We see from Fig. 2 that the tree-level corrections reduce the discretization errors to the level of 1.5% or smaller. Furthermore, the deviations from the continuum estimate are roughly independent of β in the considered range of the gauge couplings. From the figure, one can also see that the two improvement schemes give similar results. We adopt improvement scheme 1 and try to correct for the residual cutoff effects by dividing the value of the static energy calculated on the lattice with the correction factors given in Table II. The correction factors have been estimated using an iterative procedure as follows.
First we calculate the ratios of the static energy to the above continuum estimate for r/a = √ 2 and β = 7.280, r/a = √ 3 and β = 7.373, and r/a = 2 and β = 7.596, and divide the lattice results on the static energy for β = 7.825 for r/a = √ 2, √ 3 and 2 by the corresponding ratios. Furthermore, we assign a systematic error of 0.5%  to these data points and fit the corrected data points for β = 7.825 by a Coulomb plus linear plus constant form to obtain the improved continuum estimate of the static energy, which now extends to smaller distances, namely down to distance r = 0.186r 1 . With this continuum estimate we calculate the correction factor for the smallest distance r/a = 1 by dividing the β = 7.373 data at r/a = 1 by the corresponding value. Using this correction factor for the β = 7.825 static energy at r/a = 1, and assigning a systematic error of 1% to that point, we could extend the fit to the smallest lattice distance to obtain our final continuum estimate for the static energy for distances 0.12 < r/r 1 < 0.6, and also the correction factors given in Table II. We assign a systematic error of 1% for all the corrected lattice data at distance r/a = 1 and a systematic error of 0.5% for all the remaining 5 shortest distances. The corrected lattice data for the static energy with the statistical and systematic errors added in quadrature is shown in Fig. 3. The figure clearly shows that, within estimated errors, all the lattice data on the static energy agree with the continuum estimate, i.e. our final error estimates are sufficiently large to included all possible residual discretization effects.
We have also calculated the force from the lattice data at the three smallest lattice spacings. The calculation of the derivative of the static energy from a discrete set of lattice data is non-trivial, since both discretization errors and statistical fluctuations result in a non-smooth behav-β r/a = 1 r/a = √ 2 r/a = √ 3 r/a = 2 r/a = √ 5 r/a = ior of the static energy as function of the distance r, and subsequently in non-monotonic behavior of the force. To avoid problems with the lattice artifacts, we calculate the force only for distances r/a > 2, for which discretization effects are expected to be small. To obtain a result for the force from the lattice data on the static energy which is monotonic, we fit the latter with smoothing splines. The smoothing spline is determined by minimizing the χ 2 plus the integral of the second derivative of the fit function in the considered interval times a real parameter λ. The spline fit is performed using the R statistical package [7]. The errors on the spline are calculated using standard bootstrap method, i.e. for each data point on the static energy we generate a set of synthetic data according to a Gaussian distribution, with the width and mean given by the statistical error and central value of the lattice data, perform the spline fit on the synthetic data, and finally perform the statistical analysis of the resulting splines. Furthermore, we consider several values of the smoothing parameter λ in the standard range 0 < λ < 1. We find that, within statistical errors, there is little dependence on the smoothing parameter. For the distances where the dependence on λ was larger than the statistical errors, the final error estimate was increased to accommodate the difference between different splines. In Fig. 4  arated by a/8 in the interval 0.32 < r/r 1 < 0.7, with the corresponding statistical errors. As one can see from the figure, the values for the force obtained for different β agree well within the estimated errors, confirming our expectation that cutoff effects are small for r/a > 2. The data points for the force shown in Fig. 4 will be used in the analysis presented in Sec. IV C 1.

III. PERTURBATIVE EXPRESSIONS
In this section we detail the perturbative expressions that are used in our analyses.
Let us start by recalling that our goal is to compare with the lattice data for the static energy. An important point to remember, to that effect, is that the normalization of the static energy computed on the lattice is not physical, but only its slope is. At the same time, care needs to be taken with the so-called renormalon singularities, present in the perturbative expressions. At the end of the day, the renormalon has no practical effect in the extraction, i.e. the comparisons, and any fits to the lattice data involved in the analyses, can always be done in a way such that they are not affected by the presence of the renormalon. This is again due to the fact that the overall normalization of the energy is affected by the renormalon, but the slope, which is what matters for the extraction of α s , is renormalon free. It is nev-ertheless convenient to keep track of the presence of the renormalon in intermediate expressions, and to always work with a scheme where the renormalon is subtracted explicitly. We can schematically write the static energy as E 0 (r) = V s (r, ν, µ) + δ U S (r, ν, µ) + RS(ρ) + Λ s (ρ), (4) where V s (r, ν, µ) is the static potential, δ U S (r, ν, µ) represents the ultrasoft contributions to the energy (which first appear at three-loop order, see Eq. (9) below), RS(ρ) generically represents a term that cancels the leading renormalon singularity in V s (r, ν, µ), and Λ s (ρ) contains the residual mass term that is necessary to define the static limit of QCD. The ν dependence in V s and δ U S cancels order by order in perturbation theory. We mention that µ is the ultrasoft factorization scale, but it does not play a relevant role for the discussion in this paragraph. For concreteness, we adopt the RS scheme [8]. The expression for V s (r, ν, µ) is given as a series in α s (ν) and contains ln(rν) terms, starting at order α 2 s (ν) ln(rν). The expression for RS(ρ) is given as a series in α s (ρ). In order to cancel the renormalon order by order, α s (ρ) needs to be expanded in terms of α s (ν) (or vice versa); when one does that, RS(ρ) is given as well as a series in α s (ν) and contains ln(ρ/ν) terms, starting at order α 3 s (ν) ln(ρ/ν). If one sets ν = 1/r in Eq. (4), one resums the ln(rν) logarithms in V s (r, ν, µ), but on the other hand the RS-subtraction term is not constant. Since the normalization of the renormalon singularity, which enters in RS(ρ), can only be computed approximately, that would introduce additional uncertainties in the perturbative expression, that would hinder any precise α s extraction. One can also set ν constant in Eq. (4), in that case the RS-subtraction term is itself constant, and gets absorbed in the additive constant needed to compare with lattice data (see Eq. (8) below), but one has ln(rν) terms left in V s (r, ν, µ). If one only compares with lattice data in a restricted distance range, one is free to set ν at the center of that distance range and the logarithms never become large. This last option is the approach that was taken in Ref. [5]. Alternatively, one can also take a derivative of Eq. (4), and compute the force. We can take the derivative before resumming the ln(rν) terms, i.e. from Eq. (4) with constant ν, so that the RS-subtraction term, and Λ s , completely disappear. That is the procedure we follow here.
By taking a derivative, with respect to r, of Eq. (4) we obtain the force, where α E (r, ν) denotes the derivative of α E (r, ν) with respect to r, and where we defined At this point the expression for the force F (r, ν) is given as a series in α s (ν) and contains ln(rν) terms, starting at order α 2 s (ν) ln(rν). We can now resum these logarithms by setting ν = 1/r, i.e. we re-organize the perturbative expansion for the force as a series in α s (1/r). Having obtained this expression, one can either directly compare F (r, ν = 1/r) with lattice data for the force, or numerically integrate F (r, ν = 1/r) and compare with lattice data for the energy, that is, to numerically compute where the lower limit of integration, r * , is irrelevant because the constant contribution coming from it gets absorbed in the additive constant used to compare with lattice data; i.e. when we compare with lattice data we need to plot where r ref is the reference distance where we make the perturbative expression coincide with the lattice data, and E latt.
is the value of the static energy computed on the lattice at that distance. As it will be further explained in Sec. IV, we will use the comparison with the lattice data for the energy, through Eqs. (7)-(8), for our main analyses, and leave the direct comparison of F (r, ν = 1/r) with the lattice data for the force as a cross check of our results.
The perturbative expansion of F (r, ν = 1/r) is explicitly given by where tree-level accuracy corresponds to taking the first line of the equation, and i-loop order to taking terms up to order α i s in the square brackets. The coefficients a i in the equation above have been calculated over the years, and are currently known up to three-loop order, i.e. i = 3 [9][10][11][12][13][14][15][16][17]. The existence of the ln α s terms was noticed a long time ago [18], and the a L 3 coefficient was computed in Refs. [19,20]; the corresponding coefficient at the next order, a L 4 , is also known [21]. These ln α s terms, socalled ultrasoft logarithms, have also been resummed, to sub-leading accuracy; we present the ultrasoft-resummed expressions in Sec. III A. The values of theã i , a L i coefficients have been collected in Ref. [22], and we will not copy the expressions here. We refer to that reference, which uses exactly the same notation we employ here, for the concrete expressions; β i are the coefficients of the beta function and their expressions, as well as those of the color factors, are also collected there.
As mentioned above, we will use Eqs. (7)-(9) for our analyses. Note again that all reference to the renormalon is gone in those expressions. These are essentially the same expressions that are used in related analyses by, for instance, the ALPHA collaboration [23,24], where no explicit mention to renormalons is made at all. Keeping track of the presence of the renormalon is just a convenient way to identify which quantities have a wellbehaved perturbative expansion, and are adequate for a precise α s extraction. Recall that, in our previous analysis [5] we used the perturbative expression coming from Eq. (4) with ν constant. In that case the renormalon got absorbed in the additive constant that is added to E 0 (r) to compare with lattice, and did not affect the fits either. The advantage of the method pursued in the present analysis, with respect to the one in Ref. [5], is that ln(rν) terms are completely absent. When one compares with lattice data in a restricted r range, one can always choose ν at the middle of the range, and the logs do not become large, but since in the present analysis we reach shorter distances than in Ref. [5], it is better suited to use an expression that completely avoids all ln(rν) terms, as we do now. All these different options, in treating the perturbative expressions when comparing with lattice data, were already explained in Ref. [8]; the method we use here, through Eqs. (7)-(9), was mentioned but not really employed, though, presumably because the aim of that paper was mostly to study the presence of the renormalon, rather than a precise determination of a parameter.
We also mention that we explicitly checked that finite strange-quark mass effects can be neglected when comparing with lattice data, therefore we will safely ignore them in the present analysis. Finally we note that we use four-loop accuracy for the running of α s everywhere.

A. Resummation of the ultrasoft logarithms
The ultrasoft logarithms, i.e. the ln α s terms, in the static energy, were resummed at leading order in Ref. [25] and at sub-leading order in Ref. [26]. We also employ the ultrasoft-resummed expressions in our analyses here.
These ultrasoft terms can be conveniently resummed by solving Renormalization Group (RG) equations in the effective theory potential Non-Relativistic QCD (pN-RQCD) [27,28]. The RG equations in pNRQCD have been written at the level of the matching coefficients, i.e. V s and Λ s in Eq. (4). Therefore, to obtain the expressions we need to compare with lattice data, what we do is the following: We first resum the ultrasoft logarithms for the potential, i.e. V s , and Λ s , by solving the RG equations in pNRQCD from an initial scale µ 0 = ν, to a scale µ. We then take the derivative of the static energy with the ultrasoft-resummed expressions to obtain the force. Fi-nally, we re-organize the perturbative expansion for the force as a series in α s (1/r).
After doing that, the perturbative expression for the force is given by where we refer again to Ref. [22] for a concrete expression of the coefficient η 0 . The terms that are displayed in Eq. (10) account for all the contributions that are needed up to order α 4+n s log n α s (n ≥ 0), which is what we call next-to-next-to-next-to-leading logarithmic (N 3 LL) accuracy; accordingly, next-to-next-to-leading logarithmic (N 2 LL) accuracy includes the contributions up to order α 3+n s log n α s , and corresponds to the first, second, third, and fifth lines in Eq. (10). Note that the, schemedependent, constant K 2 that was present at N 3 LL accuracy in the analysis of Ref. [5] is absent in the scheme we use in the present analysis. Recall that µ is the ultrasoft factorization scale, which takes a natural value µ ∼ (C A α s )/(2r), and that the apparent µ dependence in Eq. (10) is of higher order in α s .

IV. ANALYSES FOR α S EXTRACTION
Our aim is to compare the lattice data with the perturbative expressions, as presented in the previous sections, with the goal of obtaining a determination of α s . The general concept of performing such an analysis was suggested long ago [29,30], but it is only the recent progress in the evaluation of the static energy, as summarized in the previous sections, that allowed for a practical realization [5]. The basic idea we follow here is very simple: the perturbative expressions depend on Λ MS , where Λ MS is the QCD scale (which fixes α s ), and we can use the comparison to find the values of r 1 Λ MS that are allowed by lattice data. Here, r 1 is the reference scale used in the lattice computation, its value has also been independently computed on the lattice, and using it we obtain Λ MS in physical units. To find the values of r 1 Λ MS allowed by lattice data, one goes to the short-distance region, where it is expected that perturbation theory is enough to describe the lattice results. The guiding principle we follow, to find the allowed α s values, is that the agreement of the theoretical predictions with lattice data should improve when the perturbative order of the computation is increased. At the same time, the error that one assigns to the result should be such that it reflects the uncertainties due to unknown higher perturbative orders.
As explained in the previous sections, to perform such an analysis, it is important to remember that the normalization of the lattice result is not physical, but only the slope is. Therefore one must either normalize the energy to a certain value at a given distance, as specified by Eq. (8), or compare directly with the force. Additionally, since we have lattice results for several different lattice spacings, see Sec. II, one can either put all the lattice data together, or perform an analysis for each lattice spacing separately and then take the average. Our previous analysis in Ref. [5] employed the procedure devised in Ref. [31] to implement the guiding principle of improved agreement with increasing perturbative orders, used all lattice data together, and compared with lattice according to Eq. (8).
One of the crucial aspects in this kind of analyses is always to know whether the current lattice data has really reached the purely perturbative regime, and with enough precision to perform the extraction, or not. It is not an easy task to undoubtedly state this point. In that sense, the facts that the analyses of Ref. [5] showed that the agreement with lattice indeed improved when increasing the perturbative order, and that the perturbative curves were able to describe the lattice data quite well, can be seen as positive evidence that this was the case. Thanks to the enlarged lattice data set that we have in the present paper, which reaches shorter distances, we can now perform a more detailed analysis of this issue.

A. Procedure to extract αs
With the main motivation of further testing if the lattice data has reached the purely perturbative regime, we modify the procedure we used in Ref. [5] for the extraction of α s , and proceed as we describe next.
First of all, we use now the data for each value of the lattice spacing separately, and at the end perform an average of the different obtained values of α s . The reason for proceeding this way is that, when one wants to put all lattice data together, one needs to normalize the results calculated at different lattice spacings to a common value at a certain distance, due to the additive ultraviolet renormalization of the static energy. There is some uncertainty related to this normalization procedure, and in fact, the errors on the lattice data due to the normalization are larger than the lattice systematic errors. On the other hand, if one uses the data for each lattice spacing separately, there is no need to normalize it, and the data errors are, therefore, smaller. The downside of this option is that, of course, one has less data in each analysis. With the enlarged data set we use in this paper, we have enough data at short distances to perform analyses for each lattice spacing separately. This was not quite the case with the data set used in Ref. [5].
The procedure used in Ref. [5] was devised to provide a faithful estimate of the central value and error of the result [31]. Here we modify it to be able to further test if we have reached the purely perturbative regime, while maintaining the reliability of the outcome. At the same time, we try to keep the whole procedure as simple and straightforward as possible. In Ref. [5] we used the distance region r < 0.75r 1 for all our analyses and fits. We now use r < 0.75r 1 as the largest distance range we reach, but also perform analyses with points at shorter distance ranges only. We proceed as follows: 1. Perform fits to the lattice data for the static energy, according to Eqs. (7)- (8), at different orders of perturbative accuracy. The parameter of the fits is Λ MS .
3. Use the ranges where the reduced χ 2 of the fits either decreases (or remains almost constant, i.e. it does not increase by more than one unit) when increasing the perturbative order, or is smaller than 1.
4. To estimate the perturbative uncertainty of the result, we repeat the fits in two ways: (i) by varying the scale in the perturbative expansion, from ν = 1/r to ν = √ 2/r and ν = 1/( √ 2r), and (ii) by adding/subtracting a term ±(C F /r 2 )α n+2 s to the expression at n loops.
As mentioned, we use the data for each lattice spacing separately. This means that, in order to have a significant number of data points for each of the distance ranges above, we only use the β = 7.373, 7.596, 7.825 sets. To illustrate how many points enter in each range, we show these data sets, with the distance ranges marked, in Fig. 5; note that this figure just aims at showing how many points we have, and, therefore, everything is displayed in units of the scale r 1 and the different data sets are normalized, but no error associated to that is shown. The actual fits are performed with the data for each set in units of the corresponding lattice spacing, a, and then the result is translated to units of r 1 using the corresponding value of r 1 /a and its error, as given in Tab. I. We always use fits that consider the lattice data for a given β to be uncorrelated. This is sufficient since we mostly consider the static energy at off-axis values of r, for which the off-diagonal elements of the correlation matrix are small. Note that the two sets of fits in point 4 of the list above both provide an estimate of the uncertainty due to unknown higher-order perturbative terms; we will take whichever of the two gives a larger uncertainty as the perturbative-error estimate. In addition, the fits using different distance ranges probe α s at different scales, and agreement among their results would be an indication that perturbative uncertainties are properly estimated.

B. Outcomes of the analyses
We present in this section the outcomes of the fits described above. We compare with lattice data for the static energy according to Eq. (8). For that we need to choose which lattice point defines r ref . In principle one would choose r ref as the shortest distance where lattice data is available, since the perturbative expressions should be more reliable there. In practice though, as discussed in detail in Sec. II, the lattice points at shorter distances have larger systematic discretization errors. Since when we normalize the perturbative curves to r ref we effectively put the error of that lattice point to zero, it would not be appropriate to use one of the points with larger systematic errors to define r ref . Therefore, we do not use any of the first six points at each lattice spacing to define r ref . Also, to make sure that the results of the fits are not sensitive to the point we choose to define r ref , one should perform analyses with different choices of r ref . According to these discussions, we choose the seventh, eighth, or ninth point at each lattice spacing to define r ref ; if for some data set at some distance range there are not enough points to do that, we choose the point preceding it in distance instead, i.e. if, for example, in a given distance range there were only six points we would choose the sixth point to define r ref .
r 0.75r 1 r 0.75r 1 r 0.7r 1 r 0.7r 1 r 0.65r 1 r 0.65r 1 r 0.6r 1 r 0.6r 1 r 0.55r 1 r 0.55r 1 r 0.5r 1 r 0.5r 1 r 0.45r 1 r 0. We perform the fits at tree-level, one-, two-, and threeloop accuracy. The reduced χ 2 of the fits, for the β = 7.373, 7.596, 7.825 sets, in each of our distance ranges are shown in Figs. 6-8. From these figures we see that in the largest distance range we considered, i.e. r < 0.75r 1 , the reduced χ 2 s do not decrease when we increase the perturbative order. When we reduce the distance range, the χ 2 curves move in the correct direction in order to satisfy the requirements in Sec. IV A. For r < 0.6r 1 and shorter distance ranges, in all cases (with the only r 0.75r 1 r 0.75r 1 r 0.7r 1 r 0.7r 1 r 0.65r 1 r 0.65r 1 r 0.6r 1 r 0.6r 1 r 0.55r 1 r 0.55r 1 r 0.5r 1 r 0.5r 1 r 0.45r 1 r 0. For illustration, we show the results we obtain for aΛ MS from the r < 0.5r 1 distance range, at different orders of perturbative accuracy, in Fig. 9. The figure is for N ref = 7, but the other cases look quite similar. The error bars in the figure are obtained by repeating the fits, with the scale variation and addition of higher-order term, as specified in the previous section. For each perturbative order, the point to the left corresponds to the ν variation and the point to the right to the addition of a generic higher-order term. (We note that, for some numerical r 0.75r 1 r 0.75r 1 r 0.7r 1 r 0.7r 1 r 0.65r 1 r 0.65r 1 r 0.6r 1 r 0.6r 1 r 0.55r 1 r 0.55r 1 r 0.5r 1 r 0.5r 1 r 0.45r 1 r 0. coincidence, the downward error bars from ν variation at one loop are very small). Both error bars reflect the perturbative uncertainty at a given order, and we take the largest of the two as our perturbative-uncertainty estimate. As we can clearly see in the figure, the central value at any given order is always nicely contained within the error bars of the previous order. This fact gives us confidence that the perturbative uncertainty is faithfully estimated. Note also the big error reduction in the determination of aΛ MS that the knowledge of high-order terms brings about. We show in Fig. 10 the results we obtain for r 1 Λ MS for the r < 0.6r 1 and shorter distance ranges, which are the ones that passed our criteria in Sec. IV A, at threeloop accuracy. We converted the results for aΛ MS for each lattice spacing to r 1 Λ MS using the corresponding value of r 1 /a and its error, see Sec. II; we added the perturbative error and the error in r 1 /a in quadrature in the error bars of Fig. 10. We note that the numbers obtained for the different distance ranges are perfectly compatible with each other. We will use the numbers for r < 0.5r 1 in our final results, since the fits in this distance range always have reduced χ 2 s around 1. The effect of changing N ref produces very small variations in the result, compare the three panels in Fig. 10 result for aΛ MS . Then we can perform a weighted average of the results at the three different lattice spacings to obtain our final number for r 1 Λ MS . Finally we convert it to physical units using the value of r 1 given in Sec. II. All these numbers and our final result for α s are collected in Sec. IV D below. In Fig. 10 we also displayed, for comparison, our previous result in Ref. [5] as the pink band. We can see that the present results are perfectly compatible with Ref. [5], and have a smaller error.

Statistical errors of the fit parameters
In addition to the perturbative errors due to unknown higher-order terms (and the error due to the uncertainty in the values of r 1 /a) that we showed above, there is also an uncertainty in the values of aΛ MS due to the fit parameter errors. We estimate these errors, which we call statistical, by taking the values of Λ MS at one χ 2 unit above the minimum. Figure 11 shows these statistical errors for the different distance ranges we consider (black -darker-error bars), together with the perturbative errors (magenta -lighter-error bars) for reference. The figure is for N ref = 7, but the other cases look quite similar. We can see that the statistical errors are always smaller than the perturbative ones. They are not completely negligible, though, as was the case in our previous analysis in Ref. [5]. We include these statistical errors in our final tables and results in Sec. IV D. We also note that the statistical errors for the r < 0.5r 1 ranges are comparable in size with the errors due to the uncertainty in r 1 /a. The final tables in Sec. IV D detail each one of the errors individually, and one can easily compare their relative sizes.

Ultrasoft resummation
As mentioned in Sec. III A, the ultrasoft, i.e. ln α s (1/r), terms have been resummed at sub-leading accuracy. The numerical importance of these terms, and the necessity or not to include their resummation, depends on the distance r we are studying. It is clear that their size grows when we go to shorter distances, since α s (1/r) → 0. In this section we want to study their numerical size in the range of distances where we compare with lattice. Note that, since in principle we do not know if it is necessary or not to include ultrasoft resummation (and at which order), it is more suitable to perform the study in the previous section (testing if the lattice data has reached or not the perturbative regime) before including the ultrasoft resummed expressions, as we have presented.
To visualize the importance of the ultrasoft terms, we plot in Fig. 12 the a L 3 ln(C A α s /2) term (solid bluedarker-line), together with the rest of the three loop term for the force, i.e.ã 3 − 6ã 2 β 0 − 4ã 1 β 1 − 2β 2 (dashed blue -darker-line), see Eq. (9). From the figure we see that the ultrasoft log is much smaller than the non-logarithmic term for r ∼ 0.75r 1 , but it grows quite rapidly as r decreases; it becomes comparable to it for r ∼ 0.2r 1 , and it is larger than the non-logarithmic part for the shortest distances where we have lattice data. We also notice that the ultrasoft logarithm and the non-logarithmic part have opposite signs, and largely compensate each other for distances around 0.2r 1 . Therefore, although the size of the ultrasoft logarithm at three loops is never much larger than the non-logarithmic part, it is certainly suitable to include leading-ultrasoft resummation in our analyses. Consequently, results including leading-ultrasoft resummation are also collected in Sec. IV D below. One would not expect, though, a big difference between the results at three loop with and without ultrasoft resummation.
In principle one could also argue that the leading ultrasoft resummation should be included together with the two-loop term, which would correspond to counting ln α s ∼ 1/α s . Figure 12 also shows the two-loop term a 2 − 4ã 1 β 0 − 2β 1 (dashed magenta -lighter-line), and compares it with (α s /4π)a L 3 ln(C A α s /2) (solid magenta -lighter-line). We can see that the solid magenta line is never comparable in size with the two loop term. Therefore, in our present comparison with lattice, it is better not to include the leading ultrasoft resummation with the two-loop term. Instead it is more consistent to include it with the three-loop term, as we did in the previous paragraph. One could also consider including the subleading ultrasoft resummation in the analysis, since one of the reasons that motivated its computation was that the α 5 s ln α s terms, computed in Ref. [21], were found to be quite large. Nevertheless, given the behavior of the three-loop terms shown in Fig. 12, it can not be excluded that a similar compensation, in the range of distances we are studying in this paper, between the four-loop α 5 s ln α s terms and the unknownã 4 coefficient takes place. This means that if we were to include sub-leading ultrasoft resummation in the analysis, we would have to allow for generous perturbative uncertainties for the four-loop terms, to reflect for the aforementioned possible cancellations, and therefore the extraction of α s would not really benefit much.
We have, in any case, performed the analyses with leading-ultrasoft resummation included along with the two-loop term (N 2 LL accuracy), with leading-ultrasoft resummation included along with the three-loop term, and with sub-leading-ultrasoft resummation included along with the three-loop term (N 3 LL accuracy). These studies allow us to verify the expectations described in the previous paragraph. What we find is that: (i) the N 2 LL fits tend to give similar but a bit larger χ 2 values than the corresponding two-loop fits, (ii) the fits with leading-ultrasoft resummation at three loops tend to give χ 2 values that are very close to the three-loop ones, and (iii) the N 3 LL fits tend to give larger χ 2 values than the other fits with three-loop accuracy. This behavior is consistent with the observations in the previous paragraphs about the size of the ultrasoft terms. For these reasons, and as already mentioned above, we only present the numbers with leading-ultrasoft resummation included along with the three-loop terms in our final results in Sec. IV D. We use the value µ = 1.26r −1 1 ∼ 0.8 GeV, for the ultrasoft factorization scale, in all our analyses. We have checked that variations of µ only produce small effects on the results, as could be expected since the impact of the ultrasoft resummation itself was already small, as we just discussed.

C. Additional cross checks
Before presenting and collecting our final numbers in Sec. IV D below, we describe in this section several cross checks that we have performed. They serve to further verify the reliability of our procedure and results.
A first natural question is what would have happened if we had used our previous procedure in Ref. [5] with the current enlarged lattice data set. We have repeated that analysis, and obtained r 1 Λ MS = 0.48 ± 0.05 (at threeloop plus leading-ultrasoft-resummation accuracy). This number is perfectly compatible with the result quoted in Ref. [5], and has an error of similar size. The fact that the error did not really decrease, with respect to Ref. [5], may be seen as an indication that, to really benefit from the extended data set and the shorter distances reached, one should really employ expressions that completely avoid all ln(rν) terms, as we have done now. The previous figure is also compatible with the results obtained in the present paper with our modified analysis.

Comparison with lattice data for the force
One can also obtain lattice data for the force, and directly compare the perturbative expressions for the force with it. The problem in obtaining the lattice data for the force is that one needs to perform a numerical derivative of the lattice result for the energy, which is a non-trivial task. There are several uncertainties related to this numerical derivation, see Sec. II for details. Since we need to introduce smoothing splines, and use synthetic data, to compute the derivative and obtain the force and its error, the resulting data we have for the force does not fluctuate independently. For this reason the χ 2 analyses that we perform in our study may be seen as not completely adequate, or a bit unreliable. This fact led us to use the comparison with the energy, which does not suffer from this shortcoming, in our main analyses. It must nevertheless be checked what is the outcome if one compares with the force directly, which is what we do here.
To compare with the lattice data for the force, we can put the three data sets that we are using together, since the issue about the normalization, that affected the energy, is absent when we take the derivative. Apart from that, we follow the same procedure described in the pre-r 0.75r 1 r 0.75r 1 r 0.7r 1 r 0.7r 1 r 0.65r 1 r 0.65r 1 r 0.6r 1 r 0.6r 1 r 0.55r 1 r 0.55r 1 r 0.5r 1 r 0.5r 1 r 0.45r 1 r 0. vious sections, except that now we fit the perturbative expressions for r 2 F (r, 1/r) to the lattice data for the force. In Fig. 13 we show the values of the reduced χ 2 s for the different distance ranges. As we can see from the figure, for distances r < 0.6r 1 , which are the ones that passed our criteria when doing the analyses with the energy in the previous sections, the χ 2 s are below one, confirming therefore that we can perform the extraction of α s in these distance ranges. Actually, from Fig. 13, it would seem that an even larger distance range, something like r < 0.7r 1 , could be used. The results we obtain for r 1 Λ MS at three-loop accuracy, for r < 0.6r 1 and shorter distance ranges, are shown in Fig. 14. As we can see, the results are in very good agreement with what we obtained in the previous sections, compare with Fig. 10. This agreement is an important cross check of the result, since, we recall again that, it is only the slope of E 0 (r), which is what the force encodes, that matters for the α s extraction. We could have used this force data for our main analysis, but we preferred to leave it as a cross check because of the issues discussed above.

Analyses without lattice points with larger systematic uncertainties
The lattice points for E 0 (r) at the shortest distances have larger discretization errors, which is part of the reason why we could not obtain reliable lattice data for the force at these distances. The way we estimate these errors was described in detail in Sec. II. As explained there, our error estimates for these points should correctly reflect the size of the residual discretization errors in the data. Nevertheless, the points at larger distances, where the discretization errors are negligible, could be seen as more reliable than the ones at the shortest distances. With that in mind, we have repeated our analyses but dropping the first six points for each lattice spacing, which are the ones with larger systematic errors. We use N ref = 7 for this analysis. We find that the χ 2 curves present a behavior similar to that in Sec. IV B, in particular the region r < 0.6r 1 is suitable for the α s extraction. We show the results for r 1 Λ MS at three-loop accuracy obtained with this reduced data set in Fig. 15. Note that we can only perform the fits when there are at least three data points. In particular no fits with the β = 7.373 data set are possible in the r < 0.6r 1 region, see Fig. 5. The results in Fig. 15 are compatible and in good agreement with those of the previous sections. This shows that the inclusion of the points with larger systematic errors does not distort the result for α s , and can be seen as an indication that our estimation of the discretization errors is reliable.

Analyses with only points at shorter distances
One can also take the opposite point of view, with respect to the preceding section, and take for granted that our estimation of the systematic errors correctly reflects the residual discretization uncertainty. Then, the points at shortest distances should be the best ones for the α s extraction, since this is the most perturbative region. With this view in mind, we repeated the analyses but now including only the second to seventh point at each lattice spacing; note that we chose to still omit the first point at each lattice spacing, which corresponds to a distance r = a, since for this point the reliability of the error estimate may be more questionable. We use N ref = 2 for this analysis. There is only one distance range for each data set to perform fits in this case, i.e. that spanned by the second to seventh points. We find that the χ 2 s in this case are always below 1, which is perhaps not too surprising, since the errors for these points are larger than for the rest. We show in Fig. 16 the values of r 1 Λ MS obtained at three-loop accuracy in this case. As one can see from the figure, the results are compatible with those in the previous sections. Note that, in particular, they are compatible with the results in Sec. IV C 2 above, where we did the opposite than here, and omitted the first six points at each lattice spacing.

Analyses with extended distance ranges
For completeness, in this section we present the results one obtains when using larger distance ranges in the fits, up to r < 0.75r 1 . Recall that distances r < 0.6r 1 are the ones that passed our χ 2 criteria, and were therefore deemed as suitable for the α s extraction. The point of showing here the results from larger distance ranges is to illustrate that nothing dramatic happens beyond that point. Figure 17 shows the results for r 1 Λ MS at threeloop accuracy, in all the distance ranges we have analyzed. As one can see from the figure, the fits that use distances larger than 0.6r 1 give results for r 1 Λ MS that are compatible with those used in our main analysis. The error bars, which, remember, come from unknown higher- perturbative orders, are larger in the extended distance ranges. This may be attributed to the fact that those fits involve lower-energy scales and therefore larger values of α s .

Analyses of possible influence from non-perturbative-condensate terms
As discussed in previous sections, our χ 2 criteria made manifest that we can safely use perturbation theory to describe the lattice data in the r < 0.6r 1 distance ranges. Consequently, we used a purely perturbative expression in all our fits, and neglected any possible, parametrically suppressed, non-perturbative contributions to the static energy at these distances. Nevertheless, one might ask whether or not the presence of some non-perturbative term, which should necessarily be small, could distort the outcome of our fits in a significant way. To quantitatively address this question, we repeated the fits adding a monomial term, with a coefficient to be fitted, to our perturbative expression at three-loop accuracy. We considered r 3 and r 2 monomials, which could be associated with gluon and quark local condensates, and also an r monomial. Our new fits contain therefore two parameters, which are Λ MS and the coefficient of the monomial term. From the outcome of these fits, we do not find any evidence for a significant non-perturbative term in any of the cases. Furthermore, the values we obtain for Λ MS , at the r < 0.6r 1 distance ranges, in these modified fits are perfectly consistent with the outcome of our previous fits. For illustration, we show the results for Λ MS from the fits with an r 3 monomial, together with the outcome of our default fits, in Fig. 18. The red stars in the figure are the results from the fit including the r 3 term, and the points with error bars are the results from our default fits. We can clearly see from the figure that adding a monomial term, that could encode non-perturbative effects, does not distort our results for Λ MS . This fact provides further evidence that our extraction of α s is robust. It is worth stressing again that the goal of this section was to check that our perturbative analyses at short distances are stable and robust. A different question is how should one include non-perturbative effects to the perturbative expressions, in the form of local condensates, non-local condensates, etc., in order to extent their range of validity to larger distances. This is a very interesting question in itself, but it does not fall within the scope of the present paper, which is obtaining a precise determination of α s from the short-distance region of the static energy. We leave these studies for future work.

D. Final result for αs
In this section we collect the numbers we obtained from the analyses described in Secs. IV A, IV B and give our final result for α s . Table III summarizes the results we obtain at threeloop accuracy, for the r < 0.5r 1 fit range. For each value of the lattice spacing, we present the result for aΛ MS , for each of the three points we use to define N ref , in the second to fourth columns in the table. As described in the previous sections, we take the range spanned by these three N ref choices to get our result for aΛ MS at each lattice spacing. This is shown in the fifth column of the table, the central value here is the average of the three previous columns. We then convert the result to r 1 units using the values of r 1 /a given in Sec. II, and display it in the sixth column of the table. In all cases, the first error is the perturbative one, the second error is the one from statistical uncertainties in the fit, and, for the sixth column, the third error corresponds to the one coming from the value of r 1 /a. We added the errors in quadrature on the right-hand side of the sixth column. Finally, we take a weighted average, with the inverse of the total error as the weight, of the results for the three different β values to obtain our final number for r 1 Λ MS .
We added the errors linearly in this average, since the error at each β value is dominated by the perturbative uncertainties, which are common to all data sets. This final number is shown on the bottom-right corner of the table.
As explained in Sec. IV B 2, we also include here the results we obtain at three-loop plus leading-ultrasoftresummation accuracy. These results are shown in Tab. IV, with the same format as the previous table. As we can see, the final results for r 1 Λ MS in Tabs. III and IV are quite similar, as we expected from the discussion in Sec. IV B 2. We take the number at three-loop plus leading ultrasoft resummation accuracy as our best and final result. The final figure for r 1 Λ MS , which we recall is for n f = 3 light-quark flavors and uses lattice data in the 1.3 GeV -5.1 GeV energy range, is therefore We now convert this result to physical units by using the value r 1 = 0.3106 ± 0.0008 ± 0.0014 ± 0.0004 fm = 0.3106 ± 0.0017 fm, obtained in Ref. [32] from the pion decay constant f π (where we added the errors in quadrature on the right-hand side). We obtain where the first error corresponds to the one in Eq. (11), and the second to the value of r 1 , we added the two errors in quadrature on the right-hand side. This value of Λ MS gives α s (1.5 GeV, n f = 3) = 0.336 +0.012 −0.008 , which corresponds to α s (M Z , n f = 5) = 0.1166 +0.0012 −0.0008 , where we used the Mathematica package RunDec [33] to obtain the numbers in Eqs. (13)-(14) 1 . Equation (14) constitutes the main result of the paper. This number is perfectly compatible with our previous result in Ref. [5], which was α s (M Z ) = 0.1156 +0.0021 −0.0022 , and supersedes it.

V. DISCUSSION AND COMPARISON WITH OTHER WORKS
Having obtained our result for α s , we can now show how well the perturbative expressions describe the lattice data, when this number is used. To illustrate this,   we show, in the upper panel of Fig. 19, the lattice data for the smaller lattice spacing we used, i.e. β = 7.825, together with the perturbative prediction at three loop plus leading-ultrasoft resummation accuracy. The black curve corresponds to our central value for Λ MS , and the (barely visible in that scale) grey band is obtained by varying Λ MS within our obtained range. Note that the plot is in units of r 1 , but no error associated to normalization of the lattice data is shown, since that error is irrelevant when only one lattice set is used and the perturbative expression is normalized to it, as we do in this plot. To see more clearly how well the theoretical expression describes the lattice data, we also show, in the lower panel of the figure, the result of subtracting the perturbative prediction from the lattice data. The error bars in that plot are obtained by adding in quadrature the errors of the lattice data and the band in the theoretical expression (which was due to the variation of Λ MS ). Fig. 19 clearly shows that the perturbative expression can perfectly describe the lattice data in that distance range. (We recall that it is more difficult to reliably assess discretization errors for the point at shortest distance, see Sec. II.) In Fig. 20 we put together the data for all the lattice spacings we have, including those used in Ref. [5], i.e. from β = 6.664 to β = 7.825, and compare them with the perturbative expressions at different orders of accuracy. The uncertainties due to the normalization of the lattice data to a common scale are now included in the error bars, as it is appropriate when putting together data from different lattice spacings. To visualize more easily the differences among the different perturbative orders, we also show in Fig. 21 the ratio of the lattice data over the theoretical prediction. Finally, Fig. 22 shows, for completeness, the direct comparison of the perturbative expressions with the lattice data for the force presented in Sec. II.
It is worth stressing at this point that the knowledge of the static energy at three-loop accuracy was crucial in order to obtain our precise extraction of α s in Eq. (14). Although the fits at lower orders of perturbative accuracy can already give low χ 2 values, see Figs. 6-8, the uncertainty that one should associate to a value of α s extracted at lower orders is much larger. This was clearly illustrated in Fig. 9, where the dramatic reduction of the error bars when increasing the perturbative accuracy of the prediction is manifest. We also recall that the different perturbative orders always refers to the corrections to the potential, i.e. theã i coefficients in Eq. (9), but we always use four-loop accuracy for the running of α s . If one does not do that, and uses the running of α s at lower accuracies, the χ 2 values resulting from lower orders in perturbation theory are much higher.
Let us now compare our result for Λ MS with other recent extractions of the strong coupling. The present analysis, together with our preceding paper [5], consti-     tutes, at present, the only extraction of Λ MS from the QCD static energy with at least three flavors; therefore, the only one that can be used to obtain α s (M Z ). Other analyses aiming at extracting Λ MS from the static energy with less than three flavors include Refs. [23,30,31,[34][35][36][37]. In particular, let us mention that the n f = 2 analysis of Ref. [23] concludes that smaller lattice spacings than those currently available to them would be needed to extract Λ MS . This reference uses Wilson fermions. The data for the force that they use contains only three points below r = 0.75r 1 , with none below r = 0.5r 1 . Our analyses show that we would not have been able to obtain Λ MS with that amount of data. In this sense, we do agree with the findings of Ref. [23]. Regarding lattice α s extractions from other observables, the FLAG collaboration recently presented, in Ref. [38], a comprehensive and critical review of all the available α s lattice determinations, and provided an average. We show this lattice We include: The lattice determinations of the HPQCD [3], PACS-CS [41] (we take the number quoted in Ref. [38]), and ETM [42] collaborations. Non-lattice determinations from τ decays in Refs. [43,44]; from thrust in e + e − collisions in Refs. [4,45]; the recent H1 re-analysis of Ref. [46]; and the PDF-fit ABM13 result of Ref. [47] (note that the error bars in this case do not include effects from higher unknown perturbative orders). We also show the lattice average by the FLAG collaboration [38] and the PDG average excluding lattice results [2] -black points-.
average, together with our new result, in Fig. 23 (note that the FLAG average includes the result from Ref. [5]). We also show in the figure a few other individual lattice determinations of α s , a few selected recent non-lattice determinations, and the PDG average excluding lattice results [2]. Further determinations of α s , as well as discussions about them, can be found, for instance, in the summary reports of recent dedicated workshops [39,40].

VI. SUMMARY OF RESULTS AND CONCLUSIONS
We have improved our previous extraction of α s , in Ref. [5], from the comparison of lattice data with perturbative expressions for the static energy of a heavy quark-antiquark pair. This has been possible because a considerable amount of new lattice data at shorter distances has become available [6], which has allowed us to carry out an extra number of cross checks and hence to considerably reduce the systematic errors. In particular we have been able to correct for cutoff effects in the shorter-distance points, to analyze the dependence on the fit range, and to carry out separate analyses for different lattice spacings. Thus we could, for instance, discard points which are not in the perturbative regime, points which suffer from large cutoff effects and, very importantly, avoid the lattice normalization errors that dominated our previous extraction. On the other hand, we have used improved perturbative expressions, in which not only the first renormalon is avoided, and the ultrasoft logarithms resummed, but also the soft logarithms are summed up. This appears to be necessary because the lattice data covers a relative large range of distance values now. Our final result reads which corresponds to α s (M Z , n f = 5) = 0.1166 +0.0012 −0.0008 .
This updated result reduces the errors by roughly a factor of two with respect to our previous extraction. It displays a higher central value, which is, nevertheless, perfectly compatible with our previous result. The errors of the α s extraction presented here can in principle be reduced by just incorporating more lattice data at shorter distances, with no further modification of the perturbative expressions, which are already at the three-loop level. Notice that with the lattice data available at present, there is still little sensitivity to the ultrasoft resummation, and hence we do not expect much sensitivity to the yet unknown four-loop contribution. We have also checked that there is no sensitivity to other possible non-perturbative effects, like for instance those due to gluon or quark condensates.
We conclude that the method first outlined in Ref. [31], and further developed in the present paper, is not only able to produce competitive extractions of α s when confronted with realistic lattice data [5], but also to consistently improve on the outcome upon the incorporation of new shorter-distance data. For the future, to further corroborate the result obtained here, it would be important to determine α s from the static energy calculated using other lattice actions. ACKNOWLEDGMENTS X.G.T. thanks Rainer Sommer for a useful conversation during a visit to DESY Zeuthen, which partially triggered the present updated analysis. This work was supported in part by U.S. Department of Energy under Contract No. DE-AC02-98CH10886. The work of X.G.T. is supported by the Swiss National Science Foundation (SNF) under the Sinergia grant number CRSII2 141847 1. The work of N.B. and A.V. is supported in part by DFG and NSFC (CRC110) and they acknowledge financial support from the DFG cluster of excellence "Origin and structure of the universe" (www.universe-cluster.de). J.S. is supported by the CPAN CSD2007-00042 Consolider-Ingenio 2010 program (Spain), the 2009SGR502 CUR grant (Catalonia) and the FPA2010-16963 project (Spain). X.G.T. is grateful to the Mainz Institute for Theoretical Physics (MITP) for its hospitality and its partial support during the completion of this work.