On renormalons of static QCD potential at $u=1/2$ and $3/2$

We investigate the $u=1/2$ [$\mathcal{O}(\Lambda_{\rm QCD})$] and $u=3/2$ [$\mathcal{O}(\Lambda_{\rm QCD}^3)$] renormalons in the static QCD potential in position space and momentum space using the OPE of the potential-NRQCD effective field theory. This is an old problem and we provide a formal formulation to analyze it. In particular we present detailed examinations of the $u=3/2$ renormalons. We clarify how the $u=3/2$ renormalon is suppressed in the momentum-space potential in relation with the Wilson coefficient $V_A(r)$. We also point out that it is not straightforward to subtract the IR renormalon and IR divergences simultaneously in the multipole expansion. Numerical analyses are given, which clarify the current status of our knowledge on the perturbative series. The analysis gives a positive reasoning to the method for subtracting renormalons used in recent $\alpha_s(M_Z)$ determination from the QCD potential.


Introduction
For a long time the static QCD potential V QCD (r) has been studied extensively, in order to understand the nature of the strong interaction between a heavy quark and antiquark pair. In the past decades computation of V QCD (r) in perturbative QCD has been advanced significantly [1,2,3,4,5,6,8,9,10,11,12,13,14]. In association with many theoretical developments, V QCD (r) has become an indispensable theoretical tool to describe not only the properties of the heavy quarkonium states but also for precision determinations of the fundamental parameters of QCD such as the heavy quark masses m c , m b , m t [16,17,18,19,20,21,22,23,24,25] and the strong coupling constant α s [26,27].
Before around 1998, the prediction of V QCD (r) in perturbative QCD was not successful and was plagued by the so-called renormalon problem. As it turned out, convergence of the perturbative series was fairly poor, such that a meaningful prediction could not be obtained in the distance regions relevant to the charmonium and bottomonium states. This is caused by the growth of α s in the infrared (IR) region and is characterized by the singularities in the Borel transform of the perturbative series (the singularities are called renormalons) [28]. Then it was discovered that the leading renormalon of order Λ QCD in V QCD (r) is canceled against that of the quark pole mass M in the combination of the total energy of the static quark pair 2M +V QCD (r), which led to a dramatic improvement of convergence of the perturbative series [29,30,31].
Up to now, although there exists no rigorous proof on existence of renormalons in QCD observables, there exist standard arguments based on the operator product expansion (OPE) and renormalization group (RG) equations which show that their existence is consistent and plausible theoretically [28]. This is reinforced by a number of evidences in actual computation of perturbative series of QCD observables, thanks to recent technological developments in multiloop calculations. There also exist examinations of the nature of renormalons using many approximate estimates of higher-order terms of perturbative series at various levels of rigor. See, for instance, Ref. [32].
In many analyses of renormalons in the static QCD potential, analyses of perturbative computation in momentum space play important roles [28]. For V QCD (r), it is often assumed that there are no renormalons in its Fourier transform V QCD (q) (the potential in momentum space) or that renormalons in V QCD (q) are negligible at the current level of accuracy. In fact, in Refs. [30] and [31], absence of the order Λ QCD /q 3 renormalon in V QCD (q) (corresponding to the u = 1/2 renormalon) is shown at the one-loop and two-loop levels, respectively. 1 Also, the u = 3/2 renormalon cancellation within the multipole expansion was shown [7] based on the assumption of absence of the corresponding renormalon in V QCD (q). Nevertheless, it can be the case that renormalons arise from a deep level of loop integrals in the computation of V QCD (q) and that they are simply not detected in the currently known several terms of the perturbative series.
A direct motivation of our study comes from necessity for a justification for the assumption used in a recent determination of α s from V QCD (r) [27]. There, the first two renormalons of order Λ QCD and Λ 3 QCD r 2 are subtracted from V QCD (r), in order to extend validity range of the OPE of V QCD (r) to larger r, and it is assumed that the corresponding renormalons are negligible in V QCD (q). This problem is also linked with how we renormalize the IR divergences in the potential which arise from three-loops and beyond [33,6,8,13,14,9]. In this paper we analyze the order Λ QCD and Λ 3 QCD r 2 renormalons in V QCD (r), on the basis of the standard argument by the OPE and RG equations. The discussion is to a large extent based on general features of QCD, independent of ad hoc approximations such as the large-β 0 approximation. We refine our understanding by looking into the detailed structure of the OPE within the potential-NRQCD (pNRQCD) effective field theory (EFT) [34]. In particular we elucidate the accurate structure of the u = 3/2 renormalon. Subsequently we discuss the size of the renormalon uncertainties for u = 1/2 and 3/2 in V QCD (q) with a method which does not rely on diagrammatic analysis, providing a different perspective from, e.g., Refs. [31,7]. We also believe that an argument such as the one we provide in Sec. 3 is necessary to clarify treatment of the IR part of V QCD (q). In the latter part of this paper, we test our understanding by performing numerical analyses of the normalization constants of the renormalons in the perturbative series of V QCD (r) and V QCD (q). We treat two kinds of perturbative series: one is the fixed-order perturbative series currently known and the other includes higher-order terms estimated by RG. We estimate the normalization constants of the renormalons from these two perturbative series by using Lee's method [35] and also by an analytic formula which we derive (for the latter series). The former includes updates of the analyses by Pineda [36].
The paper is organized as follows. In Sec. 2 we briefly review the standard argument on renormalons for a general QCD observable. In Sec. 3 we scrutinize the structure of renormalons in the static QCD potential. Secs. 4-8 present numerical analyses on the normalization constants of the renormalons. In Sec. 4 we study the O(Λ QCD ) renormalon of V QCD (r) from fixed-order perturbative series, followed by a study of its cancellation with that of the pole mass in Sec. 5. We compare these results with that obtained by an integral formula in Sec. 6. We study the O(Λ 3 QCD r 2 ) renormalon of V QCD (r) in Sec. 7. Finally we test the corresponding renormalons in V QCD (q) in Sec. 8. Conclusions are given in Sec. 9. In App. A we explain theoretical aspects of IR cancellation at O(r 0 ) of the multipole expansion. In App. B we present details of the derivation of a formula for the normalization of renormalons in V QCD (r).

