Determination of (cid:11) ( M z ) from an hyperasymptotic approximation to the energy of a static quark-antiquark pair

: We give the hyperasymptotic expansion of the energy of a static quark-antiquark pair with a precision that includes the e(cid:11)ects of the subleading renormalon. The terminants associated to the (cid:12)rst and second renormalon are incorporated in the analysis when necessary. In particular, we determine the normalization of the leading renormalon of the force and, consequently, of the subleading renormalon of the static potential. We obtain Z F 3 ( n f = 3) = 2 Z V 3 ( n f = 3) = 0 : 37(17). The precision we reach in strict perturbation theory is next-to-next-to-next-to-leading logarithmic resummed order both for the static potential and for the force. We (cid:12)nd that the resummation of large logarithms and the inclusion of the leading terminants associated to the renormalons are compulsory to get accurate determinations of (cid:3) MS when (cid:12)tting to short-distance lattice data of the static energy. We obtain (cid:3) ( n f =3) MS = 338(12) MeV and (cid:11) ( M z ) = 0 : 1181(9). We have also found strong consistency checks that the ultrasoft correction to the static energy can be computed at weak coupling in the energy range we have studied.


Introduction
The static potential or, more precisely, the energy of a static quark-antiquark pair separated by a distance r, is one of the objects most accurately studied by lattice simulations. This is due to its relevance in order to understand the dynamics of QCD. On the one hand, it is a necessary ingredient in a Schroedinger-like description of the Heavy Quarkonium dynamics. On the other hand, a linear growing behavior at long distance is signaled as a proof of confinement. Moreover, throughout the last years, lattice simulations with dynamical fermions have improved their predictions at short distances, see for instance [1][2][3][4][5][6][7][8][9]. In addition, the accuracy of the perturbative prediction of the static potential has also improved significantly over the years [10][11][12][13][14][15][16][17][18][19][20][21][22]. The precision reached nowadays is nextto-next-to-next-to-leading order (NNNLO) for fixed order computation and next-to-nextto-next-to-leading logarithmic (NNNLL) order for renormalization group (RG) improved computations. 1 The combination of these two items: high order perturbation theory and lattice data at short distances, potentially allows quantitative comparison between perturbation theory and lattice simulations. Nevertheless, naive comparison between lattice data and perturbative results may lead to strong disagreement depending on how the perturbative expansion is implemented in practice. Such different behavior can be understood [23] on the basis JHEP09(2020)016 of a renormalon based picture. Overall, analyses that correctly implement the renormalon cancellation (even if it is not mentioned explicitly) show reasonable agreement between lattice data and the corresponding perturbative expressions. This was something that was observed around twenty years ago, see for instance [18,[23][24][25][26][27]. Nowadays, the more recent unquenched data and the knowledge of higher orders in perturbation theory have allowed to obtain competitive determinations of Λ (n f =3) MS , like those in [4,7,9,28,29]. Nevertheless, none of them have implemented the resummation of large ultrasoft logarithms to NNNLL order. This is something that we will do in this paper. In particular, expressions for the force with NNNLL precision are given for the first time.
A most relevant, and novel aspect, of this paper will be the use of hyperasymptotic expansions to deal with the asymptotic behavior associated to the renormalons. Hyperasymptotic expansions were first introduced in the context of the asymptotic solutions to ordinary differential equations [30,31]. We use here its generalization to quantum field theories with a running coupling constant and renormalons developed in [32][33][34]. Such an approach allows us to compute observables with a well-defined power counting and reach exponential accuracy. In [32], the general formalism was explained in great detail, and the static potential in the large β 0 approximation was used as a test-case and analyzed in great detail as well. In [33] ultraviolet renormalons were explicitly included in the formalism, and the pole mass, and other related objects, were studied in great detail. The general counting for the truncation of the hyperasymptotic approximation was also given in this reference (see also [34] where a summary of such counting is given). In this paper, we give the hyperasymptotic expansion of the energy of a static quark-antiquark pair with a precision that includes the effects of the subleading renormalon. The terminants associated to the first and second renormalon are incorporated in the analysis when necessary. In particular, we determine the normalization of the leading renormalon of the force and, consequently, of the subleading renormalon of the static potential. Note that the incorporation of the effects associated to the subleading renormalon of the static potential (which is the leading one of the force) is novel compared with the analyses of [4,7,9]. On the other hand [28,29] use a different method to handle the renormalon singularities, and the subleading renormalon has also been studied in the situation when α/r Λ QCD .

Hyperasymptotic expansion of the static energy
The energy of a static quark and a static antiquark in a colour singlet configuration separated by a distance r, E(r), admits an operator product expansion using pNRQCD [35,36]: E(r) = V (r; ν us ) + δE us (r; ν us ) .
(2.1) V (r; ν us ) only encodes the physics associated to the scale 1/r. δE us (r; ν us ) encodes the physics associated to scales smaller than 1/r and appears at O(r 2 ) in the multipole expansion.
As explained at length in [32] we need to give meaning to the different terms of the operator product expansion with nonperturbative power accuracy. In particular this implies that we have to regularize the perturbative series expansion of V . We define it using the JHEP09(2020)016 principal value (PV) prescription, i.e. the average of the Borel integrals infinitesimally above and below the cut in the real plane. See [37] and the appendix of [32] for explicit expressions. We then have to use the same prescription for δE us . Therefore, we have E PV (r) = V PV (r; ν us ) + δE PV us (r; ν us ) . (2. 2) The asymptotic expansion of V PV coincides, by construction, with the asymptotic expansion of V and reads 3) The natural scale for this expansion is ν s ∼ 1/r. The coefficients n a n ν s r; ν us ν s (2.4) were computed in [10][11][12][14][15][16] for n = 0, 1, 2, 3. For ease of reference we quote them in appendix A. The perturbative expansions of the pole mass, and of the static potential in the large β 0 approximation, are characterized by having a single scale (leaving aside Λ QCD ): the heavy quark mass renormalized in the MS scheme, m, for the case of the pole mass, and 1/r for the case of the static potential in the large β 0 approximation. In other words, they are infrared finite at any finite order in perturbation theory. This is not so for V PV . This reflects in that V PV is logarithmically infrared divergent. Such behavior first appears in V 3 (see eq. (A.1)), which endures a linear ln(ν us ) dependence that is not absorbed in the renormalization of the strong coupling in V PV . This logarithmic behavior is cancelled instead by δE us . Therefore, both V PV (ν us ) and δE us (ν us ), unlike the pole mass, are renormalization scale and scheme dependent in perturbation theory. Such scale and scheme dependence cancels in the sum in eq. (2.1).
Powers of α(ν s ) ln νs νus (we take ν s ∼ 1/r) can be resummed using RG techniques. They have been computed in [17] with leading logarithmic (LL) accuracy and in [18] with next-to-leading logarithmic (NLL) accuracy (see also [19]). This produces the following correction to V PV : δV RG (r; ν s , ν us ) = r 2 (∆V ) 3 G(ν s ; ν us ) , (2.5) where (V o is the potential of a static quark-antiquark pair in a colour octet configuration, and
(2.8) δE us (r; ν us ) encodes the physics associated to scales smaller than 1/r and can be computed using the multipole expansion. At O(r 2 ) its explicit expression reads (in the Euclidean) This quantity has a different behavior depending on the relative size between Λ QCD and ∆V ∼ α(ν s )/r. If Λ QCD α(ν s )/r the above expression can be approximated to On the other hand, if Λ QCD ∆V , δE us can be computed at weak coupling as an expansion in powers of α(ν us ). It has the following scaling δE us (r; ν us ) ∼ r 2 (∆V ) 3 H(α(ν us )) , (2.11) where ∆V is generated dynamically and H(α(ν us )) admits a perturbative expansion in powers of α(ν us ) (up to logarithms). δE us (ν us ) is known to order r 2 (∆V ) 3 α 2 (ν us ) ∼ 1 r α 5 in the MS scheme: 2 where we take the next-to-leading order (NLO) expression from [21]. This is one order more than we have for V PV : ∼ 1 r α 4 . We remind again that δE us is scheme dependent. Finally, if both scales, Λ QCD and ∆V , are similar in size, δE us (r; ν us ) is an unknown function of the ratio of these two scales.
We now focus on the hyperasymptotic approximation to V PV . The leading asymptotic behavior of V n is known (and up to a factor minus two is equal to the asymptotic behavior of the pole mass [39]). It reads JHEP09(2020)016 . The coefficients c k are pure functions of the β-function coefficients, as first shown in [40] for the case of the pole mass. They can be found in [41][42][43]. At low orders they read (c 0 = 1) where the s n are defined in [32]. We construct the hyperasymptotic expansion of V PV along the lines of [32,33]. It does not have ultraviolet renormalons, whereas the leading infrared ones are located at dimension 1 and 3 (i.e. at u = 1/2 and at u = 3/2 in the Borel plane). The termination of the perturbative series associated to these renormalons produces nonperturbative power contributions of order Λ QCD and r 2 Λ 3 QCD respectively: V n α n+1 (ν s ) ; (2.17) and Approximating V PV by V P corresponds to achieve superasymptotic approximation. In the generic labeling (D, N ) of the truncations of the hyperasymptotic approximation defined in [33,34], it corresponds to (0, N P ) precision. The next order in the hyperasymptotic approximation is labeled as (1,0) and means adding Ω V 1 /r to V P . Its explicit expression reads (in [32] Ω V 1 was named Ω V , and Z V 1 was named Z V )

JHEP09(2020)016
Note that the u = 1/2 renormalon does not cancel with the renormalons in δE us . Eq. (2.19) exactly corresponds to −2Ω m defined in [33], since the renormalon of twice the pole mass cancels with the renormalon of the static potential (in other words Z V 1 = −2Z m ). Only after the inclusion of the pole mass in real observables, this renormalon cancels. This is the reason we have to specify the prescription used to regularize the perturbative sum of E. Had we included 2m OS , the summation scheme dependence would have disappeared.
We now move beyond (1,0) hyperasymptotic precision by including the third term in eq. (2.16). One then generically reaches hyperasymptotic precision (1,N ). For N large one would start to be sensitive to the next renormalon, which then has to be considered. We can use eq. (2.10) to determine the renormalon structure of the subleading infrared renormalon. Due to the fact that V A = 1 with NLL accuracy, we can determine the first two terms of the asymptotic expansion of the perturbative series associated to the u = 3/2 renormalon of the static potential (we can do similarly for the force, see eq. (4.1)). The best way to quantify the asymptotic behaviour of the perturbative series is by performing its Borel transform: (2.24) The Borel transform will have a singularity, due to the dimension d = 3 non-local condensate, at u = d/2 = 3/2: where b 1 = ds 1 . This singularity produces the following asymptotic behavior 3 This asymptotic behavior has associated the terminant Ω V 3 , which reads .
If we have enough terms in perturbation theory to be sensitive to the u = 3/2 renormalon, we may reach order (3, 0) precision in the hyperasymptotic counting for V PV (see [33,34] for a detailed account of the hyperasymptotic counting we use in this paper). This is the maximal accuracy we will test in this paper.

JHEP09(2020)016
In principle, we will work under the hypothesis that we have enough terms in perturbation theory to be sensitive to the u = 3/2 renormalon and that Λ QCD α/r so that δE PV us can be computed in perturbation theory, though we will also test predictions where δE PV us is modeled by a Λ 3 QCD r 2 term. In any case, the final expression we obtain reads 31) and N max = 3. We will confront the above theoretical expression with nonperturbative evaluations of E(r). E(r) can be determined accurately using Montecarlo simulations in the lattice of the quantity where W is the rectangular Wilson loop with edges x 1 = (T /2, r/2), x 2 = (T /2, −r/2), y 1 = (−T /2, r/2) and y 2 = (−T /2, −r/2). The symbol means the average over the massless gluons and the light quarks. This quantity is linearly divergent in 1/a by an r independent constant. Therefore E PV (r) and E latt (r) are equal up to an additive rindependent constant (and up to O(a) lattice artefacts which can be r dependent). Leaving aside these lattice artifacts, E PV (r) and E latt (r) have the same r dependence.

