The QCD static energy at NNNLL

We compute the static energy of QCD at short distances at next-to-next-to-next-to leading-logarithmic accuracy in terms of the three-loop singlet potential. By comparing our results with lattice data we extract the value of the unknown piece of the three-loop singlet potential.


I. INTRODUCTION
The energy between a static quark and a static antiquark separated at a distance r, which is referred to as the static energy, is a basic object to understand the dynamics of QCD [1].
When calculated at short distances in perturbation theory, the virtual emission of ultrasoft gluons (i.e. gluons with energy and momentum smaller than 1/r that can change the color state of the quark-antiquark pair from singlet to octet) produce infrared divergences. These, after resummation of certain diagrams, induce in the static energy logarithms, ln α s (1/r) [2].
In an effective field theory framework [3,4], the separation of scales in the problem is made explicit: the static energy becomes the sum of a matching coefficient, the static potential, which encodes all contributions from the scale 1/r, and of the contributions coming from ultrasoft gluons, which start at three loops. Logarithms are remnants of the cancellation between infrared divergences of the potential and ultraviolet divergences of the ultrasoft contributions [5]. These logarithms may be potentially large when r is very small, so that α s (1/r) ln α s (1/r) terms in the static energy may need to be resummed. The resummation of the ultrasoft leading logarithms (LLs) 1 was performed in [6]. The anomalous dimension of the ultrasoft next-to-leading logarithm (NLL) 2 was calculated in [7] 3 and turned out to be quite large, showing that the resummation of the ultrasoft NLLs should also be addressed.
In this paper, we will consider this resummation.
It has been argued in [9,10,11,12] that the proper consideration, and cancellation, of the renormalon singularities is crucial to obtain a good convergence of the perturbative series for the static potential in the short-distance region. The detailed analysis of the possible influence of ultrasoft effects in the renormalon structure of the potential will not be presented in this paper. We will just follow the analysis of [9,10,11,12], which essentially take advantage of the fact that a constant term may be added to the potential, a freedom that remains even if ultrasoft effects are taken into account. As we will see later, this seems to be enough to obtain a convergent perturbative series for the static potential in the 1 Leading ultrasoft logarithms contribute to the static energy at order α 3+n s ln n α s for n ≥ 0, i.e. at nextto-next-to-leading logarithmic, N 2 LL, order. 2 Next-to-leading ultrasoft logarithms contribute to the static energy at order α 4+n s ln n α s for n ≥ 0, i.e. at next-to-next-to-next-to-leading logarithmic, N 3 LL, order. 3 The calculation relies on the next-to-leading order (NLO) calculation of the chromoelectric correlator done in [8].
short-distance region.
The static potential is a basic ingredient of heavy-quarkonium physics [13]. In particular, its perturbative evaluation at higher orders is relevant to describe the top-quark pair production process near threshold. This process is expected to allow the extraction of the top quark mass to a high precision, and hence a remarkable effort is being made to calculate it at next-to-next-to-next-to-leading order (N 3 LO) [14,15,16,17]. The past experience with the next-to-next-to-leading order (N 2 LO) results [18] indicates that both renormalon cancellation and the logarithmic resummation [19,20] are necessary for accurate determinations of the position of the pole and of the shape of the cross section respectively. Our results will be relevant for the N 3 LL calculation of this process.
The paper is structured as follows. In the next section, we introduce residual mass terms in potential Non-Relativistic QCD (pNRQCD) and summarize the current status of the perturbative calculations for the static potential. In section III, we present and solve the renormalization group equations for the static pNRQCD Lagrangian, at next-toleading order, and briefly describe the renormalon subtracted scheme that we use. Section IV presents a comparison of our results for the static energy with lattice data, and the numerical extraction of the missing piece of the three-loop static potential. We conclude in section V.