Structure of renormalons
Let us first review briefly the structure of renormalons in QCD observables [28].
Consider a general RG-invariant dimensionless observable X(Q) with a typical energy scale Q. Its perturbative expansion is given by µ denotes the renormalization scale in the MS scheme. It satisfies the RG equation with the beta function given by The first two coefficients of the beta function are given explicitly by It is conjectured that for many observables the coefficients of the perturbative series grow factorially, d n ∼ n!, for large n. To quantify uncertainties induced by this property, the Borel transform of X PT , defined by is studied. Renormalons of X PT refer to the singularities of B X (t) located on the real axis in the complex t-plane. We assume the form of the Borel transform in the vicinity of each renormalon singularity at t = b 0 /u as 2 with parameters N u , ν u and c k 's. Formally we can reconstruct X PT from its Borel transform B X (t) by the inverse Borel transform given by the integral However, if there are singularities (renormalons) on the positive real axis, the integral is ill defined. We can regularize the integral by deforming the integral contour to the upper or lower half plane: We can define the ambiguity induced by the renormalon from the discontinuity of the corresponding singularity, and the singularity at t = u/b 0 with u > 0 of eq. (6) gives where we have used The parameters u, ν u , c k andc k in eq. (6) or eq. (10) can usually be determined from the OPE. In the context of the OPE in 1/Q, X PT is identified with the Wilson coefficient of the leading identity operator. Let us denote by O u the lowest dimension (dimension 2u) renormalized operator responsible for cancellation of the renormalon in X PT . 3 The OPE reads We assume that the leading ambiguity induced by the renormalon of C X 1 as given in eq. (10) is canceled by the second term of the OPE. Then, the Q-dependence of the renormalon uncertainty of C X 1 should coincide with that of the second term in the OPE, which can be detected as follows. Suppose that the Wilson coefficient satisfies the RG equation, This RG equation specifies the Q-dependence of C X Ou (Q/µ, α s (µ)) as where const. denotes a Q-independent (but µ-dependent) constant. Now the Q-dependence of the second term of the OPE (12) is made explicit. Requiring the same Q-dependence for the renormalon uncertainty of X PT , using eq. (12) with eq. (15), we obtain for eq. (6) or (10) The factor ∞ k=0c k α s (Q) k in eq. (10) should be proportional to [1+O(α s (Q))]C X Ou (1, α s (Q)) in eq. (15). Therefore, c k 's andc k 's can be determined one by one from smaller k in terms of b n 's, γ n 's and f n 's from smaller n. The overall normalization N u cannot be determined from this argument. We note that, in the case that C X Ou O u is independent of Q, γ 0 = 0 andc k = 0 for k ≥ 1.

Renormalons in the static QCD potential
In this section we investigate theoretical aspects of renormalons of the static QCD potential, focusing on the u = 1/2 and 3/2 renormalons, on the basis of the above general understanding. Part of the argument given in this section has already been discussed in [6]. (See also [37].) We refine the discussion and present new observations. In particular, main part of the discussion on the u = 3/2 renormalon is new.

Basics of static QCD potential
The static QCD potential is defined from an expectation value of the Wilson loop as where C is a rectangular loop of spatial extent r and time extent T . P stands for the path-ordered product along the contour C. For convenience we define the dimensionless potential as v(r) = rV QCD (r). It is conjectured that renormalon singularities are located in the Borel transform of v(r) at t = 1 2b 0 , 3 2b 0 , etc. (i.e., u = 1 2 , 3 2 , etc.). In calculation of the static QCD potential, we have two different scales. One is the soft scale 1/r, which is the inverse of the distance between a static QQ pair. The other is the ultrasoft (US) scale, which is set by the energy difference of the color singlet and octet states of the static QQ pair, where The pNRQCD EFT describes dynamics in which the QQ system emits or absorbs US gluons whose energies are comparable to or smaller than the energy differences of different QQ states [7]. Accordingly, the factorization scale µ f (= cut off scale of pNRQCD) is chosen to satisfy ∆V ≪ µ f ≪ 1/r. The OPE of the static QCD potential V QCD (r) in r can be performed within the pNRQCD EFT in the static limit based on the scale hierarchy 1/r ≫ ∆V : 4 V QCD (r) = V S (r) + δE US (r) + . . . , The leading term V S (r) denotes the singlet potential, which is the Wilson coefficient of the bilinear singlet field operator S † S in the context of pNRQCD. In eq. (22), V S (r) is 4 This is equivalent to 1 2 C A α s (1/r) ≪ 1, which holds at sufficiently small r. multiplied by S | S † S | S = 0 | 1 | 0 = 1. The second term is the O(r 2 ) term in the multipole expansion, given by where represents the color electric field; the color string for the adjoint representation is given by ϕ adj (t, t ′ ) = T exp ig t t ′ dτ A c 0 (τ, X) T c adj . δE US is generated by insertions of the operators g O a † r · E a S and g S † r · E a O a , and V A (r) denotes the Wilson coefficient of these operators. Note that eqs. (22) and (23) are exact to all orders in α s . V S (r) coincides with the naive expansion of V QCD (r) in α s : To see this, we adopt the energy integral representation of δE US (r), which can be obtained with Fourier transform of the correlation function, Then, if we naively expand δE US in α s before loop integrations, the US scale ∆V = O(α s ) disappears from the propagator denominator in eq. (25), and the integrals become scaleless and vanish. (The same applies to beyond O(r 2 ) terms.) We can rephrase this in the computation of V QCD (r) in expansion in α s , by applying expansion-by-regions technique to loop integrals [38]. We can separate contributions from the UV scale 1/r and the US scale (≪ 1/r), where the latter contributions vanish to all orders in α s since they are given by scaleless integrals. We investigate theoretical aspects of renormalons at u = 1/2 (Sec. 3.2) and u = 3/2 (Sec. 3.3) based on the above general understanding, and in particular determine some of the parameters in eq. (6) assuming this expansion form around the singularities. Before this, let us comment on the IR divergences present in the perturbative result of V S (r). The naive perturbative expansion of V QCD (r) includes IR divergences at and beyond order α 4 s [33,6,8,13,14,9], hence so does V S (r). The IR divergences of V S (r) have their counterparts in the OPE at order r 2 or beyond in the multipole expansion. Indeed, δE US contains UV divergences if we compute it in double expansion in α s and log α s consistently with the philosophy of pNRQCD, that is, keeping ∆V ( > ∼ k g ) in the propagator denominator. At O(r 2 ), the UV divergences of δE US and IR divergences of V S (r) cancel in V QCD (r), reflecting the µ f -independence of V QCD (r). In the subsequent argument, we implicitly assume a certain regularization prescription for making these divergences finite to discuss renormalons in the perturbative series whose each expansion coefficient is finite. We will propose explicit regularization (renormalization) schemes and also discuss their relevance to the renormalon structure (Sec. 3.5) after the renormalon structure is clarified.

