Superasymptotic and Hyperasymptotic approximation to the Operator Product Expansion

Given an observable and its operator product expansion (OPE), we present expressions that carefully disentangle truncated sums of the perturbative series in powers of $\alpha$ from the non-perturbative (NP) corrections. This splitting is done with NP power accuracy. Analytic control of the splitting is achieved and the organization of the different terms is done along an super/hyper-asymptotic expansion. As a test we apply the methods to the static potential in the large $\beta_0$. We see the superasymptotic and hyperasymptotic structure of the observable in full glory.


IV. Conclusion
neglect in what follows (as they yield smaller NP corrections than those we consider in this paper), such relation can be quantified using the operator product expansion (OPE) of the observable for large Q. The allowed operators determine the allowed corrections in powers of Λ QCD (up to logarithms), and, therefore, the large order behavior of the perturbative expansion, since the latter can be related with singularities in the Borel plane (located in the positive real axis), which mix with the NP corrections. To these singularities (and the associated asymptotic perturbative expansion) we generically refer to as infrared renormalons [1].
On a more general scenario one can consider more than one large scale: Q 1 Q 2 Λ QCD . Then the use of the OPE and the factorization between the different scales makes the perturbative expansions associated to each scale to be asymptotic. In some cases one has renormalon singularities associated to the scales Q 1 and Q 2 that cancel among themselves.
This is indeed the case for the leading renormalon singularity of the pole mass and the static potential, as first found in [2], and later in [3,4]. We name these renormalon singularities spurious.
So, in general, we want to: 1. Predict observables with e −A 2π 2. Avoid spurious renormalon problems.
In this paper we focus on 1), though our results will be relevant for 2) too.
Besides its intrinsic theoretical interest, the asymptotic behavior of perturbative expansions in QCD is starting to be seen in a series of observables, in particular, in heavy quark physics. In this case, in order to handle the renormalon problem associated to the pole mass, different threshold masses have been introduced [4][5][6][7][8][9]. Some of these threshold masses introduce (explicitly or implicitly) an scale ν f that acts as as an infrared cutoff. Such infrared cutoff kills the renormalon behavior of the perturbative series producing a convergent perturbative series and introducing a linear power-like dependence in ν f . In practice these threshold masses work quite well. The error associated to the fact that we have this linear cutoff is typically small (see, for instance, [10][11][12]). Still, it is not optimal conceptually 1 .
Other of these threshold masses use approximate expressions for the Borel transform of the pole mass that partially incorporate the renormalon singularities in the Borel plane. The inverse of the Borel transform (which we will name Borel sum or Borel integral in the following) is then ill defined. This requires using some prescription to regulate the Borel integral.
In this last case the perturbative series is typically abandoned and one directly works with the Borel integral expression. In this approach it is not quantified what is the error made by using (the unavoidably) approximated expressions for the Borel transform.
This discussion leads us to consider an alternative method that is also often used to tame the asymptotic behavior of the perturbative series: truncating the perturbative sum at the minimal term. In mathematical literature, such approximation is often named superasymptotic approximation of the original function (see [13]), which is a name we will also use in the following. This procedure has been used since long (see [14], or [15], for references), mainly in the context of solutions to one-dimensional differential equations. Nevertheless, in that context, renormalons do not show up, nor it does the issue of scheme/scale dependence.
In the context of four dimensional quantum gauge field theories, truncation of the perturbative sum in different formulations or using approximated expressions for the Borel integrals has also been considered since the early days of OPE/renormalon analyses to determine ob- 1 In the same way that there is nothing conceptually wrong in using cutoff regularization in perturbative computations, but regularizations that kill spurious power-like divergences, like dimensional regularization, and preserve more symmetries are much more convenient.
servables with NP accuracy (see for instance, [15][16][17][18][19][20][21][22]). However, it was not possible to make quantitative analyses beyond the large-β 0 approximation, since the existing perturbative series were only known to low orders. More recently, perturbative expansions have been obtained to high enough orders for some observables in the lattice scheme [23][24][25][26]. This has allowed to quantitatively use perturbative sums truncated at the minimal term and successfully determine the gluon condensate andΛ in the quenched approximation [27]. This success motivates trying to improve this approach, and to revisit with it observables already computed in the MS scheme, even if only few coefficients are known, since in the MS scheme (and in particular in heavy quark physics) renormalon dominance shows up at relative low orders.
Whereas, by construction, the superasymptotic approximation does not explicitly introduce the factorization scale ν f , the dependence on the renormalization scale ν remains to be assessed. Therefore, to push this method forward we need to get a quantitative understanding of the error on the truncation of the sum and of its remaining scheme and scale dependence. Similarly, the NP power corrections are potentially dependent on how the divergent perturbative series is regulated and on the renormalization scheme/scale used to define the strong coupling: α X (µ). A major point of this paper is to be able to control (in an analytic way) the dependence of the power corrections in this generalized scheme dependence. We will only then be able to add NP power corrections to the pertubative series in a systematic way, since the mixing between the perturbative series and the leading NP terms (or between the perturbative series associated to the scales Q 1 and Q 2 ) makes impossible to determine them independently. An unambiguous definition of the NP power corrections requires defining the perturbative series with power accuracy. Such combined expansion of perturbative series and NP terms will be called hyperasymptotic expansion as in [13]. Organizing the computation in this way allows us to precisely state the parametric accuracy of the result at each step.
The mixing between perturbative and NP effects may hinder estimating the real size of the NP effects. This happens when using threshold masses. In this case the problem is not severe. A more extreme example of this problem appears in lattice regularization. The gluon condensate is, up to a factor, the expectation value of the plaquette: where p 0 = 4π/3. For β = 3/(2πα) = 6.65 we have G 2 latt ∼ 3.3 × 10 4 r −4 0 , whereas the NP gluon condensate is ∼ 3.2 r −4 0 [27]. We see that the perturbative contribution overwhelms the NP contribution by orders of magnitude. Therefore, it is convenient to devise schemes where one has extracted as much information as possible from perturbation theory in such a way that the remaining NP object has a minimal mixing with perturbation theory. This scheme would provide a natural place to estimate the real size of the NP corrections without the distortion due to perturbative effects. We believe that in this scheme one could get a better understanding of the real structure (size) of the NP effects. This could be important, once the precision increases and to set an standard for the future.
It is also our aim to relate the hyperasymptotic expansion with the previously mentioned methods used to handle the pole mass renormalon. This will allow us to parametrically quantify the error those methods have, in particular those using approximate expressions for the Borel transform.
Finally, it is also worth mentioning that truncating the perturbative series at the minimal term can be motivated in the context of factorization of scales and effective field theories, where one wants to factor out the physics associated to Q from the physics associated to Λ QCD . The point is that in n-loop diagrams, new scales are effectively generated. These scales are proportional to Q, but are modulated by small factors ∼ e − n k , where k is an integer. The dominant contribution to the n-loop diagram is not then due to Q but to Qe − n k min , where k min is the smallest possible k for the process at hand. For the case of the pole mass k min = 1 and Q = m. In this case, for small n, we still have that me −n Λ QCD .
Nevertheless, for n ∼ 2π β 0 α , we have me −n ∼ Λ QCD . Doing perturbation theory for n > ∼ 2π β 0 α would simply mean treating me −n as much bigger than Λ QCD , which is incorrect.
The structure of the paper will be as follows. In Sec. II we discuss the general case. In Sec. III we discuss the QCD static potential in the large β 0 approximation. We use this quantity as a toy-model NP observable to test our methods.
Real QCD examples will be discussed in followup papers.