II. POTENTIALS AND RESIDUAL MASS TERMS IN PNRQCD
The general form of the dimension 6 operators in the pNRQCD Lagrangian is where S is the singlet and O a the octet fields. The coefficients c s and c o have dimension 1.
Let us recall that, in order to define Heavy Quark Effective Theory (HQET) beyond perturbation theory, or even in perturbation theory when regularizations with an explicit scale (cut-off) are used, one needs to introduce a residual mass term δm Q in the Lagrangian with m Q the heavy-quark mass and ψ the heavy-quark field. We may associate to δm Q the size of the typical hadronic scale Λ QCD . This residual mass term will be inherited by the pNRQCD Lagrangian. In the paper, we consider the weak-coupling regime of pNRQCD, at leading order in the 1/m Q expansion; in this situation, the residual mass term is absorbed in the coefficients c s and c o above. Therefore, it is useful to split them in a part that is proportional to 1/r, which corresponds to the singlet and octet potential 4 , and a part that is proportional to Λ QCD : where 4 The singlet potential is often referred to as the static potential, a terminology which we also adopt in the paper. Recall that it coincides with the static energy up to two loops but differs from it beyond that order.
µ is the ultrasoft factorization scale. The color factors are defined as C F = (N 2 c − 1)/(2N c ), C A = N c where N c is the number of colors; n f is the number of (massless) flavors; γ E is the Euler constant. The strong coupling constant α s is in the MS scheme. The beta function is defined as where , and explicit expressions of β 2 and β 3 may be found, for instance, in [22,23].
The one-loop coefficientã 1 was calculated in [24,25], the two-loop singlet coefficientã 2,s in [26,27,28,29] and the two-loop octet coefficientã 2,o in [30]. The logarithmic piece of the third-order correction was calculated in [4,5,19,31], whereas the non-logarithmic piecẽ a 3 s,o has not been completely calculated yet. The fermionic contributions ofã 3 s has been presented very recently in [32], where the computation of the n f independent piece is reported to be in progress. A Padé estimate ofã 3,s gives: s (n f = 5) = −20.5 [33]. The double logarithmic coefficient a L2 4 may be obtained from [6,7] and the logarithmic coefficient a L 4 was obtained in [7] 5 . Λ s,o stands for Λ s,o (r, µ) 6 .
At order r 0 in the multipole expansion, the dimension 6 operators of pNRQCD do not have an anomalous dimension and, therefore, the renormalization group equations for the coefficients Λ s and Λ o will have the same structure as in the HQET the renormalization group equation for the coefficient of the operator ψ † ψ has. At next-to-leading order in the multipole expansion, the pNRQCD Lagrangian reads where L light is the part of the Lagrangian involving gluons and light quarks, which coincides with the QCD one, are matching coefficients associated with the O(r) operators of the pNRQCD Lagrangian 5 Only the coefficient of the singlet potential was obtained there, it will be shown later in the paper that it coincides with the coefficient of the octet potential. 6 To simplify the notation, we will often suppress the dependence on r and just write Λ s,o (µ). and the dots stand for higher-order terms in the multipole expansion. Ultrasoft gluons cause transitions between singlet and octet fields and generate an ultrasoft anomalous dimension for the dimension 6 operators. In particular, this modifies the renormalization group (RG) equations for Λ s and Λ o .

III. RENORMALIZATION GROUP
The general structure of the renormalization group equations of pNRQCD in the static case has been discussed in [6], where the complete N 2 LL order was calculated. Here, we will calculate the complete N 3 LL order. It has been proved in [7] that in order to perform the calculation one needs not to consider higher orders in the multipole expansion beyond those already contributing to the N 2 LL calculation. Hence, the structure of the RG equations remains the same as in [6], but the anomalous dimensions need to be calculated to one order more in the ultrasoft loops. The RG equations read where the anomalous dimensions γ s (α s ), γ o (α s ), γ A (α s ) and γ B (α s ) are needed at order α 2 s . Strictly speaking the equations above hold for and one stays at linear order in Λ s,o . If quadratic or cubic terms in Λ s,o are included, additional counterterms in the potential are needed to absorb the ultraviolet divergences of the ultrasoft calculation.