Hyperasymptotic expansion of the force
We now consider the force, which we define as the derivative of the static potential: where the present known values of f n read (we take ν s = x s /r, and for n ≥ 3, f n also depends on ν us : f n ν s r; νs νus ) F PV admits a strict perturbative expansion in powers of α(ν s ) (the dependence in ν us is hidden in the coefficients f n ). On the other hand F PV (r; ν us ) is not RG invariant, since it is dependent on ν us .

JHEP09(2020)016
We also consider the quantity where (at leading order and assuming Λ QCD ∆V ) d dr δE PV us (r; ν us ) = C F r(∆V ) 3 α(ν us ) 9π 6 ln ∆V ν us + 6 ln 2 + 1 . If we neglect renormalons, F (r) is now known with N 3 LO precision. Adding all the terms the resulting expression is equal to eq. (10) of [4]. Notice, though, that the definition of the coefficient a 3 is different here and in that paper. This is compensated by the last term in eq. (3.3) which is also different. We can also make a RG improved version for the force and for F (r): The pure perturbative expression for F RG (r) is known at the NNNLL level. Our expression corrects eq. (11) of [4] at the NNNLL level. We emphasize though that such eq. (11) was not used for analyses in this reference but rather the same expression we have obtained here. 4 Finally, let us emphasize that there is no unique way to decompose F RG (r). In the above result we have taken ν s and ν us to be r-independent in the derivative, though to resum the large ultrasoft logarithms we have to take ν s ∼ 1/r and ν us ∼ C A α(νs)

2r
. We can make this dependence explicit: ν s = x s /r and ν us = x us . This introduces an explicit dependence on r in ν s and ν us . We show how the different terms of F RG (r) look like in this situation in appendix B.
We now introduce in the discussion renormalon effects. F PV does not have the renormalon at u = 1/2. The leading renormalon is located at u = 3/2. On the other hand, F (r) can be considered to be an observable. The uncancelled renormalon at u = 1/2 that exists in E PV (r) vanishes in F after taking the derivative with respect to r. All other renormalons of V PV (r; ν us ) cancel with the analogous renormalons of δE PV us (r; ν us ) in E PV (r) and consequently the same cancellation takes place in F (r).