SCHEME DEPENDENCE
The generic form of the OPE of a dimensionless observable is the following: . On the other hand S(α X (Q)) can be computed as a Taylor expansion in powers of α X (Q). This series is assumed to be asymptotic. Up to some anomalous dimension, C O,d (α X (Q)) can also be computed as a Taylor expansion in powers of α X (Q) and the generated series is also assumed to be asymptotic (we will not explicitly ellaborate much on this fact though, since this leads us to consider subleading corrections in the OPE expansion, which can be handled in an analogous way). Then, the observable is often represented in the following way: where the dots stand for terms suppressed by higher powers of Λ QCD /Q, β 0 = 11 3 C A − 4 3 T F N f , γ is the anomalous dimension of the operator O d , and and so on. Obviously, the three equalities in Eq. (3) are symbolic representations of the observable, as the perturbative series are asymptotic. We need to define first the perturbative series with Λ QCD power accuracy. This definition of the perturbative series will, in turn, define unambiguously the NP power correction. In realistic cases, the only information that we will have of the OPE of the observable will be: 1. The exact knowledge of the coefficients p n up to n = N , where N 1 is large enough such that p n is well approximated by its asymptotic behavior 2 ; 2. The knowledge of the structure of the leading NP power corrections: values of d, γ and the very first few terms of p (d) n ; 3. The knowledge of the asymptotic behavior of p n (which relies on the previous item and demanding consistency to the OPE): 4. The knowledge of the µ dependence of p n and p '(d) n dictated by the renormalization group invariance. In realistic cases, only to some order.
Therefore, we will devise definition methods that only use this information. This naturally leads us to consider perturbative series truncated at the minimal term N * (or close by): Note that N * depends on µ and on the renormalization scheme X used to define the strong coupling constant: α X (µ).
Therefore, we define After truncating, S T (Q) depends on small variations of N around N * , on µ, and on the scheme X (in this paper we will consider perturbative expansions either in the lattice or in the MS scheme but the expressions are valid for general renormalization schemes). Overall, we generically label all the summation scheme dependence by T . The ambiguity (freedom) of the truncated perturbative series is "of-the-order" of the power correction. By this we β 0 α X (µ) ). We also assume that the truncated sum S T for N ∼ N * to be asymptotic to the full result in the following way where d is the dimension of the leading NP term of the OPE.
As the observable is summation scheme independent, the T -scheme dependence of S T (Q) should cancel with the scheme dependence of the NP power corrections. We would like to determine Eq. (9) with higher precision. As we have mentioned, one can get right the dimension d of the NP power correction by approaching N to N * . It is more complicated to fix the overall coefficient (and its structure in powers of α and ln α) that modulates the NP power correction. This will heavily depend on the freedom in truncating the perturbative series. It also needs some extra information in the relation between the perturbative series and the observable.
In order to quantify this difference, we first search for generalized summation schemes of the perturbative sum that are T -scheme independent, i.e. that they are independent of µ, X and N . The Borel integral of the Borel transform is a natural candidate. In our case the inverse of the Borel transform needs regularization, as it has singularities in the real axis at positive values of the integration variable. Here we take the principal value (PV) prescription of the perturbative expansion: where one takes the arithmetic average of the integral above and below the real axis and For values of t larger than the radius of convergence of this series, we take the analytic continuation of this function. For instance (this function will be useful later on), where u ≡ β 0 t/(4π). For this and related equations we collect a useful set of equalities in Appendix A. Now, our first task is to show that S PV is indeed T -scheme independent.
In the large β 0 approximation the PV Borel integral can be shown to be factorization-scale and scheme independent (actually in the large β 0 approximation both things are the same) [28]. Beyond the large β 0 approximation things are more complicated. Nevertheless, we can still show the factorization and scheme independence of S PV under some assumptions.
We first consider the renormalization scale dependence. We restrict the discussion to the inclusion of β 1 to the running of α. Then, renormalization scale independence gives the following relation between coefficients [28]: Using these relations we can deduce that This is a total derivative and vanishes. It is possible to include β 2 to the running of α.
Renormalization scale independence of the perturbative series now gives the following relation between the coefficients of the perturbative expansion Though much lengthier expressions show up, it is still possible to deduce that These are total derivatives. The behavior of g i (u) for small u is g i (u) ∼ u α with α > 0. For where h(u) does not grow exponentially. Therefore, g(0) = 0 and g(∞) = 0, proving the renormalization scale independence of S PV . The inclusion of higher order terms seems to produce also total derivatives that vanish. Note that our conclusion disagrees with [28].
We now turn to the scheme dependence. Given the perturbative series in a given scheme: we consider a general change of scheme: The independence on the coefficients d i of Eq. (18) produces the following relation Note also that Overall we get For simplification one could work in schemes that make higher order terms (the dots) to vanish. Then, in a similar computation to the one we did to get the scale dependence, we get the PV integral does not change under these variations, as we get total derivatives that vanish: Therefore, S PV is T -scheme independent.
We do not enter in this paper into global definitions of the Borel integral of the observable itself, which may not exist [1]. For the purposes of this paper, it is enough that we can define the Borel transform of the perturbative series and its Borel sum (with the PV prescription). We then assume that difference between the Borel sum regulated using the PV prescription and the complete NP result obtained from full QCD can be absorbed in the NP terms of the OPE. An analytic proof (of disproof) of that is tantamount to given a NP proof of the OPE in QCD, which is, at present, beyond reach. Since we assume that such generalized resummation scheme preserves the structure of the NP OPE, the difference with the observable has to exactly scale as the NP corrections of the OPE: where the last term refers to higher order terms in the OPE (d > d). K (PV) X is independent of µ and Q. We also demand K (PV) X to transform as Λ −d X under changes of scheme of the strong coupling, α X , i.e. the combination K PV X Λ d X is scheme independent. Indeed, since the structure of the NP OPE should be preserved, alternative generalized summation schemes should be different from the S PV by a term exactly proportional to the µ and scheme independent quantity Note also that the exponent γ and the O(α X (Q)) can be determined by RG analyses. In some cases RG analysis says that there is no O(α X (Q)) corrections (the Wilson coefficient is identically 1). This indeed would be the case of the B-meson mass.
S PV has the handicap, though, that it needs the full analytic structure of the Borel transform in the Borel plane, i.e. it requires the knowledge of the perturbative series to all orders. This can make them unpractical 3 . Remarkably enough, however, this problem can be bypassed by relating S PV with truncated versions of the perturbative series. This is the strategy we follow: devising truncated sums that we can relate with the PV result. This allows us to control the scheme dependence and error of using S T . Quite remarkable, this approach also allows us to quantify the error of using approximate expressions for S PV , since we do not know the complete perturbative series.
For fixed µ, the N → ∞ limit of S T (Q) diverges, since the perturbative series is divergent.
Therefore, if we want to keep µ finite, we have to keep N finite as well. Alternatively, if we want to take N → ∞, then we should as well send µ → ∞. Therefore, we explore two possibilities. One is to take µ ∼ Q in N ∼ N * , the other is to take µ → ∞ (correlated with N ∼ N * → ∞): 1) N and µ ∼ Q large but finite: 2) N → ∞ and µ → ∞ in a correlated way. We consider two options: where c > 0 but c is arbitrary otherwise. Note that in case 1), c can partially simulate changes on the scale and scheme of α X .
We will study case 1) and 2) in the following two subsections.
A. N large and µ ∼ Q Λ QCD . Eq. (27). Case 1) We first study option 1). Now the truncated sum reads (N P is defined in Eq. (27)) in realistic scenarios we will only have approximated evaluations of the leading singularity in the Borel plane.
We want to estimate what is the leading contribution to the difference between the PV sum and its truncated sum. This difference is dominated by the leading renormalon. Therefore, we focus on the contribution associated to it: where I is defined in Eq. (12). The finite sum stands for the contribution to S P (Q) associated to the leading renormalon. Ω is the terminant [14] of the asymptotic series when we truncate at α N +1 : where and ∆Ω admits the following integral (but not a Borel integral) representation With these definitions Ω has the desired asymptotic expansion: Even if Eq. (33) is not in a Borel integral form, this integral is amenable for a saddle approximation analysis (still, note also that we can evaluate it numerically exactly). We consider the integral where D b (x) is defined in Appendix A. Setting (to avoid considering non-integer values of N , for a given value of µ we will restrict to values of c that ensures that N P is integer) the integral H has the following expansion (this result is obtained by explicit computation and checked with an alternative computation using the recursion relations one can find in [14]) and Ω reads or in terms of Λ QCD , where Let us note that Eq. (38) also has a factor α γ (µ) besides the prefactor α(µ). In this paper we will only consider situations where γ = 0. To properly account for this factor one has to perform a resummation of ln(µ/Q) terms that effectively transform α γ (µ) into α γ (Q) in Eq. (38). For one example of a case with γ = 0 where this is done, see, for instance, [29].
In the large β 0 it is possible to write Ω in a Borel integral form. It reads After integration we obtain (η Obviously this result coincides with the full result when setting b 1 = β 1 = · · · = 0.
Subleading NP renormalons give subleading power corrections. A function with a finite radius of convergence in the α plane yields a Borel transform that is an analytic function in the whole complex u plane. Such function generates corrections smaller than any NP Since Ω gives the leading NP correction to S P V we can write Overall we obtain Note that with this method we do not expect a bad behavior when we take c → 0: The result is smooth, unlike what will happen with method 2B). Remarkable enough, this result quantifies the error of determinations of NP corrections obtained by truncating the sum at (or around) the minimal term, which is of O( α(µ)Λ d QCD ) irrespectively of the scale and scheme (in particular this applies to the analysis in [27]). We now can do better, as we now can compute these subleading terms that before went into the error. Therefore, we can increase the precision with which the genuine NP term can be determined.
If the precision of the computation is high enough one may consider going beyond the leading power accuracy and include the first correction to the above equations in the hyperasymptotic expansion. It would read where N P stands for the power in α where the perturbative series will mix with the subleading renormalon and Ω can be easily deduced from Eq. (37) adapting dimension and anomalous dimension to the next renormalon .
The truncated sum S P (Q; µ) depends on µ but not S PV . This has the important consequence that we can determine the µ dependence of S P (Q; µ) with Λ QCD power accuracy, and also to control the scheme dependence. We obtain We will typically take µ = kQ, where k is a constant of order 1 to avoid large factors.
Note also that the Taylor expansion in powers of α of the last term in Eq. (48) starts at n = N P +1. This effectively transform this term in a NP power correction. Moreover, the fact that the leading renormalon is subtracted from the perturbative series expansions further suppress this contribution. A naive estimate can be obtained by saturating the coefficients by the next renormalon. For the case of the static potential, the next renormalon is located at u = 3/2. This produces that the series roughly scales as which is obviously subleading, but still more important than the next NP correction. We will visualize the size of the different terms of the hyperasymptotic expansion in more detail in Sec. III B for the case of the static potential in the large β 0 approximation.
The correction associated to an analytic function in the whole complex Borel plane (of order α N ∼ e −#N ln(N ) ) is smaller than any NP correction (of order e −# N , where # is finite and bigger the further away the renormalon singularity is from the origin). Still, one can also worry about the role played by the logs generated in the perturbative computation: ln(µ/Q) . Assuming they are large, the leading contribution to the order α N is of O(α N ln N (µ/Q)). Since it is still 1/N ! suppressed compared with the renormalon contributions, it can be written as e −#N ln(N/(µ/Q)) . Obviously if k is made parametrically big it could jeopardize the hierarchy of the corrections that we have here. Therefore, we will always keep k parametrically of O(1).
We now illustrate the above general discussion using the particular case of the heavy quark mass (we neglect ultraviolet renormalons). We then have m PV =m[1 + S P (m; µ) + Ω m (µ) + · · · ], wherem ≡ m MS , we set d = 1, and also set the Wilson coefficient of the nonperturbative correction to 1 [30] in S P and Ω m .
We now compare our analysis with existing threshold masses. We focus on the RS mass [6] and relatives 4 . The RS mass is defined in the following way: where where one typically takes N = N max ≡ the maximal number of coefficients of the perturbative expansion that are known exactly (we assume that N max is not that high that we have to worry about subleading renormalon). In order to lessen the ν f scale dependence, the RS'≡ RS (1) was also defined: It is obvious that one could generalize to RS (n) where the subtraction starts at order α n+1 : Nevertheless, we can not increase n arbitrarily, otherwise the renormalon is not canceled.
Moreover, the value of n for which there is no cancellation of the renormalon will depend on µ. Therefore, when including higher orders one should do it with care once approaching to the minimal term. Another issue is the ν f dependence. To connect with the approach used in this paper we should take ν f = µ. Note that then r RS (n) In the original applications of the RS schemes this could be a problem, since the natural scale in the pole mass is different from the natural scale in the static potential 5 . To connect with the approach used in this paper, we control the scale dependence by fixing n = N = N P (µ).
This smoothly connect the RS schemes with the schemes where the series is truncated at the minimal term. One can then add Ω m and higher orders terms in the hyperasymptotic expansion of m PV .
We now consider the threshold mass named m BR , defined in [7] (see also [32]). The author directly works with the Borel transform and then regulate the Borel integral using the PV prescription. The complete expression of the Borel transform is not known. Therefore, in practice, an approximated expression is used that agrees with the known terms of the pole mass perturbative expansion till N = N max = 2 (the known coefficients at that time) and incorporates the leading singularity in the Borel plane. The author also makes a conformal mapping of the Borel transform. The µ dependence of m BR was usually fixed to µ = m, except in [33]. To make a quantitative comparison with our analysis, we leave aside the conformal mapping and make explicit the µ scale dependence in m BR . The key point then is the comparison of N (= 2) with N P . If N < N P (µ) there is power-like µ dependence that gets uncancelled with the contribution of S P . In other words m m Note that this produces an strong (linear) renormalization scale dependence (r (as) n ∼ µ) that is missed if one sets µ =m. This problem is potentially more severe in top physics (see for instance [34]), since one includes orders in perturbation theory beyond those presently known if the perturbative expansion is made with α(m t ). 5 If the scales are widely separated, this problem could be overcome using the resummation of logarithms of ν f , as first worked out in [31].
Overall, the only problematic situation would be if N < N P . For N ≥ N P , m BR and m PV are equal within the approximation used, and our analysis reorganizes the result within an hyperasyptotic expansion. This allows us to quantitatively control the µ dependence, and to parametrically state the error, of the result (for a given truncation) with NP power accuracy using an hyperasymptotic counting.
We can also connect our results with m MRS , defined in [9], in the following way (the expression of J can be found in Eq. (2.17) of [9]).
In this definition, µ has been fixed tom. By doing so we can not estimate the error associated to the µ dependence of m MRS . Therefore, we introduce it and generalize the definition of m MRS in the following way [we could indeed write a more general definition by putting a different scale for the renormalon term: m MRS (ν f ) = m RS (ν f ) + J (ν f ). This would still achieve renormalon cancellation]: which makes explicit the µ scale dependence of the definition. In principle one could think that, since it is related with RS mass, this would make a linear dependence in µ appear.
Remarkably enough this is not the case. We can relate this expression with the quantities defined above. Indeed the difference between m (N ) This quantity diverges in the large β 0 limit, which makes not possible to take the large MRS (µ) (alternative definitions were then proposed in [9]). The possibility to subtract this term from the PV regulated Borel integral was also considered in [32], though with a different (but related) motivation. In this respect, we note that subtracting this quantity from the PV result has been criticized in [35], on the basis of analytic properties of the observable. Nevertheless, this discussion is not directly relevant for us 6 , as adding or subtracting this term would just be equivalent to a change of resummation scheme that can be absorbed in the genuine NP power correction. Note though that this difference is parametrically bigger thanmΩ m , since the latter scales like O( √ αΛ QCD ). In any case, since the difference with the PV result is a scale/scheme independent quantity proportional to Λ QCD , the comparison with our analysis runs in complete parallel to the previous discussion of m (N ) BR with respect to N . Again, problems will show up if N < N P , but for N ≥ N P , m MRS is equal to m PV within the accuracy of the computation, except for Eq. (61). Therefore, it can be written in terms of a modified version of the hyperasymptotic expansion discussed in this section.
A more extensive discussion and a quantitative analysis for the case of the top, bottom and charm quark masses will be carried out in [36].
As promising as method 1) is, it is worth to explore alternatives that yield results that are explicitly N (and therefore µ) independent. They may also lead to a better analytic understanding of the observable. This can be achieved by taking µ and N going to infinity in a correlated way. The simplest possibility one may consider is taking the limit as in 2A) in Eq. (28).
The case 2A) was studied in the large β 0 limit in [37,38] for the case of the static 6 It would be if we were able to relate the PV Borel integral with a NP definition of the observable. potential (a more general case, including subleading corrections to the running of α, was also considered in [38]). It was observed that S T was logarithmically divergent in N and the proportionality coefficient found. Nevertheless, it was not possible to get a direct connection of this coefficient with the normalization of the leading renormalon in the Borel plane. This problem has been solved in [39], where it has been shown how to relate the coefficient of the ln N term with the normalization of the renormalon. This analysis has also been done for the Adler function. Unfortunately, the validity of these findings is restricted to the large β 0 approximation.
Beyond the large β 0 approximation only the static potential has been studied [37,38].
Remarkable enough the ln N (and an associated ln(ln rΛ QCD )) behavior survives, albeit with different coefficients. This may point to certain universality (beyond large β 0 ) of this result.
Unfortunately, it is not known now how to relate such coefficient with the normalization of the renormalon. This would be very useful for analyses beyond the large β 0 .
We will discuss all this in more detail in Sec. III C 1 where we study the static potential in the large β 0 approximation in this limit.
We have seen that S T was logarithmic divergent in N when taking the limit 2A). It was also not possible to connect S T with its Borel sum. We now consider the limit 2B). In this case one truncates before reaching the minimum, i.e. for N < N * = d2π β 0 α X (µ) . This will yield a finite result. The other point we address is the relation of S T in the limit 2B) with its Borel sum.
For some specific models of sign alternating perturbative series, it was soon realized that the N → ∞ limit of their associated truncated sums could be related with a modified version of the Borel integral [40,41], if such N → ∞ limit is performed in an specific way. For instance, it was shown that where with Later work generalized this result to more general series expansions, even to some that show a non-sign alternating series (but assuming that their Borel transform has a finite radius of convergence), and for arbitrary χ (as far as it satisfies some conditions). Their results can be summarized in the following equation: where and we require χ to be such that ∞ j=0 p j (τ 0 ) j! t j is analytic for |t| < 4π β 0 χ . Therefore, we can indeed sum the Borel series unambiguously inside the disc.
This was originally proven in [28,42] by brute force computation. It was also proven using a different method (integration in the complex plane) in [43] (in this last reference the O(1/ √ N ) corrections were also computed). In both cases the running of the strong coupling is restricted to follow the large β 0 approximation.
Whereas the above result applies to arbitrary perturbative series (with the qualifications mentioned above), the running of α is constrained to follow the large β 0 approximation. This is an important constraint if we want to consider the case of QCD, where the perturbative expansion of the beta function is not a monomial but has more terms. One can bypass this constraint if Eq. (2) in [43] is understood as a change of scheme instead of a change of a renormalization scale. It is also possible to generalize the derivation of [43] for a strong coupling with a general beta function. In this generalization new 1/N terms are generated.
Alternatively, one can slightly modify how the µ → ∞ is taken in Eq. (28). Instead of case 2B) one can take The difference with N A vanishes when µ → ∞. With this modified scaling it is possible to show that Eq (5) in [43] holds taking k = N χ with χ = d/(1 − c α(Q)). The derivation is then analogous to the derivation in [43]. Overall, we are then able to obtain (taking the µ → ∞ limit according to 2B) of Eq. (28)) beyond the large β 0 approximation, where 1 χ < d 2 . In particular, we will take 1/χ close to d/2, and parameterise it in the following way: where c > 0 (this is the reason we took c > 0 in Eq. (28)). The reason for the sign of c is that we have to approach to the closest singularity to the origin in the Borel plane from the left. Indeed, in [28,42], in the context of the large β 0 approximation, it was shown that in order the integral to be well defined one needed 1 χ < d 2 . It was also noticed that by taking the limit 1 χ → d 2 the correct exponent (of the NP power correction) is obtained, i.e. the difference is of the order of the leading NP term of the OPE. Nevertheless, one does not get the right prefactor. This was quantified in [43], where it was first shown that using Eq. (69), and expanding in α, the ambiguity is of the order of the higher order condensate with the right α dependence of the prefactor.
The leading renormalon (the singularity in the Borel plane closest to the origin) gives the main contribution to the difference between S P V and S A : This yields where Subleading corrections to the leading renormalon are of the form (1 + n > 0) This gives O(α 1+n ) corrections.
All subleading renormalons potentially contribute to the same order: This contribution is O(α 1+db−γ ) suppressed with respect the leading term. This is a problem if one wants to obtain subleading corrections to the leading NP term, as one would need to know the normalization coefficient of all subleading renormalons.
An issue observed in [43], in the context of the large β 0 approximation, was that when 1/χ → d/2, i.e. when the integrand approaches the singularity of the Borel transform, the truncated PV integral diverges, and it is not a good approximation of the PV integral (for instance see Figs. 2 and 3 in [43]). Therefore, it is better not to make the combination c α(Q) very small. We study this problem in the example we will consider in the following section.
This observation also makes that we can not use the results obtained in this section to the case 2A) obtained in the previous section, as it means setting χ = 2/d, i.e. exactly at the singularity in the Borel plane. (yet it would be very interesting a dedicated study to see if the analysis of this section can be generalized to the case χ = 2/d).