A. Anomalous dimensions
We calculate the anomalous dimensions in a regularization scheme in which the gluons, light quarks and center of mass motion are taken in D dimensions but the potentials in the ultrasoft loops are kept in three dimensions. We renormalize using the MS scheme in (relative) coordinate space.
In this scheme, the anomalous dimension γ s was obtained in [7]. It is −2 α s ×∂Z (1) /∂α s , where Z (1) denotes the coefficient of the 1/ǫ (1/ǫ = 1/ǫ − γ E + ln 4π, D = 4 − 2ǫ) pole of the ultrasoft contribution to the static energy; we have with In a similar way, γ o can be obtained from the 1/ǫ poles of the ultrasoft contribution to the self-energy of the octet field. It is important to recall that, although we may obtain it by matching gauge-dependent Green's functions, the self-energy of the octet field is a gauge invariant quantity in perturbation theory, for the same reason as the pole mass is. At the order we are interested in, it can be obtained from the following expression to be taken in T a W T b stands for a T a and a T b insertion at the time T /2 and −T /2 respectively in the space sides of a rectangular Wilson loop, and φ(t, t ′ ) adj ab is a Wilson line in the adjoint representation. The Feynman diagrams involved in the evaluation of this quantity in a covariant (or Coulomb) gauge do not coincide with the ones needed for γ s (compare Fig. 1 with Fig. 4 of [7]), and would require extra calculations. Fortunately there is an argument which makes the explicit calculation unnecessary. If we take the A 0 = 0 gauge, the number of diagrams to be evaluated collapses to a few (the octet field does not emit gluons anymore), which, in addition, are the same for γ s and γ o . In this gauge, both anomalous dimensions are related by trivial color factors. Since both anomalous dimensions are gauge invariant, it turns out that we can read γ o from the known result for γ s : At order α s , this is confirmed by the explicit calculation of [4].
Let us turn now to the evaluation of γ A and γ B . We use the fact that the µ dependence of V A and V B can be also obtained from the infrared logarithms in the matching calculation between HQET and pNRQCD. There is only one diagram which is infrared divergent in the matching calculation: it is displayed in Fig. 2. However, the divergence turns out to be linear and produces no logarithms. Then these anomalous dimensions remain zero also at next-to-leading order:

B. Solution of the renormalization group equations
The renormalization group equations for V s and V o can be read from (14) using that V s,o ≫ Λ s,o and neglecting the latter. They are given by The solutions of these equations, at the order we are interested in, are V A (µ) = V A (1/r) = 1, (1/r) = 1 + O(α 2 s ), see [7]), and where The renormalization group equations for Λ s and Λ o can also be obtained from (14) where we have already approximated V o (µ) and V s (µ) by V o (1/r) and V s (1/r) (the µ dependence of V o (µ) and V s (µ) enters at N 3 LO, which is beyond the accuracy of (24)). The solutions of the renormalization group equations read, where N s , N o are two arbitrary scale-invariant dimensionless constants and Λ is an arbitrary scale-invariant quantity of dimension one.
The integration constants N s and N o are fixed by the initial conditions, Λ s (1/r) and Λ o (1/r), of the solutions of the RG equations. In turn, the initial conditions are fixed by matching a suitable Green's function in QCD with the corresponding one in pNRQCD.
Note that if at the matching scale Λ s (1/r) = Λ o (1/r) = 2δm Q , as it happens in MS-type schemes, then Λ s (µ) = Λ o (µ) = 2δm Q for any µ. We will see, in the following sections, the convenience to use an RS (renormalon subtracted) scheme. In the RS scheme Λ s (1/r) and Λ o (1/r) are different constants (which also differ from 2δm Q ) and hence Λ s (µ) and Λ o (µ) evolve in a non-trivial way according to the RG equations above.