JHEP09(2020)016
The perturbative expansions of V and δE us are series in powers of α evaluated at different scales: α(ν s ) and α(ν us ) respectively. The same thing applies to F and d dr δE us . This makes that there is no renormalon cancellation order by order in α. If the perturbative series reaches orders high enough to be sensitive to the u = 3/2 renormalon, we should indeed incorporate the associated nonperturbative contribution to the PV summation, i.e. the corresponding terminants. The complete expression then reads where and d dr δE P us (r; ν us ) = Nus n=0 p n α n+1 (ν us ) (3.10) are the superasymptotic approximations of F PV and d dr δE PV us (r; ν us ). For them we have N F = 3 2π β 0 α(νs) (1 − c F α(ν s )) and N us = 3 2π β 0 α(νus) (1 − c us α(ν us )). We will usually take N F = 3 and fine tune the coefficient c F accordingly. For N us we know less orders of the perturbative expansion, and the scale ν us is small. Therefore, we will usually take N us = 0 and fine tune the coefficient c us accordingly.
As N F is small, there is some degree of fine-tuning between the relative size of 1 and c F α. This is not a problem because to determine whether the expansion in the terminant in eq. (2.27) makes sense, we do not have to look to c F but rather to the size of the complete correction. In other words, to see whether the subleading term K 1 α X is smaller than one (except if, for some reason, the leading order is anomalously small). In the range of ν s we cover with the data set I we use in section 5, this term is in the range −0.055 <K (P ) X,1 α X (ν s ) < 0.11. Therefore, it is safely small. Even for the ultrasoft case, where the situation is potentially worse, we find −0.29 <K (P ) X,1 α X (ν us ) < 0.0006 for the range of scales of the data set I. We find again that the correction is safely smaller than one. As an extreme test, we could even set this subleading correction to zero (setK (P ) X,1 = 0). We find that the fit shifts by 2 MeV only. On top of that one could be worried about a possible asymptotic character of the weak coupling expansion of eq. (2.27). At present we cannot make definite statements about the possible asymptotic nature of this expansion (see also the discussion in appendix A of [33] and [37]) but from the previous numerical discussion we do not see signs of divergence of the weak coupling expansion.
The last items of eq. (3.8) are the terminants of the soft and ultrasoft perturbative series. For them we have 1 Notice that Ω F 3 (ν s ) and Ω F 3 (ν us ) are different, not only because of the different renormalization scale each of them uses, but also because we truncate the perturbative expansion at different orders in F and d dr δE us . Therefore, even if we set ν us = ν s , the terminants will not cancel each other in general (they would only do if we truncate the perturbative expansions of F and d dr δE us to the same order).

JHEP09(2020)016 4 Normalization of the u = 3/2 renormalon
The hyperasymptotic expressions derived in the previous sections are completely determined up to the normalizations of the u = 1/2 and u = 3/2 renormalons: respectively. They can only be computed approximately. For Z V 1 we use the value determined in [43] using the static potential: Z V 1,MS = −1.1251(520). The direct determination of Z V 3 from V is complicated because of the u = 1/2 renormalon. On the other hand, the force is not contaminated by the u = 1/2 renormalon. Therefore, it is an ideal place where to see if the perturbative series, as we know it at present, is sensitive to the subleading (infrared) renormalon, which is located at u = 3/2 in the Borel plane.
The asymptotic behavior of the coefficients of the force read: and By considering the ratio of the exact and asymptotic expression we can obtain an approximate determination of the normalization of the u = 3/2 renormalon. We obtain We determine these numbers taking the central value of the ratio of f 3 /f (as) 3 Z F 3 at the scale of minimal sensitivity (see figure 1), which are x ≡ µr = 1.30, and x = 1.52 for n f = 0 and 3 respectively. For the error estimate we explore different possibilities: we vary µr by multiplying and dividing the central value by √ 2. This is the first error quoted in eqs. (4.3) and (4.4). We also consider the difference between f 2 /f (as) 2 Z F 3 at the scale of minimal sensitivity. This is the second error quoted in eqs. (4.3) and (4.4). We also estimate the importance of subleading 1/n corrections by considering the difference of including the 1/n term or not in eq. (4.1). This is the third error quoted in eqs. (4.3) and (4.4). Finally, we also explore the importance of the ultrasoft associated terms (as they should not affect, or little, the determination of the normalization of the renormalon). The error associated to ultrasoft effects is estimated by eliminating the last term in the second line of a 3 in eq. (A.1) and the last term in f 3 in eq. (3.2). The variation is indeed small, as we show in the last error item in eqs. (4.3) and (4.4). The first and second error (and to some extent the third) are somewhat redundant, as they both measure the fact that n = 3 is still finite. Still, for the total error, we combine all them in quadrature and make the variation symmetric around the central value. This indeed yields a conservative estimate of the error, as we can see in figure 1. In figures 1 (a) and 1 (c), we can see the dependence of Z F 3 , i.e. of f n /f (as) n Z F 3 , with respect µr for different values of n. Around the scale of minimal sensitivity they are inside the error band, even for a coefficient as low as We profit to give determinations for other values of n f using the same error JHEP09(2020)016  analysis. They can be found in table 1, where we also give estimates of the higher order coefficients of the perturbative series of the force.
Our conclusions are different from those in ref. [44], where it was concluded that it was not possible to determine the normalization of the u = 3/2 renormalon. The authors used the function following the method proposed in refs. [45,46], and first quantitatively applied to the JHEP09(2020)016 leading renormalon of the pole mass and static potential in ref. [42]. D (N ) F (u) is singular but bounded at the first IR renormalon. Therefore, we can estimate Z F 3 from the first coefficients of the series in u, using We plot the predictions for different orders N in the figures 1 (b) and 1 (d) for n f = 3 and n f = 0 respectively. We observe that the convergence is worse than for the determination of Z F 3 using f n /f (as) n Z F 3 . This was also observed very clearly in [47] for the energy of an static source. In that case, and this case here, we observe convergence but at a slower pace. Actually, the NNLO and NNNLO predictions are well inside the error band of our predictions in eqs. (4.3) and (4.4), though less precise, as the variation between different orders is bigger than in the previous case. Compared with the analysis in [44], we make the analysis at larger values of x, but close to one, where we find stability. For the method using eq. (4.6) stability is found for x ∼ 2. For this method, working with x = 1 does not yield a convergent series. This may explain the conclusions reached in [44].