D. Strategy
In summary, we have two alternative expressions (Eqs. (46) and (71)) to determine S PV (Q) with Λ QCD power-like precision. Remarkably enough, we can achieve such precision even though we do not know the complete perturbative series expansion. The reason is that we can relate S PV (Q) with the truncated sum of the perturbative series for both methods. We also obtain an analytic expression for the leading power correction that accounts for the difference between the truncated sum and the PV result. One important feature of this result is that, in both cases, the leading power correction can be determined if the strength and structure of the leading singularity in the Borel plane is known. This result is also true beyond the large β 0 approximation. Such results are scheme independent.
There are important differences between both methods beyond the above general properties. The first one is that the method 2B) (the "µ → ∞ method") yields a finite NP correction in the limit Q → ∞. This is not so for the method 1) (the "µ = Q method").
For the latter, the leading NP correction gets multiplied by the small factor α(Q), which vanishes (albeit weakly) in the Q → ∞ limit. In principle, this makes the second method better. Nevertheless, one should also keep in mind that, in order to profit from this property, one needs to have physical data for as large as possible Q. Since in both cases the leading corrections are known analytically this could not make a practical difference. A numerical analysis can check which one is better. A more serious problem with the "µ → ∞ method" is that, in order to take the µ → ∞ limit, one needs the running of α with higher and higher precision. In the large β 0 limit, the running of α is known exactly, so this is not a problem, but it will be once we move beyond this approximation. One also needs higher and higher order coefficients of the perturbative expansion as one takes the µ → ∞ limit. Again in the large β 0 limit the coefficients can be generated to any arbitrary finite order 7 but not beyond the large β 0 limit. In the real case, the most we will have is the asymptotic behavior of the high order coefficients.
Another important issue is that with the "µ = Q method" we are potentially capable In general it is impossible to obtain closed results for the PV regulated perturbative sum on which to test the above results. This is only possible in the large β 0 approximation for a few cases. Here, we use one of them as a laboratory to check the methods we will apply to physical cases. The question here is to quantify the difference between the PV result (which we take as a "fake" NP data), and the truncated perturbative expansions (for large values of N ). Obviously such comparison is made in the short distance limit where the OPE 7 For the static potential this is indeed so, but even for the pole mass this is numerically demanding.
should apply. In Sec. III, we check our formulas (in the large β 0 approximation) for the case of the static potential. This example will allow us to quantify (in practice) when the complete result is well approximated by Eqs. (46) and (71). In particular, we try to answer the following questions: How large Q has to be in both cases 8 , how large µ has to be for Eq. (71) to hold. We also study the dependence of the answer to the scale/scheme used for the strong coupling (we use lattice and MS scheme).
The method that leads to Eq. (71) requires µ → ∞. Formally, this means that we need all the coefficients p n . As in realistic cases we do not have this information, we check the dependence on approximating the exact perturbative coefficients (starting at different orders) by their asymptotic expansions in the large β 0 approximation. In this case we will be able to see the error introduced by considering different orders from which one approximates the coefficients by the asymptotic behavior. What we will not be able to test in the large β 0 approximation is the dependence on the higher order coefficients of the beta function, which are needed for Eq. (71) (since we need to run α(µ) to µ = ∞). This is relegated to subsequent work.
Note that all the scheme dependence (in the broad sense: T ={N , X, µ}) has disappeared up to terms beyond the accuracy we achieve. We also obtain expressions for the difference between different truncation schemes.
Overall, we express the observable in the following two alternative ways up to exponentially suppressed terms. Note that Ω scales like O( α X (Q) Both methods have Λ QCD power accuracy but with the method 1) we have enough theoretical precision to determine the subleading O(α X ) corrections or even subleading terms in the OPE (hyperasymptotic) expansion (provided the "experimental" data is precise enough).