C. The Renormalon Subtracted scheme
The discussion in the previous sections is independent of the renormalization scheme used for the matching calculation between HQET and static pNRQCD at the scale 1/r.
The coefficients d k are given in terms of the coefficients of the beta function. Since the beta function is known up to four loops only, all d k for k ≥ 3 are unknown; the known terms are d 0 = 1 , Hence, it is practically unfeasible to subtract Eq. (27). This is not a real problem because only subtracting the k = 0 term, which corresponds to the non-integrable piece in the Borel transform, is necessary in order to obtain a series that is Borel summable. In the following, we will subtract all the known terms in Eq. (27), i.e. up to k = 2, as was done in the original proposal of the RS [34].
Furthermore, since the potential is given as an expansion of α s (1/r), in order to achieve a successful renormalon subtraction at every order in α s , it is important to expand α s (ρ) in terms of α s (1/r) (or viceversa). We chose to expand α s (1/r) in terms of α s (ρ) in Eq. (6) instead of doing the reverse in Eq. (27), because the uncertainty in the normalization constants R s,o of the renormalon singularities is then largely absorbed in the arbitrary additive constant needed to compare with lattice data (as it will be described in the next section).
This expansion generates ln rρ terms, which will be kept from becoming large by choosing ρ = 1/ r = 3.25/r 0 , r being the central value of the range where we compare with lattice data and r 0 being the reference scale used in the lattice computation (see the next section).
At the order we are working, we only need to keep terms up to order α 4 s (ρ).

IV. COMPARISON WITH LATTICE RESULTS
In this section we will compare our results with the (n f = 0) lattice data of Ref. [35].
This will allow us to extract a value for the three-loop coefficientã 3 s .