u = 1/2 renormalon
Let us clarify the current understanding on the u = 1/2 renormalon. The leading IR renormalon of V S (r) is located at u = 1/2, and the induced ambiguity is known to be independent of r and proportional to Λ MS . In fact, the r-independent constant part of V QCD (r) in eq. (17) is not well defined. This is inherent in the self-energy type contributions Σ to each static color charge. These contributions vanish in perturbative computation in dimensional regularization, since they are given by scaleless integrals. Hence, they are not included in the computation of V S (r), which consists of the potentialenergy type contributions (represented by diagrams with cross talks between the two static charges). See Fig. 1. The IR contributions to the self-energies 2Σ cancel against the IR contributions to V S (r). This is represented in pNRQCD by the absence of O(r 0 ) interactions of the singlet field and US gluon field and is a consequence of the fact that in the IR limit gauge field couples to the total charge (= 0 for | S ); further explanation is given in App. A. On the other hand, Σ is UV divergent, and in dimensional regularization simply Σ is set equal to zero. Thus, more precisely, V S (r) should be written as 2Σ+V S (r) in eq. (22), but Σ is omitted in accordance with the usual convention.
A standard way to confirm cancellation of the r-independent IR contributions to V S (r) with the self-energy type contributions is to show cancellation of the u = 1/2 renormalons in the combination 2m pole + V S (r) [29,30,31]. By construction of pNRQCD as a low energy EFT, IR contributions to Σ and m pole are common. Both Σ and m pole are RG invariant, hence ambiguities induced by the leading renormalons both correspond to u = 1/2 and proportional to Λ MS . This reasoning determines the parameters in eq. (10) and consequently in eq. (6) as

u = 3/2 renormalon
To clarify the structure of the u = 3/2 renormalon, the r-dependence of the u = 3/2 renormalon uncertainty should be revealed (as done in Sec. 2). To this end, we focus on δE US (r) [for instance, the expression of eq. (25)], which cancels the corresponding renormalon uncertainty of V S (r). At this stage, we note that computation of V S (r) does not include the US scale, and thus the u = 3/2 renormalon uncertainty is independent of ∆V . This reasoning and the expression, for instance, of eq. (25) tell us that the rdependence of the u = 3/2 renormalon uncertainty is given solely by ∼ From the explicit NLO result V A (r) = 1 + O(α 2 s ) [10], we see that γ 0 = 0. Thus, we determine the parameter ν 3/2 in the Borel transform in eq. (6) as [cf. eq. (16)] We also clarify that the u = 3/2 renormalon uncertainty is given by eq. (10) with The parametersc k can be parametrized by e i 's, γ i 's, and b i 's. With the NLO result of V A ,c 1 is explicitly obtained as Thus, a correction factor to the r-dependence of the renormalon uncertainty r 2 is given by 1 + O(α s (1/r)) with the abovec 1 term.

Renormalon in momentum-space potential
We now discuss the renormalon uncertainty in the momentum-space potential. Suppose that we have the ambiguity in the position-space potential due to the renormalon at with For u = 1/2, eq. (32) is exact, whereas the correction factor 1 + O(α s (1/r)) arises for u = 3/2. The momentum-space potential is obtained by the Fourier transform of V QCD (r): From the Fourier transform of v ± , we can obtain the renormalon uncertainty in the V -scheme coupling constant in momentum space: 5 Here we are concerned with the logarithms associated with the UV divergences of V A (r) in the full theory (or with respect to the soft scale). This RG equation with respect to µ is different from the RG equations with respect to µ f considered in Refs. [10,11,41], which are associated with the IR divergences with respect to the soft scale. See discussion in Sec. 3.5. 6 We implicitly assume that the u = 3/2 renormalon uncertainty is RG invariant as we assume eq. (10), which is indeed RG invariant, although RG invariance of V S (r) may be violated by the IR divergences (or IR logarithms). This assumption is justified when we adopt explicit schemes to remove the IR divergences from V S (r) such that the redefined V S is RG invariant; see Sec. 3.5.

Renormalization scheme
So far, we did not specify how to renormalize V S (r) and δE US (r), which contain the IR divergences and UV divergences, respectively, Here, we define two schemes to remove the divergences.