Fit of α
We will now compare our theoretical expressions for the static energy with recent lattice data obtained with n f = 3 active flavours. The static energy computed in the lattice can be equated with the theoretical expressions we have up to a constant. Therefore, we will always use the following equality: In principle the analysis should not depend on the value of r ref we use in this equation. In practice there will be some dependence (we will check this dependence later). By default we will take the value r ref = r min , i.e. the point at the shortest distances we use.
is renormalon free. Nevertheless, there are different ways to implement this cancellation which are not equally efficient. The use of F (r) seems optimal in this respect. On the one hand the leading renormalon identically vanishes. The subleading renormalons of the static potential cancel with those of the ultrasoft energy. This cancellation takes place order by order in α if both quantities are expanded in powers of α evaluated at the same scale, and if both perturbative expansions are truncated at the same order. This is something that we will not do, as the perturbative expansion of the static potential and the ultrasoft energy are known to different orders, and the natural energy scales in F PV and d dr δE PV us (r; ν us ) are different. This last issue also reflects in that setting ν s = ν us misses the resummation of large logarithms associated to the ultrasoft scale. These can be important, and not incorporating them can jeopardize the convergence of the perturbative series. Therefore, we also perform the resummation of logarithms, and our default expression will be the RG improved expression. In this case the precision we have is LL, NLL, NNLL and NNNLL if we do not include renormalon effects. From the analysis performed in section 4, we have seen that perturbation theory of the force JHEP09(2020)016 has reached high enough orders to be sensitive to the d = 3 renormalon. Thus, to the NNNLL expression, we will add the terminants associated to the d = 3 renormalon. We will name this approximation NNNLL hyp . In the hyperasymptotic counting, this means that the maximal accuracy that we will seek in F will be (3, 0). We will take as default that the asymptotic behavior associated to the d = 3 renormalon is reached for N = 3 for the force. For the ultrasoft term we will then take N = 0, since we will only incorporate one term of the perturbative expansion of the ultrasoft term. The different order at which we truncate the two perturbative expansions, and the different scale they depend on, make the terminants associated to F (r) and d dr δE us to be different. Obviously, the difference of determinations using NNNLL or NNNLL hyp will allow us to see the impact of including the terminants. 5 In order to compare with analyses where the resummation of logarithms is not incorporated, we will also perform computations with ν s = ν us . In this case, if we neglect renormalons, we can compute the observable with LO, NLO, NNLO and NNNLO accuracy (note that LO=LL and NLO=NLL), accordingly to the order in α we truncate the perturbative series of F (r). To make the connection smooth with the RG improved expression, we also incorporate at NNNLO the terminants associated to the u = 3/2 renormalons of F and d dr δE us . Note that for this last term we will still use N = 0, even though the renormalization scale of α is bigger. We will discuss this issue later.
To compare with the lattice results, we have to integrate over r: and to fit the outcome with the lattice data to determine Λ QCD . As we have mentioned above, our default fit is made using the RG improved expressions. For the central values of the renormalization scales, we take (ν s , ν us ) = (1/r, C A α(ν s )/(2r)). In the following we perform such fit and quantify the different sources of error. Note that F RG (r ) has to be introduced in the integral order by order in α in order to implement the renormalon cancellation in eq. (5.2), as it was first explained in [23]. Eq. (5.2) was originally used in [24], but its use for competitive determinations of Λ MS was first made in [4].
Dependence on the data points. For the fits, we use the lattice data of [9] (supplemented with the information given in [6]), which has made an updated error analysis of the data of [5]. Of these data points we only consider those obtained with β = 8.4, as they correspond to the shortest distances available: 1/a 8.3 GeV. In this ensemble the strange quark mass has been fine tuned to its physical value, and the pion mass gets the value 320 MeV in the continuum. This is only statistically significant 6 for r > 0.4r 1 ∼ 1/1.6 GeV −1 (see [8]).
In the fits, we will approximate the light quark masses to zero. The uncertainty associated with fixing the physical units of the parameter r 1 was seen in [9] to be comparatively 5 Alternatively, we have also performed fits changing the order at which we start including the terminant in the static potential from three to two. We indeed find the variation to be small. For a fixed scale ν, the change of the order at which we start including the terminant is implemented by changing the value of cF . We then find that our fits have an small dependence on cF . 6 We thank J.H. Weber for informing us of this.