A. Setting the scales and parameters
We choose, as anticipated in the previous section, ρ = 1/ r = 3.25/r 0 (the reference scale r 0 , used in the lattice computation, has a value of about 0.5 fm, see [35] for more details; we will present all our results in units of α s at the relevant scales is determined according to [37], which uses Λ MS r 0 = 0.602(48). We will use the running of α s according to the order we are working at (for instance, for the one-loop curve we use the two-loop running for the order α s term and the one-loop running for the order α 2 s term, and so on). The three-loop coefficientã 3,s , which enters in our N 3 LL results, is unknown. We can use the Padé estimate of [33] to get an idea of its expected size (that is, use c Pred 0,n f =0 = 313 from Table 1 of that paper, which corresponds toã 3,s = 114633; the relation between c 0 , V (3) s and a 3,s is given in Eq. (29d) of [33], for easier reference we reproduce it in the appendix). To estimate the uncertainty that we should associate to that Padé value we can make use of the results of [32]. In Ref. [32], all the fermionic contributions to the three-loop coefficient c 0 are calculated, therefore the difference is known (c 0,n f stands for the coefficient c 0 calculated with n f flavors). Since Ref. [33] presents the Padé estimated values of c 0 for n f = 0, . . . , 6, we can check if those results satisfy the known values for (31) or not. This comparison is presented in Table I 7 . In addition, we can also obtain c 0 for n f = 0, which is the coefficient we need here, from the Padé estimated value of c 0 for n f = 6 and (31): 8 to be compared with the value c 0,n f =0 = 313 presented in [33]. We will take this as an indication that one should assign an uncertainty of around 30 − 40% to the Padé estimate.
As we will see later, though, the lattice data is precise enough to be used to obtain an independent extraction of the value of c 0 .
7 In [32], the coefficient of the d abcd F d abcd F color structure is given numerically, but the limited numerical precision is not yet affecting the numbers presented in Table I. 8 We chose n f = 6, because it is for n f = 6 that the Padé approximation comes closest to the three-loop RG accessible coefficient c 1 , in the notation of [

B. The static potential
Using the choices of scales and parameters described in the previous section, we obtain the singlet potential in the RS scheme represented in Fig. 3 (to alleviate the notation we do not explicitly indicate the dependences on µ and ρ of the different functions in the labels of the plots). In all the plots in this section, the dotted blue curve will be at tree level, the dot-dashed magenta curve will be at one loop, the dashed brown curve will be at two loop plus N 2 LL resummation and the long-dashed green curve will be at three loop (with the Padé estimated value for c 0 , i.e. c 0 = 313) plus N 3 LL resummation. As expected and in sharp contrast to what would happen in an on-shell scheme [12], we see that when we use a threshold scheme that cancels the leading renormalon, like the RS scheme, the perturbative series for the potential exhibits a convergent behavior.

C. The static energy
In order to perform a comparison with lattice data, we have to plot the static energy E 0 as a function of r: where Λ s is given by Eq. (25) and δ US contains the contributions from ultrasoft gluons. V s , K 1 and K 2 have to be understood in the RS scheme, which is where the ρ dependence comes from (at the order we are working, f will not depend on ρ, which will be dropped from it in the following). Also, at the order we are working, the renormalized expression of δ US is only needed at leading order 9 (and can be read from equation (14) of [7]): We are considering the weak-coupling regime in the static limit, defined by the hierarchy of scales given in (3). To have a definite way to organize the terms in Eq. (33) we will also which is compatible with Eq. (3). Figure 4 shows that the scale hierarchy of Eq. (3) as well as the counting above hold for r/r 0 ≤ 0.5.
The lattice data of [35] is presented as the difference between the static energy at distance r and the static energy at a reference scale r c . Therefore, what we actually need to plot, in order to compare with the lattice data, is where r min is the shortest distance at which lattice data is available, r 0 E latt. 0 (r min ) = −1.676 is given in Table 2 of [35], andK 1 is a suitable constant, which can be obtained by imposing the equation to be true at r = r min .
From the counting described above, we see that the last two terms in the rightmost part of Eq. (37) only start contributing at the three-loop level. We also see that, at three-loop level, we only need the function f (r) in Eq. (35) at leading order, i.e.
Therefore, the static energy (that we will use to compare with lattice data) at N 3 LL order is given by Eq. (37) with V s given by (21) (understood in the RS scheme), f given by (38) and δ US given by (34). The constant K 2 will be fixed by a fit to the lattice data below r/r 0 = 0.5. In all the plots, we will always display our results until r/r 0 = 0.5, which is the region where we expect perturbation theory and our hierarchy of scales to be reliable.

Analysis of the uncertainties in the static energy
As we have already mentioned, our N 3 LL results depend on the unknown three-loop coefficient c 0 (in addition to the constant K 2 ). The static energy also suffers from uncertainties due to the Λ MS parameter, used to determine α s , and from uncertainties due to the neglecting of α 5 s /r and higher-order terms. Since we want to use the lattice data to extract c 0 , we need to make sure that the static energy is more sensitive to variations of c 0 (and K 2 , which will be also fitted to the data) than to the errors due to Λ MS and the higher-order terms.
According to the discussion in section IV A, we will let c 0 vary by 35% around the Padé value c 0 = 313. K 2 will be varied from -5 to 5. Λ MS will be varied according to the range quoted in Ref. [37], Λ MS r 0 = 0.602(48). Finally we assess the impact of neglecting the four-loop, α 5 s /r, terms by simply adding the term ±10 C F α 5 s /r to the three-loop curve (±10 is intended to be a rough estimate of the size of the four-loop coefficient in the RS scheme).
All those variations are shown in figure 5 as the green bands. We have used the Padé value c 0 = 313 for the plots where we vary K 2 , Λ MS and the higher-order terms. K 2 = 0 is taken for the plots where we vary c 0 , Λ MS and the higher-order terms. We can clearly see that the variations due to c 0 and K 2 produce larger bands than those due to the uncertainties in Λ MS and higher-order terms. Therefore, it will make sense to fit c 0 and K 2 to the lattice data.

D. Lattice comparison and extraction of c 0
We are now ready to compare with lattice data. As we have seen, the N 3 LL expression for the static energy depends on two unknown parameters, K 2 and c 0 , which we will obtain from a fit to the lattice data points below r/r 0 = 0.5. The result of this two parameters fit (Λ MS r 0 = 0.602 and higher-order terms are set to zero) gives the values To make the values of χ 2 easier to visualize, we also present, in figure 7, separate scans over values of c 0 and K 2 , i.e. we scan over c 0 or K 2 , fit the other parameter for each point, and present the ratio of χ 2 s as a function of c 0 or K 2 .
In figure 8, we show the N 3 LL curve, with the best fit values for K 2 and c 0 given in (39), together with the lattice data. We also show in the plot the tree-level, one-loop and N 2 LL curves for the static energy. Several comments are in order. First, we can see that the agreement with lattice is improved when we go from tree level, to one loop, to N 2 LL. We also note that the N 3 LL curve describes very well the data. The fit of the N 3 LL curve to the lattice data was not constrained to give a value of K 2 compatible with the counting, but this turns out to be the case (recall that our power counting requires |K 2 | ∼ α 2 s /r ∼ Λ ∼ 0.6, in units of r 0 ). This gives us much confidence in the consistency of our analysis. Note  also that the best fit value for c 0 is smaller than the Padé estimate of [33], but in better agreement with (32). We would like to emphasize that the values of c 0 in Eqs. (32) and (39) are obtained by completely independent procedures. Finally, let us also note that, when we use the value of c 0 in Eq. (39), the convergence of the perturbative series for the potential seems to be slightly improved, with respect to using the Padé value c 0 =313 (see the solid black curve in figure 3).
In figure 9 we present the bands obtained by varying c 0 and K 2 according to the ranges described above, and the bands induced by the variations in Λ MS and higher-order terms.
All the bands are obtained by keeping the rest of the parameters at their central or best fit values.
Let us comment at this point on the scheme, factorization scale and implementation of the RS scheme dependences of the extracted value of c 0 (or equivalentlyã 3,s ). Concerning the scheme dependence, recall that the scheme of Ref. [7] has been used to factorize the ultrasoft contributions, which differs from standard MS scheme (we note that ultrasoft contributions were ignored in the Padé estimate of Ref. [33]). However, if we add toã 3,s the corresponding non-logarithmic parts coming from δ US , then we obtain the non-logarithmic piece of the static energy at three loops (which was denoted asã 3 in We plot 0 (r min ) as a function of r/r 0 and the lattice data of [35] (red points).
The dotted blue curve is at tree level, the dot-dashed magenta curve is at one loop, the dashed brown curve is at two-loop plus leading ultrasoft logarithmic resummation and the solid black curve is at three-loop plus next-to-leading ultrasoft logarithmic resummation, using the best fit values of Eq. (39).
scales ρ and µ by 10% and implementing the RS subtraction with just the d 0 term in (27) shifts the range for c 0 by a maximum of 7%.
We conclude by noting that both the best fit value from the lattice comparison and Eq.
(32) indicate a value of c 0 lower than the Padé estimate c 0 = 313. The computation of this three-loop coefficient is reported to be in progress [32]. For the sake of comparison, the static energy at N 3 LO is given by E 0 (r) = − C F α s (1/r) r 1 +ã 1 α s (1/r) 4π +ã 2 α s (1/r) 4π 2 + 16 π 2 3 C 3 A ln C A α s (1/r) 2 +ã 3 α s (1/r) 4π withã 2 =ã 2 s andã 3 in the range (1.08, 1.17) × 10 5 . respectively, represented as the gray bands. c 0 has been varied in the interval (215, 350), K 2 in the interval (-2, 2), Λ MS has been varied according to Λ MS r 0 = 0.602(48) [37] and the effect of higher-order terms has been estimated by adding the term ±10 C F α 5 s /r to the N 3 LL curve. All the bands are obtained by fixing the other parameters at their central or best fit values. The curves are the same as in figure 8. The red points are the lattice data of [35].

V. CONCLUSIONS
We have calculated the QCD static energy at short distances at N 3 LL accuracy, in terms of the three-loop singlet potential, whose coefficientã 3,s for n f = 0 is the only missing ingredient in our calculation. It is remarkable that such a higher-order calculation can be carried out analytically, which shows, once more, what invaluable tools effective field theories provide for higher-order calculations. The static energy at this order turns out to depend on two arbitrary constants, rather than one, which encode non-perturbative effects that are competing with the weak-coupling calculation at the considered accuracy. We have used the lattice data of Ref. [35] to extract the value of the unknown piece of the three-loop singlet potential. Our analysis indicates the following value ofã 3,s : a 3,s = 1.11 +0.06 −0.03 × 10 5 , we have where for n f = 0