III. THE STATIC POTENTIAL IN THE LARGE β 0 APPROXIMATION
The large β 0 approximation cannot be obtained from a well defined limit of the parameters of QCD. Still, it is useful to test techniques that can be used beyond the large β 0 approximation in a place where we know the exact solution. In this respect the static potential is an ideal object, since we have a lot of analytic control for it.
The QCD static potential is written in terms of its Fourier transform as This equation defines α v (q) in the V-scheme. In the large-β 0 approximation, we know the behavior of α v (q) as a series in powers of α X ≡ α X (µ) where L = β 0 α X 2π ln( µe −c X /2 q ). If X = MS then c MS = −5/3 (in the large β 0 approximation).
If X = V then c V = 0. If X = latt, we take the n f = 0 number for a Wilson action: c latt = −8.38807 [44], as we will only use this scheme for checking the consistency between the results obtained with different schemes. We also defineΛ = Λ X e −c X /2 and ρ =Λr. Note thatΛ is scheme independent.
Eq. (77) is ill defined but not its Borel transform. It reads [45] B which is a meromorphic function in the u complex plane.
We then define (where the single poles of the Borel transform are regulated using the PV prescription) We can also regulate Eq. (77) via We have checked that the numerical determinations of both definitions give the same. We can then use this PV prescription as a NP definition of the observable, to which to test our methods and approximations. Note that this definition is indeed scheme independent. On the other hand the result is an oscillating function of r, which violates general properties of the static potential (energy) of two static sources in the fundamental representation [46].
These state that the potential should be concavus (we should also keep in mind that we are working in the large β 0 limit, which is not a well-defined limit of QCD).
We now consider the short distance limit (r → 0) of V PV (r). In other words, we analyze its OPE. First, we study how well we can approximate V PV (r) by its perturbative expansion at weak coupling. Thus, we approximate the potential by the truncated perturbative sum: For fixed µ, the N → ∞ limit of V N diverges since the perturbative expansion is asymptotic.
Therefore, we have to be careful in the definition used for the truncated sum. For such object, we use the two definitions discussed in Sec. II (with Q = 1/r). For both of them we will need the normalization of the leading renormalon. In the large β 0 it reads It agrees with the result from the pole mass [47] after using that the renormalon of the pole mass cancels with the renormalon of the static potential [2].
We will perform computations with n f = 0 and n f = 3. In the first case we will work in lattice units (aiming to compare with quenched lattice simulations) and use Λ MS (n f = 0) = 0.602r −1 0 ≈ 238 MeV [48]. In the large β 0 approximation (with n f = 0), this yields α(M τ ) ≈ 0.29. In the second case we take Λ MS (n f = 3) = 174 MeV. This last number we fix such that it gives a reasonable value at the τ mass in the large β 0 approximation: α(M τ ) ≈ 0.3 (see for instance [49]).
We then confront V PV with the results obtained with these methods.
Applying Eq. (48) to the static potential in the large β 0 approximation, the relation between V PV and V P reads where Ω V reads for this case with K (P ) and so on. Note that in the large β 0 we identically have Λ X = µe −2π/(β 0 α X (µ)) . This makes that K (P ) X,i . A similar expression applies to Ω V ∼ α X (µ)(rΛ QCD ) 3 . By incorporating the last two terms in Eq. (85) we are sensitive to the next renormalon.
Note that subleading renormalons give Λ QCD power-suppressed corrections. The further away the singularity in the Borel plane, the more suppressed the correction is. For the next-to-leading singularity we have   In the limit case 1), Eq. (27), we do not have direct analytic control in the relation between V PV and V P (unlike what will happen in Sec. III C when using the limit case 2), Eq. (28)).
Nevertheless, we can numerically compute both and check that their difference complies with the theoretical expectations. We can study (even if in the large β 0 approximation) up to which values of r the OPE is a good approximation of V PV . Remarkably enough we can actually check more than one term of the OPE (hyperasymptotic) expansion. We also explore the scheme dependence by performing the computation in the lattice and the MS scheme (actually in the large β 0 approximation this is equivalent to a change of scale). We will do these analyses for the cases with n f = 0 and n f = 3. The first in view of comparing with quenched lattice simulations, the second to simulate a more physical scenario, for which we can draw some conclusions that could be applied beyond the large-β 0 limit. In Figs. 1, 2, We do such computation in the lattice (Fig. 1) and the MS (Fig. 2). In Fig. 3 we compare the results in the lattice and MS scheme. We observe a very nice convergent patter in all cases down to surprisingly small scales. To visualize the dependence on c for each case, we show the band generated by the smallest positive and negative possible values of c that yields integer values for N P . The size of the band generated by the different values of c (the c dependence) decreases as we introduce more terms in the hyperasymptotic expansion. This is particularly so when including Ω V (Ω V ) to its associated sum.
Let us discuss the results in more detail. We first observe that the r dependence of V PV is basically eliminated in V PV − V P , as expected. This happens both in the lattice and MS scheme. The latter shows an stronger c dependence. This is to be expected, as in the MS, we truncate at smaller orders in N . This makes the truncation error bigger. Note that the lattice scheme can be understood (in the large β 0 approximation) as the MS scheme with a larger factorization scale. As we can see in the upper panel of Fig. 3, both schemes yield consistent predictions for V PV − V P . We can draw some interesting observations out of this analysis. For V PV − V P it is better to choose a larger factorization scale, if we have enough coefficients of the perturbative expansion. This is particularly so at large distances: We can still get sound results up to very large distances in the lattice scheme.
We now turn to V PV − V P − 1 r Ω V . Adding the new correction produces a better agreement with expectations (which we remind is to get zero). After the introduction of 1 r Ω V , the MS yields more accurate results than the lattice scheme. This can already be seen in the upper panel of Fig. 3, and in greater detail in the lower panel of Fig. 3.
r Ω V shows a dependence on 1/r, which is more pronounced in the lattice than in the MS scheme. As in the large β 0 the difference between both schemes is equivalent to a change of scale, this results points to that µ = 1/r in MS is close to the natural scale and minimize higher order corrections. Note that the lattice scheme computation is equivalent to the MS scheme choosing µ latt = µ MS e − c latt 2 e c MS 2 . This gives around a factor 30!!. Once )α n+1 is incorporated in the prediction most of the difference disappears and the lattice scheme is marginally better. Nevertheless, after introducing Ω V , the MS becomes marginally better again. In any case, the difference between schemes gets smaller and smaller as we go to higher orders in the hyperasymptotic expansion, in particular at short distances.
We also want to stress that this analysis opens the window to apply perturbation theory at rather large distances. Note that in the upper panel plots in Figs. 1, 2, and 3, we have gone to very large distances.
As some concluding remarks let us emphasize the following points. The truncated sum is more or less constant with relatively large uncertainties. This is to be expected, as the next correction in magnitude is Ω V which is approximately constant (mildly modulated by α(µ)). After introducing this term the error is much smaller and we can see more structure. In particular we are sensitive to 3N P n=N P +1 (V n − V (as) n )α n+1 . Here we find (at the level of precision we have now) a sizable difference between lattice and MS. This can be expected.
)α n+1 is the object we expect to be more sensitive to the scale.
In the lattice and MS scheme, we observe a very nice convergence pattern up to (surprisingly) rather large scales. The agreement with the theoretical prediction (which is zero) is perfect at short distances. The estimated error is also expected to be small. It would be interesting to see if this also happens beyond the large β 0 .
Another interesting observation is that truncated sums are better in the lattice scheme than in the MS scheme. Nevertheless, this could be missleading. The sums are truncated at the minimal term. Therefore, one needs more terms in the lattice scheme. If the number of terms is not an issue (which could be the case with dedicated numerical stochastic perturbation theory (NSPT) [50,51] computations in the lattice scheme) then the lattice scheme looks better. But as soon as Ω is introduced in the computation MS behaves better (at least in the large β 0 approximation).
We now turn to the n f = 3 case. We note that Λ QCD for the physical case (n f = 3) is smaller than for the n f = 0 case (if one sets the physical scale according to r −1 0 ≈ 400 MeV). On top of that the running is less important. All this points to that the convergence should be even better than in the n f = 0 case (and it was quite good already there). We show our results in Figs. 4, 5 and 6 (these are the analogous of Figs. 1, 2 and 3 but with n f = 3).
These plots confirm our expectations. Down to scales as low as 667 MeV we see no sign of breakdown of the hyperasymptotic expansion. This is so in both the lattice and the MS schemes. Note that the precision we get is extremely high as we go to small scales: Using truncation (c): )α n+1 (green), one gets V PV in both schemes with a precision well below 1 MeV at scales of the order of the mass of the bottom. Using truncation (d): r Ω V , the error is astonishingly small (see Fig. 8 for an extra zoom in this region). The rest of the discussion follows parallel the one for n f = 0.
In the above numerics, we have used the exact expression for Ω V and Ω V . In full QCD, we will not know the exact expression. Therefore, it makes sense to study how well the exact result is reproduced by the semiclassical expansion obtained in Eq. (37). We compare in Table   I and II for an illustrative set of values the exact result and the truncated semiclassical expansion. We observe that the exact result is very well saturated by the first terms of the expansion computed in Eq. (37). Truncating the expansion produces differences much smaller than the typical precision of the different terms of the hyperasymptotic expansion.
An alternative, very effective, presentation of the above results can be done by plotting the relative accuracy of the prediction at each order in α and at each order of the superasymptotic expansion. We note that we have one observable for each value of r. Therefore, for illustration, we make the comparison with the observable for r = 0.1 GeV −1 , and for the theoretical prediction we take the smallest positive value of c corresponding to lattice or MS. We show the results in Fig. 7. We stress that several terms of the hyperasymptotic to reach the same precision. One important lesson one may extrapolate from this exercise is that, for a fixed order computation, the smaller the renormalization scale µ, the better.
One can obtain much better precision for an equal number of perturbative coefficients. Another observation is that the minimal term determined numerically need not to coincide with the minimal term computed using n = N P (though it should not be much different). The    The potential advantage of this method is that we can obtain analytic results that are µ independent. We profit from earlier analyses in [37,38] adapted to our case. In all cases the q integrals will be done in the complex plane along similar lines as the computation done in those references.  We first truncate the sum of the α v coupling:

MS-Scheme (n f = 3) r in GeV
Following [37,38] we can isolate the N -dependence from the leading contribution to the n )α n+1 |, where in the last sum the two first renormalons are subtracted. Jumps correspond to the inclusion of Ω V and Ω V . Full points have been computed in the MS scheme and empty points in the lattice scheme. We work with n f = 3. potential at short distances: where arctan(x) is defined in the branch [0, π), and We then have that Note that this equality allows us to write V PV in the following way (v C = v 1 (ρ) − π ρ with the notation of [37]): In this explicit representation of V PV each term scales differently in powers of ρ: the ρ 0 term is set to zero (or incorporated in v C ), and each O(ρ 2n+1 ) term is encoded in terms. However, this splitting naturally leads to define a short distance coupling: This definition has nice properties. It is an smooth function ∀ r ∈ (0, ∞), with the right short distance limit: A detailed study of this quantity can be found in [38]. Note also that in this definition the whole O(Λ QCD ) correction has been included in 4C F β 0Λ v C (r). The other thing that one could study, since we have the analytic behavior, is the behavior of α SD beyond the regime where it was originally defined, i.e. at long distances. In this respect, it is interesting to notice that the long distance limit is exactly equal to the value obtained in [52], within the context of analytic perturbation theory analyses. Nevertheless, one could as well argue that all O(ρ 2n+1 ) terms are short distances and should be incorporated in α SD . If one does so, α SD does not have an smooth limit for ρ → 0 anymore. Finally, one could also study the β function of α SD .
It has some interest to compare Eq.  Fig. 6).
r Ω V (blue bands) in the lattice and MS scheme with n f = 3 (as drawn in Fig. 6). Note that in this last figure the vertical axis is in MeV and the precision is at the level of 10 −2 MeV!.
consistent with the estimated made in Eq. (50). Either way, the convergence is extremely good. The precision is much below the MeV.
In real life we will not have such complete analytic control and must rely on the methods discussed in Sec. II. Therefore, we now apply the limit 2A) and 2B) discussed in Eq. (28) to We now take . (98) The large N limit of v 2 yields up to terms that vanish when N → ∞. Note that the N → ∞ limit of v 2 (logarithmically) diverges. Note also that when ρ → 0 the integral term tends to zero. Thus, the ρ ∼ 0 limit The difference between the PV and the truncated series can be computed by complex variable integration following similar lines as in [37,38]. We find (for large N ) For large values of N and small values of r (care should be taken when taking the r → 0 limit) the above expression simplifies to For completeness, we have also obtained the ln(N ) behavior in a different way. We follow the method recently proposed in [39]. There, a summation integral relation was found for a general observable. We applied it to the case of the first IR renormalon of the potential and pole mass. The advantage of this new method is that the ln N term can be determined if the normalization of the leading renormalon in the Borel plane is known. It would be very interesting to try to generalize this result beyond the large β 0 approximation, as well as to extend the analysis to the ln ln( 1 ρ ) term. The fact that we have certain analytic control of the result allows us to address some issues. The first one is to make explicit that truncated sums around the minimal term do not guaranty, per se, that they are finite. In particular, one can see that V N is divergent in the N → ∞. Therefore, it would be wrong to assign V N to the leading term in the hyperasymptotic expansion of V PV . On the other hand, we have analytic control on the divergence, which is found to be logarithmic in N . 9 In principle, one can subtract this ln N divergence from V N (this is completely analogous to subtracting 1/ divergences in perturbative computations using dimensional regularization) to obtain the first term of the hyperasymptotic expansion. Nevertheless, the difference does not still scale like Λ QCD . Instead one has which, at short distances, scales as Λ QCD ln ln( 1 ρ ) (this behavior is also seen beyond the large β 0 approximation in the context of the static potential [38]). Therefore, to get the proper scaling in Λ QCD of the different terms of the hyperasymptotic expansion requires that the Λ QCD ln ln( 1 ρ ) should be identified and subtracted first from V N . One then has the freedom to subtract O(Λ QCD ) finite pieces, which can be absorved in the next term of the hyperasymptotic expansion.
We do not do a numerical analysis here, as the method cannot, at present, be generalised beyond the large β 0 approximation. 9 It is worth mentioning again that this ln N behavior also appears beyond the large β 0 approximation in the context of the static potential [38].