Scheme (A)
At each order of the perturbative expansion of V S (r) in α s , we first set µ = 1/r and then drop all the poles in ǫ originating from the IR divergences. (µ denotes the renormalization scale in full QCD.) We also redefine δE US such that the sum V S (r)+δE US is unchanged, which is evaluated in double expansion in α s and log α s . The renormalized V S and δE US are both independent of µ by definition. 9 In fact, this regularization is compatible with the property used in Sec. 3.3 that V S (r) is RG invariant (see footnote 6), but it may be incompatible with the one that the r-dependence of the u = 3/2 renormalon uncertainty is given by ∼ r 2 V 2 A (r). This is because the latter reasoning [and thus the results such as eq. (30)] relies on eq. (25) and additional contribution was not considered. However, we assume that the structure of IR renormalons in V S (r) at O(r 2 ) is unchanged by this prescription to remove IR divergences. This is indeed the case in the large-β 0 approximation of δE US , in which the IR divergences and IR renormalons are clearly separated; the former is given as 7 An alternative method is to use dimensional regularization and to utilize analytical continuation with respect to the dimension D = 4 − 2ǫ. 8 The uncertainty (37) is obtained in a parallel form to eq. (10) in the sense that the part given by the series expansion in α s is specified with b i 's and γ i 's. Thus, the result sounds plausible. In particular, c 1 is given byc 1 = −γ 1 /b 0 . 9 V S (r) and δE US (r) at different µ are obtained by rewriting α s (1/r) by α s (µ).
a convergent series in α s expansion, while the latter is given as a factorially diverging series. This is shown by computing δE US in the large-β 0 approximation [40]: where .
The UV divergences and UV renormalons of δE US (r) are canceled, respectively, by the IR divergences and IR renormalons of V S (r). In addition, the US logarithms at LL [10] and NLL [11,41], associated with the IR divergences in V S (r), are known to be given by convergent series in α s expansion, computed explicitly using the RG equation of pNRQCD. Thus, to the best of our knowledge, the above assumption seems to be valid. As a result, we consider that the scheme (A) is suitable for studying the renormalon of V S (r), where the renormalon structure revealed in Sec. 3.3-3.4 based on the OPE is not expected to be modified. We note the existence of the UV renormalons at u = 3/2, 1/2, −1/2, −3/2, · · · of δE US in eq. (39). It is confirmed that the leading UV renormalon at u = 3/2 cancels the known O(r 2 ) IR renormalon of V S (r) in the large-β 0 approximation [40,42]. The subleading renormalons are also expected to be canceled against V S (r) (although it cannot be confirmed within the large-β 0 approximation) since the IR structure of V S (r) should match the UV structure of δE US . The residues of the subleading renormalons at smaller u are proportional to higher powers of ∆V (r)/µ. This leads to less powers 10 of r, which contradicts to the naive expectation that the renormalons of V S (r) beyond u = 3/2 are suppressed by higher powers of r in accordance with the multipole expansion. This feature originates from the fact that if we expand 11 the integrand of eq. (25) in ∆V , higher power singular IR behaviors ∼ (∆V (r)/k) n appear. The IR structure of V S (r) includes the same power behaviors, since the IR structure of δE US is common to that of V S (r) once the integrand is expanded in α s . The higher power singular IR behaviors generate the above more singular IR renormalons as well as higher power IR divergences. 12 13 10 Since ∆V (r) ∼ (r| log r|) −1 , the form of the renormalons are not integer powers of r. 11 Note that E a i ϕ ab adj E b j (k) is independent of ∆V . 12 The poles on the negative axis will not generate an essential problem since they are Borel summable. The pole at u = 1/2 induces a serious problem. 13 Up to date, these more singular IR renormalons have not been investigated seriously. One reason would be that they are generated only at higher loops, since they arise with higher powers of ∆V . In this connection, we note that the UV renormalons at u < 3/2 in δE US | large-β0 do not have their IR renormalon counterparts in V S | large-β0 because the order counting is different between these quantities. The former are suppressed by higher powers of ∆V ∼ O(α s ) compared to the latter. We need to go beyond the large-β 0 approximation to detect these IR renormalons in V S (r).
The above observation in particular means that V S (r) has a renormalon at u = 1/2 corresponding to the above UV renormalon of δE US . The u = 1/2 renormalon uncertainty is given by O(Λ MS r 2 ∆V (r) 2 ), as seen from eqs. (38) and (39). Thus, this renormalon is different from the familiar renormalon at u = 1/2, which induces an rindependent uncertainty. (See footnote 10.) We note that this unfamiliar renormalon at u = 1/2 can be an obstacle in estimating the normalization constants of the familiar renormalons at u = 1/2 and 3/2. This possibility is taken into account later in numerical analyses, while we also present naive analysis by simply neglecting this peculiar renormalon.