JHEP09(2020)016
small compared with other uncertainties. Therefore, we will neglect it in the following. It was also observed in this reference that the effect of the correlation of the points to the final error was small. Thus, we also neglect this source of error. The discretization errors depend on the size of the parameter r/a. They have been studied in detail in [9], where it was concluded that, for r/a ≤ √ 8, tree-level improvement was enough to bring the discretization errors down to the point that they were smaller than the statistical errors and could, in comparison, be neglected. Therefore, we will use tree-level improved data and disregard the lattice data at shortest distances (for r/a ≤ √ 8), as well as the special geometry r/a = √ 12. This corresponds to one of the methods followed in [9] to account for discretization effects. This means that the shortest distance we consider is r min = 2.827 a, which in physical units reads r min = 0.353 GeV −1 . We have also compared with the older unquenched data of [1]. Overall, we observe the same qualitative features, though the lattice errors are bigger. Therefore, we will only present quantitative analyses with the data of [9].
To test the sensitive of the fit to the data we consider different ranges of data (similarly as it was done in [4]). We consider the following ranges: Set I: 0.353 GeV We show the result in figure 2. We observe the following. The dependence of the value of Λ MS on the range of the data set is very small. Obviously, as we increase the number of points, the statistical errors get smaller. This small dependence holds irrespectively of the order in the approximation for the theoretical expression used. We only see very small differences at NNNLL and NNNLL hyp order between the value obtained from the data set I and the rest (within one sigma for the statistical error, which is the only one we display in figure 2), and basically vanishing between the data sets II, III and IV.
To see how reliable the results are, we study the reduced χ 2 obtained with each data set. For the data sets I and II, the fit yields χ 2 red ∼ 0.5 to all orders in the hyperasymptotic expansion. Therefore, there is no significant dependence on the number of data points. For the data Set III, there is a mild increase: χ 2 red ∼ 0.5 − 0.6 but still well below 1. It is when we consider data set IV, which includes points down to 1/r = 1 GeV, that we see a significant increase in the χ 2 red . The magnitude of this increase, however, depends on the order, and, even in the worst case, it is not much bigger than 1. The LL and NLL fits yield χ 2 red slightly below 1, with a slight increase when going from LL to NLL. The χ 2 red reaches the maximum, 1.46, at NNLL. Since then higher order fits improve the quality of the fit and, significantly, the NNNLL hyp fit, which adds the terminant of the u = 3/2 renormalon, yields a χ 2 red = 0.67, similar to the other data sets. The Set IV is the more sensitive to the infrared, as it goes down to 1/r ∼ 1 GeV. This may reflect in a larger sensitivity to ultrasoft associated physics, which will then need to be described more accurately. This matches with what we see for the χ 2 red with set IV: LL, NLL, NNLL and NNNLL show a bigger χ 2 red (we emphasize, though, that they are still of order 1), which then goes down to a value similar to the one obtained with the other data sets after the inclusion of the terminants. This result is not trivial. We expect more sensitivity to infrared physics with the data set IV. What is not trivial is that this larger sensitivity to the infrared can be well described by our weak-coupling analysis. Indeed it is surprising that the ultrasoft effects do not blow up in any of the fits, since the α(ν us ) is evaluated at a rather low scale. For illustration (to produce these numbers we take Λ

JHEP09(2020)016
For this last data set, the very last points with smaller energy reach a regime where α(1/r) grows faster than 1/r, so that ν us grows for them. Therefore, their inclusion in fits should be taken with caution.
Overall, by only looking at the χ 2 red , we do not have a clear signal of which data set to use and, indeed, the fits yield similar numbers and χ 2 red at NNNLL hyp . Therefore, we will use set I as it is less sensitive, in principle, to long distances, though, as we said, the χ 2 red of the fits does not give a clear signal of a deterioration of the quality of the fit (something that one would expect if our perturbative approximation were not a good approximation to the data). In this respect note that the data set IV, which is more sensitive to the ultrasoft scale, yields a good χ 2 red after the introduction of the terminants associated to the u = 3/2 renormalon. We interpret this as an indication that the ultrasoft effects can be well described by a weak-coupling computation even at scales as low as ν s ∼ 1 GeV.
Another motivation to use the data set I is that the β = 8.4 ensemble suffers from frozen topological charge in the Montecarlo evolution. It has been shown (see [8]) that the effects of frozen topology in different sectors are statistically irrelevant for r < 0.4r 1 ∼ 1/1.61 GeV −1 . Therefore, by using the data set I this problem is completely avoided. On top of that, as mentioned above, the effects due to finite light-quark masses are not statistically significant for this energy range.
In order to see the quality of our fit, we also compare our theoretical expression using the values of Λ (n f =3) MS obtained from the fit with the lattice data. It is customary to compare directly with the potential (this can be done after fixing a normalization constant K that we fix below). The comparison is very good in the whole range we compare (up to 1 GeV), as we can see in the upper panel of figure 3 (note that we plot as a function of 1/r, a plot in terms of r, as it is customarily done, is even less precise). Nevertheless, such comparisons do not allow us to see the fine details due to the dependence in powers of r of the potential. For such comparison, it is better to define and we adjust K such that most of the r dependence vanishes. We show the comparison in the lower panel of figure 3. It is remarkable that pure perturbation theory predicts very well the data down to 1 GeV. The error band perfectly encodes all the data. This means, in particular, that with the precision of our computation we do not see any trace of nonperturbative effects down to scales 1/r ∼ 1 GeV.
Dependence on ν s . We now test the sensitivity of the fit on ν s . We will mainly work with the data set I, with which we can do variations of the parameters without entering in the regime where perturbation theory breaks down. We will try to vary ν s but keeping ν us constant. Our central value for ν us is ν us = C A α(ν s )/(2r). For the data Set I this yields values around ν us = 1 GeV. Therefore, besides ν us = C A α(ν s )/(2r), ν us = 1 GeV will be the other choice we take for ν us . We observe that both fits, with (ν s , ν us ) = (1/r, C A α(ν s )/(2r)) and with (ν s , ν us ) = (1/r, 1 GeV) yield very similar results. We show the outcome in   For the variation of ν s we take ν s = x s /r within the range x s ∈ [1/ √ 2, 2]. The lower limit of ν s is chosen to avoid reaching scales too low for our weak coupling analysis to break down. This happens, both if we take (ν s , ν us ) = (1/(2r), 1 GeV) or (ν s , ν us ) = (1/(2r), C A α(ν s )/(2r)). As we increase ν s , things significantly improve, and a safe value to take as a lower limit is x s = 1/ √ 2. We show the results of the fits for the central and extreme values of the parameter in figure 4.
We now compare with the fits with ν s = 2/r. It is interesting to see the different behavior of the perturbative expansion if we work with ν s = 1/r or with ν s = 2/r. Somewhat, working with ν s = 1/r gives the right result from the beginning, and adding more terms of the perturbative expansion makes the result to oscillate around the central value. On the other hand, working with ν s = 2/r, the LL result is quite off the expected result, but then adding higher order terms of the perturbative expansion makes the prediction to converge to the same result we have obtained with ν s = 1/r. The convergence is perfect within the statistical errors, and also irrespectively of working with ν us = C A α(ν s )/(2r) or with ν us = 1 GeV. Still, the convergence pattern is not equal in these two cases for ν s = 2/r. Not taking an optimal ν s (∼ 1/r) makes the determinations with (ν s , ν us ) = (2/r, C A α(ν s )/(2r)) or with (ν s , ν us ) = (2/r, 1 GeV) to be significantly different at NNNLL. Nevertheless, this difference is nicely elliminated after the inclusion of the terminants associated to the u = 3/2 renormalon. Even more so, the inclusion of JHEP09(2020)016 the terminants associated to the u = 3/2 renormalon is also fundamental to get agreement of these fits with the fits with ν s = 1/r. A similar discussion holds for the case with ν s = 1/( √ 2r), though the overall behavior is better: the LL result is closer to the value obtained with ν s = 1/r, and the difference between the determinations with (ν s , ν us ) = (1/( √ 2r), C A α(ν s )/(2r)) or with (ν s , ν us ) = (1/( √ 2r), 1 GeV) at NNNLL is small. Finally, for (ν s , ν us ) = (1/( √ 2r), 1 GeV) and (ν s , ν us ) = (1/( √ 2r), C A α(ν s )/(2r)), we get Λ = 335 MeV respectively with the NNNLL hyp theoretical expression. We then conclude that in the range ν s ∈ (1/( √ 2r, 2/r) the result is stable at the 2 MeV level if we fix ν us =1 GeV. The spread of the result is slightly larger, at around the 3 MeV level, if we set ν us = C A α(ν s )/(2r) instead. This very tiny increase can be interpreted by the fact that ν us is not completely constant as we change x s , since ν us = C A α(x s /r)/(2r) is not exactly equal to ν us = C A α(1/r)/(2r), the value we use in our reference fit. In the next item, we will study the dependence of our fits with more extensive variations of ν us . For ν s , we conclude that, with the present level of precision reached by the theoretical expression, the dependence on ν s of the fit, of order ∼ 2 MeV, can be neglected compared with other uncertainties.
In the whole parameter range we have studied, the χ 2 red is reasonable. Therefore, all fits are equally good in this respect. The only exception is the NNNLL prediction for (ν s , ν us ) = (2/r, C A α(ν s )/(2r)), which has a χ 2 red 1.9. We find then significant that it moves away from the convergent pattern that is observed in the other fits. It is also then significant that the inclusion of the terminant brings agreement with the other fits and lowers the χ 2 red down to χ 2 red = 0.42, much below 1.
Dependence on ν us . We now test the sensitivity of the fit on ν us . In order to keep the hierarchy of scales between the soft and ultrasoft scale, we have varied them in a correlated way as a function of a single parameter x: The range we take for x is x ∈ [1/ √ 2, 2], similarly as we did in the previous item. For the ultrasoft scale, we do so because, otherwise, we reach very low values for ν us that make α(ν us ) to blow up. This happens, for instance, if we take (ν s , ν us ) = (1/(2r), C A α(ν s )/(4r)). In this case, for the data Set I, the fit significantly deteriorates with a χ 2 = 338 MeV obtained with (ν s , ν us ) = (1/r, C A α(ν s )/(2r)). This is the same difference (with the same sign) as obtained with (ν s , ν us ) = (2/r, 2C A α(ν s )/(2r)). In this respect the fit with (ν s , ν us ) = (1/r, C A α(ν s )/(2r)) can be considered a (close to the) minimum within the families of fits with (ν s , ν us ) = . We show our results in figure 5. The agreement is very good. The difference is of order 9 MeV between the x = 2 and x = 1 fits and also between the x = 1 and x = 1/ √ 2 fits. The convergence is already reached at the NNNLL level. We observe that, by correlating the soft and ultrasoft scale as done in eq. (5.4), the convergence is already reached at the NNNLL level, and the contribution of the terminant associated to the u = 3/2 is very small. Note that this was not so when we took ν us = C A α(x/r)/(2r) or ν us = 1 GeV. For illustration we show again the fits with ν us = C A α(x/r)/(2r) and with ν us = 1 GeV in figure 5. In that case the terminant contribution is crucial to get agreement between fits with different values of x s . This reflects that the terminant plays a crucial role to diminish the dependence in ν us and to get convergence to the same value irrespectively of how we correlate the soft with the ultrasoft scale. We take the largest difference between the different possibilities we have considered (∼ 9 MeV) as an estimate of higher order effects of perturbation theory. Notice that the cancellation of the ultrasoft scale dependence comes from several places. On the one hand we have the perturbative contribution of d dr δE us (the magnitude of this contribution is small), we have the contribution from the derivative of δV RG , and finally the contribution from the terminant associated to the perturbative series of d dr δE us . Note that we only include one term of the pertubative expansion in powers of α(ν us ) in d dr δE us . One may be worried then that N P = 0 is too low for the incorporation JHEP09(2020)016 of the terminant associated to the ultrasoft energy. Nevertheless, we are working at rather low scales. We will indeed see in section 5.1 for the direct comparison of the static potential that the asymptotic behavior can easily set in at basically the lowest order.
Dependence on r ref . The fit should be independent of r ref . In practice, however, the result may depend on the value of r ref used, since the range where the logarithms of r are summed is different. This error also measures the fact that the data points have some error. For the data set I, we find the largest difference between fits with different r ref to be of order 8 MeV. For the other data sets the spread is slightly smaller, except for the data set IV, which is slightly larger (∼ 9 MeV) after considering the most extreme difference. This one is obtained with the largest r ref we have in our data set, which, on the other hand, produces a rather large χ 2 red : χ 2 red 4.8.
Dependence on Z F 3 . We have also studied the dependence of our central value on Z F 3 . We find it to be very small compared with other uncertainties, since the contribution associated to Ω F 3 is small for our central value determinations. The variation does not change the last digit. Therefore, we will omit it for the final error budget. It is worth mentioning though that for other values of ν s and ν us the terminant is important and, when so, it is a crucial element to get agreement with our central value.
Estimate of higher order contributions. For the error analysis, we need to determine the error associated to our lack of knowledge of the complete perturbative series. We have several ways to estimate this error. We have studied the error produced by the variation of ν s and find it to be very small, of the order of 2 MeV. We have also studied the error produced by the variation of ν us and find it to be of around 9 MeV for the data set I. As an alternative way to estimate the error, we considered the difference between the NNNLL and NNNLL hyp , i.e. adding or subtracting the terminant. This produces a very small shift. Alternatively, we have also performed fits changing the order at which we start including the terminant in the force from three to two. We indeed find the variation to be small: ∼ 6 MeV. The fits have been performed using the running of α with 4 loop accuracy [48], as it is the analogous accuracy to the perturbative expansion of the static energy. We have also made the fit including the running of α with 5 loop accuracy [49], and find a 3 MeV difference with our central value. To consider more conservative estimates of the error, we have also looked at the difference between NNLL and NNNLL fits. For the data set I, we obtain similar numbers, marginally larger, than from the variation in ν us : ∼ 10 MeV. We take the largest of all these possibilities. We believe this yields a conservative error estimate for the higher order contributions.
Final numbers. Out of this analysis, we proceed to give our prediction, for which we use the data set I. It reads The central value is taken from the fit of the NNNLL hyp theory expression with (ν s , ν us ) = (1/r, C A α(ν s )/(2r)) to the data Set I. The first error is the statistical error of the fit. The

JHEP09(2020)016
second one is the one associated to higher order corrections. We estimate it by taking the biggest number among the different estimates for higher order corrections we have discussed above, which corresponds to the difference between the NNNLL and NNLL number. We finally consider doing the fit with different r ref . We take the largest difference. This error is a mixture of two sources: on the one hand it is partially related to our lack of knowledge of higher order logarithms, and, on the other, on the error of the lattice data point. Still, we will treat it as an additional source of error. We then combine all errors in quadrature and obtain Λ This number can be compared with other determinations of the strong coupling at around these low energies. One can, for instance, compare with determinations using the heavy quarkonium spectrum [50,51]. Those also have as a fundamental input the static potential but, at present, they suffer from larger errors than those presented here. In this respect, applying hyperasymptotic expansions to these analyses may improve the accuracy of such determinations.
Out of these numbers we can also determine α (n f =5) (M z ). We follow the preferred method advocated in [52], which has built in the error from decoupling and truncation when going from the scales we have made the fit up to the M z mass. We obtain α (n f =5) (M z ) = 0.1181(8) Λ MS (4) Mτ →Mz = 0.1181 (9) . (5.8) The first error is the error associated to the error of our determination of Λ (n f =3) MS , and the second to the transformation of this number to α (n f =5) (M z ) as described in [52]. In the last equality we have combined the errors in quadrature. Our number is perfectly consistent with the world average number [53], or with the lattice final FLAG average value [54], and with a very competitive error.
As we have mentioned above, our prediction has been obtained using the data set I, which is the one less sensitive to long distances. Nevertheless, we have performed similar error analyses for the other data sets. We find Set II Λ Notice that all the central values obtained with the different data sets are within one sigma of our preferred value. The data sets II, III, IV have smaller statistical errors, but larger errors associated to higher order in perturbation theory effects, as they suffer from a larger difference between the NNLL and NNNLL result. at LL, NLL, NNLL, NNNLL and NNNLL hyp using the data set I with (ν s , ν us ) = (1/r, C A α(ν s )/(2r)) (continuous black line). We also give the determination of Λ (n f =3) MS at LO(LL), NLO(NLL), NNLO, NNNLO and NNNLO hyp using the data set I with ν s = ν us = 2/r (dashed green line), ν s = ν us = 1/r (continuous blue line) and ν s = ν us = 1/( √ 2r) (dotted red line). The error displayed is only the statistical error of the fits. We also show the error band generated by our prediction eq. (5.6).
Comparison with fixed order computations. Fixed order computations can be obtained from the RG improved ones by setting ν s = ν us . Therefore, this approximation does not incorporate the resummation of large ultrasoft logarithms. This effect can be important. We show the results of the fixed order computation and the comparison with the RG improved result in figure 6. Let us first remind that the first two orders are equal, i.e.: LL=LO and NLO=NLL, as there are no ultrasoft logarithms to resum. The difference shows up at higher orders. For the same value of ν s , and for order NNLO and NNNLO versus NNLL and NNNLL, the fits at fixed order give significantly lower values than those that perform the resummation of logarithms. On top of that, the incorporation of the u = 3/2 terminants does not improve the convergence of the determination, unlike when resuming the large ultrasoft logarithms, where we see a very nice convergence pattern. This shows that the resummation of large logarithms appears to be compulsory to find convergence, and to cancel the scale dependence that we have in the terminants. The magnitude of the incorporation of the terminants is larger for larger ν s . This may say that using N = 0 for the ultrasoft contribution for a scale as large as ν us = 2/r could be a bad approximation. In this respect notice that as we lower ν us , (see the fits with ν us = 1/r, and, particularly, with ν us = 1/( √ 2r) in figure 6), the convergence of the fixed order computation significantly improves.
What about if ν s = constant? In principle, the optimal way to resum the large logarithms is to scale ν s with 1/r and ν us with α(ν s )/r. In practice, the range of scales we JHEP09(2020)016 have is not that large. We then consider fits with fixed ν s and ν us . We choose (ν s , ν us ) = (1/r ref , 1 GeV). For the data sets I, II, III and IV, the NNNLL hyp fits yield Λ (n f =3) MS = (339, 342, 344, 346) MeV respectively. Note that, for the data sets I and II, the result is identical (difference is indeed below 1 MeV and only gets to 1 MeV after rounding) to the central values obtained before and displayed in eq. (5.5) and eq. (5.9) respectively. For the data set III, the difference is 1 MeV, and for the data set IV, the difference is slightly more significant: 3 MeV. This agreement is very rewarding, since fits at low orders in the hyperasymptotic approximation show large differences with the analogous fits using the default scales: (ν s , ν us ) = (1/r, C A α(ν s )/(2r)). For the data set I, we show the values of Λ 3.9, which then goes below 1 as we increase the accuracy) but then steadily converge to the central value, such that, as we said, the difference for the NNNLL hyp prediction is below 1 MeV. The fixed order fits, which are also displayed in figure 7, show the same kind of behavior to the one discussed in the previous item.
Nonperturbative corrections. In all the determinations above, we have assumed that the ultrasoft scale is in the perturbative regime. In this situation, nonperturbative effects are parametrically suppressed compared with the precision obtained with our hyperasymptotic approximation. This assumption is safer if we take the points with higher energy. For the points of the data Set I, ν us moves in the range ν us = C A α(νs) 2r ∈ (1.06, 0.86) GeV, for which we consider safe to use perturbation theory at the ultrasoft scale.
If the ultrasoft scale is in the nonperturbative regime, we can say little from first principles about d dr δE us . To make an estimate, we consider the data set IV after subtracting the points of the data set I (those at smallest distances that we used in the previous section for the determination of Λ MS in the purely perturbative regime). As a test, we assume that for this set of data the ultrasoft scale is in the nonperturbative regime. To simplify the parameterization of these nonperturbative effects, we assume that we are in the regime where 1/r Λ QCD α(1/r)/r. In this situation, δE PV us = k PV Λ 3 MS r 2 (where k PV is a nonperturbative dimensionless constant), instead of being equal to eq. (2.12) plus the terminant (i.e., the − 1 r 2 Ω F 3 (ν us ) is not included in the fit, unlike in the pure perturbative case, as such contribution is inside the nonperturbative term). We first want to see how sensitive the determination of Λ MS would be to considering the ultrasoft scale to be in the nonperturbative regime, which implies that we also have to fit δE PV us . For such fit, we obtain Λ MS = 356(3) MeV (the error is only the statistical error of the fit) with a χ 2 red = 0.55 (in the fit we have fixed ν us = 1 GeV in V RG PV ). Notice that this number for Λ MS is consistent with the value obtained from the pure perturbative fit (only a little bit more than one sigma away of eq. (5.6)). For the value of k PV we obtain k PV = −0.82 (7) . (5.12) In principle, we did not care much about k PV . Nevertheless, this value of k PV rings a bell.

JHEP09(2020)016
In the perturbative regime we have that where we set N = 0. At low scales this expression is dominated by the terminant, which, we remind, has the following form 14) The dependence on ν us is mild and, effectively, the terminant scales as  4). Therefore, what the nonperturbative fit does is to effectively fit the terminant assuming that the O(α(ν us )) term of δE PV us is subdominant. Notice that eq. (5.12) is, within one statistical standard deviation, the value predicted by perturbation theory, eq. (5.16). We take this as a very strong confirmation that our weak coupling analysis is safe, and that, indeed, one can apply perturbation theory to scales as small as 1/r ∼ 1 GeV. Finally, to confirm this picture, we do the fit over the complete data set IV assuming δE PV us = k PV Λ 3 MS r 2 . The results barely change: we obtain Λ MS = 355(3) with also χ 2 red = 0.55 and k PV = −0.8. Overall, we can even take this analysis as a strong indication that the data has enough precision to be sensitive (and, to some extent fit, albeit with large errors) the value of Z V 3 . This discussion also explains why the nonperturbative fit also has a small χ 2 red , as it loosely corresponds to the perturbative expression but letting the normalization of the terminant to be a free parameter of the fit.
Comparison with earlier work. Determinations of Λ (n f =3) MS using lattice data of the static energy have been obtained in the past. Some recent determinations are those of [28,29]. They compare with a different data set including lattice data at longer distances. They work directly with the potential. The precision of the theoretical expression is NNNLO in our counting. No resummation of ultrasoft logarithms nor the incorporation of the terminants is considered, but they use an alternative method for dealing with the renormalons. In the large β 0 , it is possible to see what is the precision that corresponds to in the hyperasymptotic approximation but not beyond the large β 0 . The ultrasoft scale is assumed to be in the nonperturbative situation. Therefore, the comparison should better be done with the number we have just obtained in the previous item. If we set JHEP09(2020)016 δE PV us = k PV Λ 3 MS r 2 and fix ν s = ν us = 1/r (i.e. we work with NNNLO hyp precision, we obtain Λ (n f =3) MS = 305(2) MeV where we only put the statistical error. This number is smaller than the number obtained in [28,29].
Closer to our analysis are [4,9]. In particular from the last reference we borrow the lattice data. In these references, they use the force as the starting point and later integrate it to recover the potential, as we have done above. Their central values are obtained by fitting to the NNNLO result after adding the NNLL ultrasoft contribution. Therefore, they mix different orders according to our counting and do not include the complete NNNLL result. The number they obtain is smaller than ours. In this respect, note that our numbers with analogous NNNLO precision are also smaller.

Direct fit from the static potential
We now present an alternative determination of Λ (n f =3) MS to the one obtained in the previous section. As in the previous section, we will mainly work with the data set I, but, in this section, we fit the lattice data to eq. (5.1) (we also consider here energy differences) using eq. (2.30). In this equation we will set N P = 1 by default. This also means that 3N P = N max = 3 and we reach the next renormalon. Nevertheless, we also explore the impact of choosing different values of N P . The counting of the hyperasymptotic expansion will then be the following: the static potential at tree level corresponds to the LO (which is equal to the LL) approximation. In the hyperasymptotic counting (D, N ) it corresponds to (0,0) precision. The static potential at one-loop corresponds to the NLO (which is equal to the NLL) approximation. In the hyperasymptotic counting, it corresponds to (0,1) precision. We then add the terminant associated to the u = 1/2 renormalon. We name such approximation NLO/NLL hyp1 . In the hyperasymptotic counting, it corresponds to (1,0) precision. We then add the O(α 3 ) terms in eq. (2.30). We name such approximation NNLL hyp1 or NNLO hyp1 if the O(α 3 ) contributions of δV RG are added or not, respectively. In the hyperasymptotic counting, it corresponds to (1,1) precision. We then add the O(α 4 ) terms in eq. (2.30). We name such approximation NNNLL hyp1 or NNNLO hyp1 if the O(α 4 ) contributions of δV RG are added or not, respectively. In the hyperasymptotic counting, it corresponds to (1,2) precision. Finally, we add the terminants associated to the u = 3/2 renormalon of V and δE us . Note that each terminant depends on the order we truncate each perturbative series and the scale of α in each respective perturbative expansion. We will name such approximation NNNLL hyp2 or NNNLO hyp2 if the O(α 4 ) contributions of δV RG are added or not, respectively. In the hyperasymptotic counting, it corresponds to (3,0) precision.
We will use the fits performed in this section to reassure the results obtained in the previous section. A priori one would expect the error of this determination to be larger due to the error associated to the first renormalon. We can avoid this error completely if we set ν s = 1/r ref . We then first perform a fit setting ν s = ν us = 1/r ref , so we avoid any r dependence in the renormalization scale. Since ν s = 1/r ref , the renormalon associated to u = 1/2 exactly cancels, since the scale is the same and the perturbative series is truncated at the same order. Therefore, the associated terminants are set to zero. Thus, JHEP09(2020)016 The error displayed is only the statistical error of the fits. We also show the error band generated by our prediction in eq. (5.6). the computation is equivalent to standard perturbation theory, and, in the counting above, we have NLO=NLO hyp1 , NNLO=NNLO hyp1 , NNNLO=NNNLO hyp1 . As we said, we also put ν s = ν us . Nevertheless, in this case, the renormalon associated to u = 3/2 of the static potential and of δE us do not cancel each other. The reason is that the perturbative expansion of the static potential is truncated at N = 3, whereas the perturbative expansion of δE us is truncated at N = 0. Therefore, the associated terminants do not cancel each other. The results are identical to those found in the previous section when setting ν s and ν us constant. We show the results in figure 7. We then take ν us = 1 GeV. This still avoids any r dependence in the renormalization scale and still NLL=NLL hyp1 , but the following orders are different, but only by a constant that cancels in the energy difference, except for the NNNLL hyp2 . Again, the results are identical to those found in the previous section using F when setting ν s and ν us constant. Therefore, we reach to the same conclusions: we saw in the previous section using F that working with ν s = 1/r ref and ν us = 1 GeV produced very similar fits to those with ν s = 1/r and ν us = C A α(νs) 2r GeV at high orders in the hyperasymptotic expansion, with a difference less than 1 MeV. As the results are identical, we observe the same behavior here: using directly eq. (2.30) at NNNLL hyp2 with (ν s , ν us ) = (1/r ref , 1 GeV) yields the same result (with a difference smaller than 1 MeV) than a fit using eq. (5.2) at NNNLL hyp with (ν s , ν us ) = 1/r, C A α(νs)

2r
. We show the results in figure 7. It is also very rewarding to see the stability of this result to changing N P . This far we have set N P = 1, but even if we set N P = 3 (such that we do not include the JHEP09(2020)016 We have investigated the origin of the problem. The fact that it only shows up when we introduce r dependence in ν s and ν us induces to think that it has to do with renormalons, similarly to the discussion one can find in [23] for the leading renormalon. The leading renormalon is under control. Therefore, the issues we face should have to do with subleading renormalons, maybe also with the renormalons encoded in ∆V . Here we do not have a clear explanation. We postpone a detailed analysis to future work. What we have been able to do is to identify where the effect seems to be hidden, and it is in δV RG . If instead of using eq. (2.5), we obtain δV RG from the perturbative expansion of its derivative, eq. (3.7), which we then integrate, we find that most, if not all, of the difference cancels. We show this effect in figure 8. We will not perform a full-fledged error analysis in this case, though, as we do not have a clear understanding of how the subleading renormalons are showing up when ν s = 1/r and ν us = C A α(ν s )/(2r). Finally, even though it is not discussed in this paper, notice that, when working with ν s = 1/r, the error associated to Z V 1 shows up, and it can be potentially large.

Conclusion
In this paper, we have obtained hyperasymptotic approximations for the static energy and for the force with a precision of (D, N ) = (3, 0). Our expressions also implement the resummation of large logarithms to NNNLL precision. We have used these expressions to obtain very precise determinations of Λ MS and α (n f =5) (M z ). We have mainly used the force as the starting point of our theoretical analyses. Our final result reads The resummation of logarithms and the introduction of the terminants associated to the u = 3/2 renormalon are essential to get a very well convergent series. This, together with precise data at short distances, allows us to get accurate values for Λ MS . The lack of any of these novel elements significantly deteriorates the convergence, and consequently, the accuracy of the prediction. Our fits are based on the shortest available data from [5,9] that do not suffer from lattice artifacts. This means that we have restricted the fit to r smaller than 0.353 GeV −1 . To make the weak coupling analysis as reliable as possible, the smallest soft scale we take is 2 GeV. This corresponds to an ultrasoft scale of 0.86 GeV. The fit shows no signal of needing extra nonperturbative correction. This is so even if we relax the infrared cutoff of the soft scale down to 1 GeV. The fits are still perfectly ok. In this respect, we do not have any indication from the fit of the need of nonperturbative corrections. One strong, very nontrivial, validation that we can use perturbation theory in our fit actually comes from fits assuming that the ultrasoft scale is in the nonperturbative regime. We do fits assuming that we are in the situation where Λ QCD α(ν s )/r. In this energy regime δE PV us = k PV Λ 3 MS r 2 and, besides, Λ MS , we also have to fit k PV . Within one sigma of the statistical error of the fit, the number obtained for k PV is the same as the one predicted by the terminant of the weak coupling expression of δE PV us .

JHEP09(2020)016
The largest source of error comes from unknown higher order corrections in perturbation theory. The statistical errors of the fit are small, though the dependence on r ref , which is a mixture of lattice and theory error is large. Increasing the number of points of the data set gives a very mild tendency to increase the value till stabilizing at 343 MeV, very well inside the error we give.
Whereas the cancellation of the leading renormalons is under control, the situation is not that clear for subleading renormalons. δE us depends on ∆V , which has a renormalon located at u = 1/2. The specific way this renormalon cancels is something that should be investigated. This is what has stopped us from using the perturbative NLO expression for δE us . In this respect, the analysis of [55] could be of help. In this reference, linear power-like divergences that appear in the coefficients of the perturbative expansion due to renormalons located at u = 1/2 become logarithmic divergences that are identified in dimensional regularization. In that specific example, one could see the cancellation between different terms. Likely related with this discussion, there is another issue that has to be investigated: the existence of renormalons in δV RG . We observe a rather different behavior if we first derive it and afterwards integrate it or if we directly use eq. (2.5). This difference only appears if ν s and ν us are made to be explicitly r dependent. This is consistent with the existence of a renormalon in this object. Indeed, it resembles the situation faced by early studies of the leading renormalon [23]. In that case, depending on how one does the expansion, in powers of α(1/r) or in powers of α(ν s = constant), the perturbative series was also convergent or not. Working with the force, we observe that the cancellation of renormalons is incorporated from the start, irrespectively of working with ν s =constant or not. In the case of working with the potential, we are, at present, only safe if working with ν s and ν us being r-independent. These issues will be investigated in future work.

JHEP09(2020)016
B F RG (r) with ν s = x s /r and ν us = x us C A α(νs) 2r We have (we only discuss the pure perturbative terms) F RG (r) = F PV (ν us = ν s ) + d dr δV s,RG (r; ν s , ν us ) + d dr δE PV us (r; ν us ) , (B.1) where the coefficients of F PV (ν us = ν s ) are Note that these coefficients are equal to those in eq. It is remarkable that the differences among the different terms cancel each other such that the expression for F RG (r) we obtain here is equal to the one used in section 3.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.