Case 2B)
We now take Under these conditions, we can take the N → ∞ limit (the result does not diverge in this limit). Adapting [38] derivation to our case we obtain Therefore, we define (using the relation Eq. (104)) 10 Note that this far, the expressions for v 1 and v 3 are valid ∀ r. It is also possible, and most relevant for us, to relate the truncated sum (in the limit µ → ∞) with V PV . We obtain Again this result is valid ∀ r. We now focus on the ρ → 0 limit. This will allows us to connect with the limit 2B) of Eq. (28). Nevertheless, this connection has to be done with care. One has to take the limit r → 0 and s → 2 in a correlated way, following the limit 2B) of Eq. (28). Therefore, we take Then, the previous expression reads 2 cos( π 2 (−1 + c α(1/r))) + ln(Λ r x ) sin( π 2 (−1 + c α(1/r))) ln 2 (Λ r x ) + π 2
We can now obtain the ρ → 0 limit: where, for x ∈ R, Nicely enough Eq. (110) agrees with the prediction of Eq. (71) applied to V PV .
For future reference, we are also interested in the next correction in powers of α(1/r) of Eq. (110). We obtain Note though that Eq. (71) can not predict the O(α) correction.
We have already emphasized that obtaining the ρ → 0 limit was delicate. Let us illustrate this. If we take the ρ → 0 limit with s fixed (but close to 2), such that s < 2 , we obtain If we rephrase this discussion in terms of the c behavior, what we have is that Eq. (110) is not obtained by taking the limit c → 0 before taking the limit r → 0 of Eq. (109).
Indeed, the limit c → 0 before taking the limit r → 0 produces Eq. (114), which does not correspond to the limit 2B) we are following in this paper. As we can see from the explicit computation, both limits yield NP power corrections with the right scaling (pointing out that there is not unique procedure to get/define the NP correction). Nevertheless, the overall coefficient is different, whereas Eq. (110) diverges logarithmically in c , Eq. (114) diverges like 1/c for small c . In this paper we stick to method 2B) as it allows us to go beyond the large β 0 approximation and to relate the normalization of the power correction with the normalization of the renormalon. X Λ X for n f = 3 in the lattice and MS scheme. For each case, we generate bands by computing V A with c = 1 and c = c min . We also compare with (b) r Ω V obtained with method 1) with the bands generated for Fig. 6. X Λ X at c = 1 and c = c min generating a band. We also explore the dependence on the scheme by comparing the results in the lattice and MS scheme. We stress again that in the large β 0 approximation lattice and MS schemes just correspond to a redefinition of µ, but quite large indeed. On the other hand the final result is µ independent. Nevertheless, the way the µ → ∞ limit is taken is fixed  by N A , as defined in Eq. (28), which is dependent on µ. This explains why different results are obtained.
In Figs. 9 and 10, we also compare with results obtained using method 1), more specifically we compare with V PV −V P − 1 r Ω V , as they both have analogous power accuracy (though method 1) is parametrically more precise). For Ω V we take the exact expression but using its approximated expression does not change the discussion, as the difference is very small.
What we see is that the MS scheme yields more precise predictions than the lattice scheme, and that method 1) yields considerable better results than method 2B).
Another issue specific of method 2B) is to determine how large we need to take N (and consequently µ) of the truncated sum such that it approximates well V A . For illustrative purposes we show the convergence in Fig. 11 for n f = 3 in the lattice and MS scheme. We find that we have to go to relatively large values of µ (and N ) to get it precise. This can be a problem if one wants to go beyond the large β 0 . This problem would be less severe if one can use the asymptotic expression for the coefficients beyond certain n. Nicely enough, we find that the use of asymptotic expression for the coefficients for n > N * (∼ 4 in the MS and ∼ 9 in the lattice scheme) is very efficient and basically yields the same results as the exact result. Finally, we also remind that to approximate well V A by the truncated sum is more costly for small values of c .