Scheme (B)
We subtract IR divergences from V S (r) by adding δE US evaluated in double expansion in α s and log α s (Scheme B1). In this scheme, we do not distinguish V S (r) and δE US , and we exclusively treat the sum of them, which is regarded as a redefined V S (r). In this way we can subtract the IR divergences. Furthermore, after canceling the IR divergences, we can replace the argument of US logarithms as log(r∆V ) → log(rµ f ) (Scheme B2). Since both V QCD (r) and ∆V (r) are µ independent, the renormalized V S (r)'s are also µ independent up to O(r 2 ) (although V Finally we point out that it is not straightforward to cancel simultaneously both the IR divergences and IR renormalon at u = 3/2 of V S in the sum V S + δE US . We can observe this in the large-β 0 approximation. The renormalon uncertainty of δE US coincides with minus that of V S (r) when ∆V in δE US is not expanded in α s . If ∆V is perturbatively expanded instead, the power of α s shifts by three in the perturbative series due to ∆V = O(α s ), as seen from eq. (38), 14 and the renormalon cancellation breaks down. Hence, it would be optimal not to expand ∆V in α s for the renormalon cancellation. On the other hand, this prescription is not preferable to cancel the IR divergences. The IR divergences are cancelled when ∆V is expanded in α s as the IR and UV divergences in V S and δE US , respectively, appear at O(α 4 s ). The proposed two schemes above can remove the IR divergences from V S (r), but cannot remove the IR renormalons of V S (r). It remains a future task to develop a method for subtracting the IR renormalons completely. 15 14 The normalization of the renormalon is changed by the expansion, which ruins the renormalon cancellation. 15 Suppose that we can remove completely the u = 1/2 and 3/2 IR renormalons from V S defined in the scheme B in some way. The remaining renormalons are proportional to (r 2 Λ 4 QCD /∆V ) × (Λ QCD /∆V ) n (n ≥ 0) or that with ∆V replaced by µ f . They are obtained by expanding the correlator of eq.(23) in t. In particular the leading IR renormalon (n = 0) is given in terms of the local gluon condensate [43,44,45,7], or ∆V replaced by µ f . Thus, the leading renormalon in V (B1,B2) S (r) is located at u = 2 and suppressed by Λ QCD /∆V or Λ QCD /µ f compared to the original u = 3/2 renormalon in V S (r).

Renormalon subtraction by contour deformation
One motivation of the above investigation is to give a justification to the method used to subtract the u = 1/2 and 3/2 renormalons from V QCD (r) in a recent determination of α s (M Z ) [27]. There, it is assumed that the corresponding renormalons contained in V QCD (q) can be neglected. (The IR divergences are canceled in momentum space.) Then, using the one parameter integral form with respect to the momentum q and deforming the integral contour in the complex q-plane, the renormalons at u = 1/2 and 3/2 which stem from the original q integral are subtracted. See [47,27] for the details. As we have seen above, the normalization of the u = 1/2 renormalon in the momentum-space potential is exactly zero, while the u = 3/2 renormalon is suppressed by α s (q) 2 . While the r 2 Λ 3 MS renormalon that is generated purely by the q integral is subtracted, the suppressed renormalon inα V (q) can still contribute to the position-space potential. That is, ifα V (q) does not converge, its uncertainty will give an uncertainty to the renormalon-subtracted V S (r) constructed by the contour deformation method. It is expected to generate a renormalon of order r 2 Λ 3 MS α s (1/r) in the renormalon-subtracted V S (r), corresponding to the correction proportional toc 1 of Sec. 3.3.

Numerical study of u = 1/2 renormalon
In the rest of this paper we perform numerical analyses on the normalizations N u of renormalons to check the above observations and to see the current status of our knowledge on the perturbative series for V QCD (r) and V QCD (q). We treat two perturbative series: one is the fixed order perturbative series, and the second one includes higher-order terms estimated by RG, which is used extensively in Refs. [46,47].
In the case that the renormalon singularity of the Borel transform is given by eq. (6), we can estimate the normalization constant N u from the fixed-order result of the perturbative series. This was first proposed in Ref. [35], whose method is as follows. We consider the function (In the following analyses, ν u of eqs. (27) and (30) are used.) We can obtain the normalization constant by expanding this function in t and then substituting t = u/b 0 : as long as the corresponding renormalon is the one closest to the origin. 16 This method is fairly general and can be used with only known terms of the perturbative series. Using this method, Ref. [36] studied N 1/2 from the fixed-order result. We present an The IR divergence at NNNLO is subtracted in the scheme (A). 17 We also examine the scale dependence of the estimated normalization constant. We use the perturbative coefficients with the renormalization scale µ = s/r to estimate the normalization constant. 18 The results are shown in Fig. 2. The scale dependence decreases as we include higher-order terms. These results indicate that the series (43) shows convergence for the u = 1/2 renormalon and N 1/2 (s = 1) ≃ −1.1.
In the second method we estimate the higher-order terms by RG. Explicitly we can use [46,47] where the N k LL terms of the perturbative series [coefficients of α s (µ) n+k+1 log n (µr) for arbitrary n] can be determined using the (k + 1)-loop beta function and the fixed-order result up to k-loops. We now estimate N 1/2 from the RG improved series obtained from eq. (45) using eq. (43). Since we have an all-order perturbative series (at each order of improvement), we can obtain N LL , · · · , N N 3 LL with arbitrary precision in principle. The results from 17 If we adopt the scheme (B2), the NNNLO result is given by N 1/2 = −1.07764−0.0541542 log(2rµ f ). 18 When the scale µ = s/r is used in constructing the Borel transform, the normalization constant N (s) of the renormalon at u behaves as N (s) = N (s = 1)s u as seen from eq. (6). The s-dependence of the estimated result for s −u N (s) is expected to reduce as we include higher-order terms. We always consider N (s = 1) unless stated otherwise explicitly.
We use the scheme (A) to subtract the IR divergence in the NNNLL analysis. 19 From Table 1, which shows the convergence speed, we infer that 20-50 perturbative coefficients are needed in order to obtain the normalization constants with one-percent accuracy.

Renormalon cancellation in total energy
It is interesting to examine renormalon cancellation in the total energy (namely, V (r) + 2m pole ) from the estimated N 1/2 [36]. The leading renormalon in the Borel transform of m pole /m MS is given by It is possible that treatment of IR divergences affects the cancellation. Let us examine this. So far, we subtracted the IR divergence in the scheme (A), but now we make the three-loop coefficient finite in the scheme (B2). Then, the IR divergence is replaced by the logarithmic term like log(µ f r). In Fig. 3, we investigate renormalon cancellation while varying µ f in this logarithm in a reasonable range. This figure shows that the treatment of IR divergences can be non-negligible to the precise cancellation. We note that the numerical error on the four-loop result of the mass relation hardly affects this result.