IV. CONCLUSION
We aim to accurately describe observables characterized by having a large scale Q Λ QCD . For those it is believed that the OPE is a good approximation (we do not enter in this paper on the issue of duality violations). We want to make the most of available perturbative expressions of the observable. Our aim is to organize the computation and its associated accuracy within an hyperasymptotic expansion. For this, we carefully study the connection between truncated sums of the perturbative expansions in powers of α and the associated NP corrections. In practice, we relate those truncated sums with the Borel sum of the perturbative series regulated using the PV prescription. This object has the nice properties of being scale and scheme independent. It may also open the window to connect with studies directly aiming to the NP regime. We then hypothesize that the difference between the Borel sum and the full NP evaluation of the observable complies with the structure of the NP OPE (at least for the first terms of the NP power expansion). Such computational scheme allows us to get an hyperasymptotic expansion of the observable, and, consequently, to unambiguously state the magnitude of the different terms of the hyperasymptotic expansion.
Relating truncated sums of the perturbative expansion with NP definitions of them is not trivial in general. However, this is possible for the case of the PV prescription. We have studied two methods that achieve this goal and explored how reliable they are in practice.
We have given analytic formulas (with exponential accuracy) that relate the truncated sum with the PV-regulated Borel sum. We emphasize that these formulas are valid beyond the large-β 0 approximation.
These methods allow us to efficiently disentangle the pure perturbative term from the first NP corrections of an arbitrary observable that admits an OPE at large energies. General expressions for arbitrary observables are given (for this paper we neglect ultraviolet renormalons). Nevertheless, the accuracy we achieve for each case is different: • The method 2B) (see Eq. (76)) has the handicap that (in principle) needs the perturbative expansion of the observable and the running of α to all orders. On top of that we are only able to obtain the O(e − 2πd β 0 α(Q) α − dβ 1 2β 0 (Q)) term of the Borel sum, which then sets the precision of the analysis. On the other hand, it has the nice feature that the leading NP power correction of the Borel sum has exactly the same scaling as the NP corrections dictated by the OPE, and that the result is explicitly µ independent.
• On the other hand, method 1) (see Eq. (75)) shows to be much more powerful. At low orders it is just standard perturbation theory. At high orders (quantified by N P ) the series is truncated. This corresponds to the superasymptotic approximation. We can quantify the error committed in summations truncated at the minimal term and state the independence of the result on the scale and scheme used for the perturbative expansion to a given accuracy. This allows us to state the parametric accuracy of determinations of genuine NP power corrections obtained by subtracting the perturbative series from the full observable (the latter being obtained either from lattice simulations or directly from experiment).
We then incorporate the NP corrections to the truncated sum associated to the renormalons using the PV regularization prescription. The procedure uses the theory of terminants discussed in [14]. The scale and scheme dependence of this merging is under control in the whole process. This process is, in principle, systematically improvable. Subleading power corrections can be incorporated in the analysis, reaching hyperasymptotic accuracy. This analysis also allows us to visualize that truncating the perturbative sum at the minimal term produce, in general, terms that can not be absorbed in the NP terms of the OPE, because of prefactors proportionals to √ α.
Overall, one obtains an smooth connection between the standard (pure) perturbative computation and the OPE (hyperasymptotic) expansion that includes the NP power corrections.
With these methods it is possible to determine the leading difference between the perturbative series truncated at the minimal term with the Borel integral regulated using the PV prescription in terms of the closest singularity to the origin of the Borel transform. This is very good because it allows us to determine such leading NP correction in terms of the normalization of the leading renormalon, Z X O d , for which approximate determinations can be obtained if the perturbative series is known to high enough orders. It is also worth mentioning that the dependence on Z X O d of the hyperasymptotic approximation to the Borel sum is minimal, since it only appears in Ω. Finally note that there is no need of introducing an infrared cutoff ν f . We plan to apply these methods to general observables, but before we want to study the methods in test-objects for which the approximations are under control. In this paper we take the static potential in the large β 0 approximation, regulating the asymptotic perturbative expansion using the PV prescription, as the observable. It has nice properties: A lot of analytic control is known for it, its Borel transform is known exactly, and it does not have ultraviolet renormalons. In this case we know what the genuine NP corrections are. They are zero by construction.
Whereas the general expressions we give in this paper are valid for any scheme, for the specific analysis worked out in this paper (the static potential in the large β 0 approximation), we use two different schemes: the lattice and the MS schemes. In the large β 0 this is equivalent to a redefinition of the renormalization scale. Nevertheless, let us stress that it corresponds to a rather large change in the scale. Different values of c (see Eq. (27)) can also be understood as a change in the renormalization scale. The result is independent on the scheme and factorization scale used for the α (within the error of the computation). The scheme/scale dependence is a higher order effect. The important thing is that both schemes converge. This does not mean that all schemes converge equally fast. We observe that MS appears to be more convenient for method 2B). It is also interesting to see the dependence of the observables/methods with n f . Indeed we observe that the range of validity of the hyperasymptotic expansion is sensitive to the value of n f . Changing from n f = 0 to n f = 3 significantly enlarge the range of validity of the OPE. This is a relevant discussion when trying to determine up to which scale one can apply perturbation theory and the OPE.
Concerning how well method 1) and 2B) perform in practice for this observable, we find that both methods converge to the expected result. Method 2B) is not particularly precise though. Method 1) appears to converge faster (besides being systematically improvable).
Finally, and specific to method 2B), one issue that we address is how large the renormalization scale µ has to be such that the perturbative expansion simulates well the truncated integral in Eq. (68). For the case of the static potential in the large β 0 approximation, we observe that we have to go to relatively high scales. This makes this method not very useful.
The application of these analyses to QCD observables (beyond the large β 0 approximation) is left to forthcoming papers. We define

Acknowledgments
where x > 0. Note that this integral has a cut in the integration line starting at u = 1.
We have to define how we handle the singularity. We demand D b (−x) to be real for real and positive x. The first expression can be understood as the analytic continuation in b of the second expression (which is first defined for arbitrary positive integer values), and in the second expression we use the PV prescription. Both expressions produce the same asymptotic expansions. Finally, we obtain the following expression where (Γ(b) ≡ Γ(b, 0)) is the incomplete Gamma function. The second term in Eq. (A2) is explicitly real, not so for the first term. Note that the last term in Eq. (A2) is proportional to Λ QCD . From these expressions is difficult to take the b → 0 limit. It is more convenient to set b = 0 before computing.
is defined in [14]. Variants of that formula read (originally generated with a > 0) (A5)