Normalization by analytic formula in RG-improved method
In the RG-improved method, we derive a formula for the normalization constants of renormalons given as a one-dimensional integral. The Borel integral of the QCD potential in the RG-improved method can be written as where the contours C ∓ are displayed in Fig. 4. The details of the derivation are given in App. B. Since the integrand satisfies {f (x)} * = f (x * ), the imaginary part can be 20 Note that the perturbative series needs to be expressed in terms of the coupling of the theory with n l light quarks only, while originally the pole mass is expressed by the coupling in the theory with n l light quarks plus one heavy quark. This is needed to ensure the renormalon cancellation, since N M and N 1/2 are proportional to Λ MS and the same Λ MS should be used for both quantities. (In principle one can pursue the calculation in the different couplings if the difference in the definitions of Λ MS is properly taken into account.)

calculated by a contour integral
where the contour C is displayed in Fig. 4. By expanding sin(qr) in qr, the normalization of the renormalon at u = 1/2 is found as In the last line, we changed the integration variable to x = q/Λ MS . (Note that α V (q) is a function of q/Λ MS .) Then, from eqs. (10) and (52), we obtain with We present numerical values of N 1/2 via numerical evaluation of I 1/2 : They agree well with the estimates from the finite number of terms (46). The scheme (A) is adopted at NNNLL in accordance with Sec. 4. It is possible to calculate the normalization constants of other renormalons in a parallel manner. The normalization constant of a general renormalon at u is expressed as This expression stems from the Taylor expansion of sin(qr). In this method, the r-dependence of a renormalon uncertainty due to any half-integer renormalon at u is given exactly by r 2u+1 . In particular for u = 3/2, the correction factor of [1 + O(α s (1/r))] is not detected (as long as we work at N k LL with finite k). It is because this method relies on the assumption that α V (q) does not possess renormalon uncertainties.
For the u = 3/2 renormalon, which is the second IR renormalon, the numerical values of N 3/2 are given by

Numerical analysis of u = 3/2 renormalon
We now estimate N 3/2 from the fixed-order result. We annihilate the leading renormalon at u = 1/2, whose uncertainty is an r-independent constant, by considering the QCD force. Then we use the same method as in the u = 1/2 renormalon. We first examine the relation between the normalization constants of the potential and force. The potential v = rV has the u = 3/2 renormalon uncertainty as eq. (10) with eqs. (30) and (31), which gives the uncertainty to the dimensionless force f = r 2 dV /dr as Thus, the normalization constant of the dimensionless force N F 3/2 is related as To obtain N 3/2 , we first consider the fixed-order perturbative series of the potential without setting µ = 1/r. The derivative with respect to r gives the fixed-order result of the force. Finally we set µ = 1/r to estimate N F 3/2 and translate it to N 3/2 . We present the results: We adopt the scheme (A) to obtain the NNNLO result. 21 Although the estimate of the normalization constant may look already convergent, this seems to be a numerical accident. We examine the scale dependence of the estimated normalization constant in a parallel manner to the u = 1/2 renormalon. The result is shown in Fig. 5. We find that large dependence on the renormalization scale remains in this estimate. Let us perform a parallel estimate from the finite order result in the RG improved method. In this case, since we know N 3/2 as given by eq. (58), it would be useful to grasp how many terms are needed for a reasonable estimate. Table 2 shows the result. At NNNLL, we adopt the scheme (A). One sees that typically 20 terms are needed for a good estimate.
We examine scale dependence of the estimate of N 3/2 using finite number of terms in the RG improved scheme. Since we know the exact answer in this case, we can directly check whether mild scale dependence indicates reliability of the estimate. In Fig. 6, we examine this at NNLL. We see that at higher order the scale dependence of N 3/2 decreases and it approaches the correct value. (At further higher order, for instance from 30 terms, we obtain N 3/2 = 0.0079136 and 0.079172 for s = 1/2 and 2, respectively.)  Table 2: Estimates of the normalization constant N 3/2 from truncated perturbative series using the RG-estimates for the higher-order terms. In Sec. 3.5, we pointed out that V S (r) can have the unfamiliar renormalon at u = 1/2 associated with the IR divergences. Since the corresponding renormalon uncertainty for the static QCD potential is not an r-independent constant in contrast to the familiar renormalon at u = 1/2, this renormalon cannot be eliminated in the QCD force. Taking into account this possibility, we present another estimate for the u = 3/2 renormalon, whose method is not plagued by the two renormalons at u = 1/2. We carry out this by using a mapping from the t-plane to a new z-plane, where the u = 3/2 renormalon becomes closer to the origin than the u = 1/2 renormalon. Namely we have to change the relative distances of the two IR renormalons from the origin. 22 A possible mapping is given by b 0 t(z) = 1 2 (z + e iπ/6 ) 6 + 1 2 .
(61) 22 One may compare with Ref. [35], in which for the Adler function the closer UV renormalon at u = −1 is made farther than the IR renormalon at u = 2 by a mapping.
However, it turned out that with the above mapping convergence is too slow for practical analysis. 23 Instead of Eq. (61), we use This mapping is obtained with a similar idea to the above, but the main difference is that we first consider square of the difference from b 0 t = 2, i.e. (2 − b 0 t) 2 . The mapping (62) consists of the following steps , w(y) = y 6 , y(z) = −z + e iπ/6 . We note that the singularities of 1/(1 − 2b 0 t(z)/3) and 1/(1 − 2b 0 t(z)/5) with respect to z are not common, and the u = 5/2 renormalon does not affect the estimate of the normalization constant at u = 3/2. With this mapping t(z), we consider a function By expanding this function in z and then substituting z = −(8/7) 1/6 + e iπ/6 , we can obtain the normalization constant of the u = 3/2 renormalon. Using this mapping, we estimate the normalization constant of the u = 3/2 renormalon from the fixed-order perturbative series as We adopt the scheme (A) at NNNLO. 24 In this method, imaginary parts appear in fixed-order results, but we omit them in the above estimate since we know that the true normalization is real. The size of the imaginary parts can be used for an error estimate of the results. We also estimate N 3/2 of the RG improved series using this mapping. Table 3 shows the result. One can confirm that the estimated values converge to the results in Eq. (58). We start to obtain reasonable results with about 60 terms.  8 u = 1/2 and 3/2 renormalons in V QCD (q) Let us investigate the renormalon uncertainty of α V (q). We estimate the normalization constants of the renormalons of α V (q) at u = 1/2 and 3/2 assuming that they are the leading renormalon individually. More explicitly we assume The theoretical discussion in Sec. 3.4 shows that the normalization constants N α V ,u defined in this way should be zero both for u = 1/2 and u = 3/2, because the u = 1/2 renormalon is completely absent inα V (q), and for u = 3/2 the expansion of the Borel transform around the singularity takes a form corresponding to the α s (q) 2 suppression. 25 The estimates from the fixed-order results read We subtract the IR divergence at NNNLO in the scheme (A). If we instead use the scheme (B2), the NNNLO results are modified as and N α V ,u=3/2 = −0.535547 + 1.09662 log(2rµ f ) at NNNLO .
By taking rµ f = 0.2 as an example, we obtain N α V ,u=1/2 = −0.1723293 and N α V ,u=3/2 = −1.54037 at NNNLO. In Fig. 7 we show the estimates of N α V ,u=1/2 and N α V ,u=3/2 from the RG-improved perturbative series, in addition to the ones from the fixed order results. In these figures, we subtract the IR divergence in the scheme (A) at NNNLO. (Note that the RG-estimates of the terms beyond N 3 LO are zero for V QCD (q) since we set µ = q.) We also plot their estimates using the large-β 0 approximation for the higher-order terms (they are nonzero even beyond N 3 LO). In both cases we know that the normalization constants are zero. We see in the figures that the estimates approach zero as we include more terms. Since the normalization constants are expected to be zero (even if we do not use any approximation), this figure shows overall consistency.
Thus, in both cases the observed results are consistent with the expectation that the renormalon at u = 1/2 is absent and the u = 3/2 renormalon is suppressed. For u = 1/2, we may already observe smallness of the renormalon contribution from the known perturbative series. For u = 3/2, however, the number of terms are much too few to make any statement on the size of the renormalon. By using the formula (35) and the fact that r-dependence of V A = 1 + O(α 2 s ) is suppressed, we can make a stronger prediction on the smallness of the renormalon. We confirm validity of this formula by using the higher-order estimates by RG-improvement (trivial) or by the large-β 0 approximation.

Conclusions
We have investigated the u = 1/2 and u = 3/2 renormalons in the static QCD potential in position space and momentum space. In particular we have presented detailed examinations of the u = 3/2 renormalon for the first time. In terms of pNRQCD EFT, we have studied the renormalon of the Wilson coefficient V S (r) (and in connection with this the second term of the multipole expansion, δE US , as well).
We have determined the structure of the u = 3/2 renormalon based on the OPE (or multipole expansion) and the RG equations. Although there are non-trivial features specific to the QCD potential (originating from the fact that the multi-scales are involved), we find that the renormalon uncertainty can be parameterized (besides the overall normalization) similarly to the general case as reviewed in Sec. 2. The relevant parameters are the Wilson coefficient of the O( r) interaction V A , in particular its anomalous dimension (associated with the logs from the soft scale), and also the coefficients of the beta function. We have also clarified how the renormalon uncertainties of the positionspace potential propagate to the momentum-space potential. The u = 1/2 renormalon is completely absent in the momentum-space potential, and the u = 3/2 renormalon uncertainty is suppressed by α s (q) 2 in momentum space compared to that in position space. While the renormalon uncertainty of the momentum-space potential has been believed to be small, our result provides a quantitative insight on this issue. We have given a systematic and precise analysis of the old problem, including renormalization prescription and treatment of the IR divergences (US logarithms) based on the multipole expansion in the pNRQCD EFT.
There are some difficulties caused by the IR divergences, however. First, it is not obvious whether the renormalization of V S (r) to remove the IR divergences affects the renormalon structure detected from the OPE argument. We have proposed a way to remove the IR divergences which is likely to keep the renormalon structure unchanged based on our current knowledge. Secondly, we have pointed out that it is difficult to eliminate the IR divergences and IR renormalon at u = 3/2 of V S (r) simultaneously in the multipole expansion, i.e., V S (r)+δE US . In particular, the perturbative result for the sum given by the double expansion in α s and log(α s ) is free from the IR divergences but not from the IR renormalon. A systematic method which can subtract the IR renormalon as well needs to be developed for obtaining an accurate prediction. The contour deformation method used in Ref. [27] has an advantage in this respect (see below).
We performed numerical analyses and checked our understanding as well as the current status of our knowledge on the perturbative series of V QCD (r) and V QCD (q). With the available first four terms of the perturbative series, we find that already the normalization constants of the u = 1/2 renormalons can be estimated with moderate accuracies (consistent with the analyses [36]). On the other hand, the normalization constants of the u = 3/2 renormalons are still not reachable. According to the RG estimates of the higher-order terms (neglecting beyond NNNLL terms), it is suggested that we need 15-20 terms of the series expansion to obtain reliable estimates of the normalization constant. In the same RG method, we obtained an analytic formula for the normalization constants for half-integer renormalons [eq. (57)], which is confirmed to be valid by comparison with the estimate using Lee's method, which utilizes finite number of terms of perturbative series.
We noted the existence of a peculiar renormalon at u = 1/2, which is related to IR divergences of the static QCD potential in naive perturbation theory and induces an uncertainty of O(Λ MS r 2 ∆V 2 (r)). This can be an obstacle in estimating the normalization constants of the familiar renormalon at u = 1/2 and u = 3/2. To investigate the familiar u = 1/2 renormalon (which induces an r-independent uncertainty), it is better to study perturbative expansion of the pole mass in terms of the MS mass, which is free from IR divergence. To study the normalization of the u = 3/2 renormalon, we proposed a method using a non-trivial mapping, which is not disturbed by the renormalons at u = 1/2.
As an application, the present work clarifies the status of the method (contour deformation method) used in a recent determination of α s (M Z ) from V QCD (r) after subtracting the u = 1/2 and 3/2 renormalons [27]. (The IR divergence is canceled as well.) There, it is assumed that the corresponding renormalons contained in V QCD (q) can be neglected. As we have seen in Sec. 3.4, the normalization of the u = 1/2 renormalon in the momentum space potential is exactly zero. For the u = 3/2 renormalon, it turned out that the dominant (or leading) uncertainty ∼ r 2 Λ 3 MS , which comes from the IR region of the Fourier transform of the momentum space potential, is subtracted since this method removes the IR region. On the other hand, the subleading part ∼ r 2 Λ 3 MS α s (1/r), which comes from the uncertainty of the momemtum-space potential, is generally expected to remain. This shows how the u = 3/2 renormalon is suppressed theoretically, and the current status is that with the first four terms of the perturbative series the u = 3/2 renormalon in V QCD (q) is far from detectable, based on the detailed numerical analysis. Therefore, the method is reasonable for subtracting the renormalons from V S (r) using the currently known terms of the perturbative series.
To study renormalons beyond u = 3/2, there still remain works to be done. In particular, it has not been clarified yet which renormalons are specified from the OPE of pNRQCD EFT beyond u = 3/2. Figure 8: To take into account the couplings of a gluon to the total static currents, both self-energy and potential-energy diagram contributions need to be included, and a cancellation takes place between them in the IR limit of the gluon momentum | q| → 0.
Appendix A: IR cancellation at O(r 0 ) The cancellation of IR contributions between the self-energy 2Σ and the potential energy V S (r) is a general property of gauge theory, which can be seen as follows. A static current has only the time component, j µ a,i (x) = ±T a δ µ0 δ 3 ( x ∓ r/2), (i = Q,Q) since a static color charge has no spatial motion. Here, ± r/2 denote the positions of the static charges Q andQ. (We fix the c.m. coordinate to the origin 0.) Hence, an IR gluon, which couples to the static currents via minimal coupling A a µ ( q, t) j µ a,i (− q, t) = A a 0 ( q, t) j 0 a,i (− q), couples to the total charge of the system in the IR limit | q| → 0: Therefore, an IR gluon decouples from a static color-singlet system. Diagrammatically, however, an IR gluon can detect the total charge of the system only when both selfenergy diagrams and potential-energy diagrams are taken into account, as can be seen from Fig. 8. This means that a cancellation takes place between these two types of diagrams, since the IR gluon couples to individual diagrams but decouples from the sum of them.
On the other hand, in analogy with classical electrodynamics, gauge field couples to the total charge of the system in the lowest order [O(r 0 )] of the multipole expansion: which follows from eq. (70). Accordingly, in the pNRQCD Lagrangian (in the static limit), there is no coupling of the singlet field S and the gauge field at the lowest order of the multipole expansion [34]. Hence, the IR cancellation between the self-energy and potential-energy diagrams is explicit at O(r 0 ) in the multipole expansion (OPE) of the total energy of a static QQ pair. 26 Intuitively, IR gluons with wavelengths of order Λ −1 QCD (≫ r) cannot resolve the color charge of each particle, hence they only see the total charge of the system. More accurately, coupling of IR gluons to the system can be expressed by an expansion in r (multipole expansion) for small r, in which the zeroth multipole (=total charge) of the color-singlet QQ pair is zero.
The modern approach (after around 1998) to use the MS mass for the computation of E tot (r) = 2m pole + V S (r) for a heavy quarkonium system can be viewed as follows. The total energy of the system is computed as the sum of (i) the MS masses of Q andQ, (ii) contributions to the self-energies of Q andQ which are not included in the MS masses, and (iii) the potential energy between Q andQ. Contributions of IR gluons with wavelengths larger than r automatically cancel between (ii) and (iii) in this computation [48]. In this way we can eliminate a large part of the IR contributions from the computation of E tot (r).

Appendix B: Derivation of eq. (50)
We show some details of the derivation of eq. (50). The regularized dimensionless potential is given by The integral contour of t is rotated to the positive imaginary axis (t = is). v − can be obtained by setting ǫ → −ǫ and s → −s (or, by taking the complex conjugate). The Borel transform of α V (q) can be expressed in the integral form as 27 We approximate α V (q) by [ α V (q)] N k LL . According to our current knowledge of the RG equation at N k LL, α s (q) diverges at q = q * if the running starts from µ > q * with the initial condition α s (µ) = 1/p > 0. In this case there is a singularity on the real axis of p, due to the singularity of α s (q). At LL (k = 0), the singularity is located at p = b 0 log(µ 2 /q 2 ). At N k LL, the singularity of p is on the real axis for given values of q(> q * ), µ and k, where the relation is given by q = q * (µ, α s (µ) = 1/p; k). We are concerned with the case that q is in the vicinity of q * , where α s (q) ∼ (q −q * ) −1/(k+1) . The singularity of p can be shifted infinitesimally into the upper-half plane (hence, without crossing the contour of p integral) by a shift q → q − iǫ in the vicinity of q * .
After changing the order of the integration, we can integrate over p and s, which transforms α V (q) to α V (q). Then we are left with the q integration with the integral contour deformed into the lower-half plane in the vicinity of q * . The above iǫ-prescription for q specifies how to avoid the singularity of [ α V (q)] N k LL at q = q * compatibly with the deformation of the integral contour of t.