QCD static energy at next-to-next-to-next-to leading-logarithmic accuracy

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 was performed in [6]. The anomalous dimension of the ultrasoft next-to-leading logarithm was calculated in [7] and turned out to be quite large, showing that the resummation of the ultrasoft next-to-leading logarithms should also be addressed. In this paper, we will consider this resummation. It has been argued in [9–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–12], which essentially takes 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 short-distance region. The static potential is a basic ingredient of heavyquarkonium physics [13]. In particular, its perturbative evaluation at higher orders is relevant to describe the topquark 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 (NLO) [14–17]. The past experience with the next-tonext-to-leading order (NLO) 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 NLL calculation of this process. The paper is structured as follows: In the next section, we introduce residual mass terms in potential nonrelativistic QCD (pNRQCD) and summarize the current status of the perturbative calculations for the static potential. In Sec. III, we present and solve the renormalization group (RG) equations for the static pNRQCD Lagrangian, at next-to-leading order, and briefly describe the renormaLeading ultrasoft logarithms contribute to the static energy at order 3þn s ln s for n 0, i.e. at next-to-next-to-leading logarithmic, NLL, order. Next-to-leading ultrasoft logarithms contribute to the static energy at order 4þn s ln s for n 0, i.e. at next-to-next-to-nextto-leading logarithmic, NLL, order. The calculation relies on the next-to-leading order calculation of the chromoelectric correlator done in [8]. PHYSICAL REVIEW D 80, 034016 (2009)


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 1 was performed in [6].The anomalous dimension of the ultrasoft next-to-leading logarithm 2 was calculated in [7]  3 and turned out to be quite large, showing that the resummation of the ultrasoft next-to-leading log-arithms 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 takes 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 short-distance region.
The static potential is a basic ingredient of heavyquarkonium physics [13].In particular, its perturbative evaluation at higher orders is relevant to describe the topquark 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-tonext-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 nonrelativistic QCD (pNRQCD) and summarize the current status of the perturbative calculations for the static potential.In Sec.III, we present and solve the renormalization group (RG) equations for the static pNRQCD Lagrangian, at next-to-leading order, and briefly describe the renorma-lon 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 Sec.V.

II. POTENTIALS AND RESIDUAL MASS TERMS IN PNRQCD
The general form of the dimension six 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 one.
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 (cutoff) are used, one needs to introduce a residual mass term m Q in the Lagrangian [21] with m Q the heavy-quark mass and c 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, defined by 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 with is the ultrasoft factorization scale.The color factors are defined as , 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 nonlogarithmic piece ã3s;o has not been completely calculated yet.The fermionic contributions of ã3s have been presented very recently in [32], where the computation of the n f independent piece is reported to be in progress.A Pade ´estimate of ã3;s gives ã3;s 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; Þ. 6At order r 0 in the multipole expansion, the dimension six operators of pNRQCD do not have an anomalous dimension and, therefore, the RG equations for the coefficients Ã s and Ã o will have the same structure as in the HQET the RG equation for the coefficient of the operator c y c 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, S ¼ T a O a , E is the chromoelectric field, V A and V B are matching coefficients associated with the OðrÞ operators of the pNRQCD Lagrangian, 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 six operators.In particular, this modifies the RG equations for Ã s and Ã o .

III. RENORMALIZATION GROUP
The general structure of the RG 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

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= 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 the T ! 1 limit [4]: T a W h 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 0 Þ 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 that 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 that 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 Ultrasoft contributions to the self-energy of the octet field in a covariant gauge, at the order of interest.Symmetric graphs are understood for (c)-(i).

B. Solution of the RG equations
The RG 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 ), and where The RG equations for Ã s and Ã o can also be obtained from ( 14) by expanding c s;o about V s;o and keeping the terms linear in Ã s;o .They are given by 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 RG 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Þ ¼ 2m Q , as it happens in minimal subtraction-type schemes, then Ã s ðÞ ¼ Ã o ðÞ ¼ 2m 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 2m Q ) and hence Ã s ðÞ and Ã o ðÞ evolve in a nontrivial 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 outcome of the matching calculation only enters through the initial conditions of the RG equations.It is well known that the singlet potential V s calculated in the MS scheme displays a bad behavior as a series in s ð1=rÞ even at small values of r.This bad behavior may be ascribed to renormalon singularities that lie very close to the origin of the Borel plane.In order to treat the renormalon singularity, we shall follow the procedure described in Ref. [34], the so-called RS scheme.Under the RS scheme we understand a class of subtraction schemes that subtract from the perturbative series of V s;o in the MS scheme the nonintegrable piece at u ¼ 1=2 in the Borel transform of the potential, expanded about u ¼ 0 and integrated over u.The whole nonanalytic piece (nonintegrable and integrable) at u ¼ 1=2, expanded about u ¼ 0 and integrated over u, reads (at the scale ) 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 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 nonintegrable 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 vice versa).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 lnr terms, which will be kept from becoming large by choosing ¼ 1=hri ¼ 3:25=r 0 , hri 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 ã3s .

A. Setting the scales and parameters
We choose, as anticipated in the previous section, ¼ 1=hri ¼ 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 r 0 ).The remaining scales and parameters entering in the expressions are chosen as follows: The number of light flavors n f is set to zero.The ultrasoft scale is set to ¼ 2=r 0 .The normalizations of the u ¼ 1=2 renormalon singularities for the singlet and octet potentials, R s;o , are determined using the procedure described in [36]; one obtains [12,37] R s ¼ À1:333 þ 0:499 À 0:338 ¼ À1:172; R o ¼ 0:167 À 0:0624 þ 0:00972 ¼ 0:114: (30) s at the relevant scales is determined according to [38], 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 Pade ´estimate of [33] to get an idea of its expected size (that is, use c Pred 0;n f ¼0 ¼ 313 from Table I of that paper, which corresponds to ã3;s ¼ 114633; the relation between c 0 , V ð3Þ s and ã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 Pade ´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 Pade ´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 Pade ´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 Pade ´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 .

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 dotted-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 Pade ´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 Eq. ( 14) of [7]): K 1 and K 2 are the constants N s Ã and ðN o À N s ÞÃ in Eq. ( 25), while fðr; Þ is 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 use 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 II of [35], and K1 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

FIG. 3 (color online)
. Static potential r 0 V s ðrÞ, in the RS scheme, as a function of r=r o .The dotted blue curve is at tree level, the dotted-dashed magenta curve is at one loop, the dashed brown curve is at two loop plus leading ultrasoft logarithmic resummation, and the long-dashed green curve is at three loop (Pade ´estimate) plus next-to-leading ultrasoft logarithmic resummation.The solid black curve is also at three loop plus next-toleading ultrasoft logarithmic resummation but using c 0 ¼ 250:7, see Sec.IV D.
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 Sec.IVA, we will let c 0 vary by 35% around the Pade ´value c 0 ¼ 313.K 2 will be varied from À5 to 5. Ã MS will be varied according to the range quoted in Ref. [38], Ã 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 AE10C F 5 s =r to the three-loop curve ( AE 10 is intended to be a rough estimate of the size of the four-loop coefficient in the RS scheme).All those , and the one-loop running of s ¼ s ð1=rÞ is used.Impact of the variation of c 0 , K 2 , Ã MS and effect of higher-order terms in the N 3 LL expression for the static energy, respectively, represented as the green bands.c 0 has been varied by 35% around the Pade ´value c 0 ¼ 313.K 2 has been varied from À5 to 5. Ã MS has been varied according to Ã MS r 0 ¼ 0:602ð48Þ [38], and the effect of higher-order terms has been estimated by adding the term AE10C F 5 s =r to the N 3 LL curve.c 0 ¼ 313 and K 2 ¼ 0 have been used as central values.The red points are the lattice data of [35].
variations are shown in Fig. 5 as the green bands.We have used the Pade ´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 K 2 ¼ À1:0465; c 0 ¼ 250:7; (39) the 2 =d:o:f: of the fit is 0.07 (we have 4 degrees of freedom).Those are the best fit values, but to obtain an allowed range of values for c 0 , according to the lattice data, we will analyze the 2 of the N 3 LL curve for different values of the c 0 and K 2 parameters (with Ã MS r 0 ¼ 0:602 and higherorder terms set to zero).We choose the allowed range of values for c 0 by identifying the region of the c 0 -K 2 parameter space where the reduced 2 of the N 3 LL curve is better than that of the N 2 LL curve, and, at the same time, K 2 retains a reasonable power counting value.The reduced 2 of the N 2 LL curve is 3383 (in this case we have 6 degrees of freedom).We will allow values of jK 2 j up to jK 2 j ¼ 2, to be conservatively consistent with our power counting.Figure 6 presents a scatter plot which shows the values of the ratio ð 2 =d:o:f:Þ N 3 LL =ð 2 =d:o:f:Þ N 2 LL with different values of the K 2 and c 0 parameters.In that plot, lighter (darker) colors correspond to higher (lower) values of the ratio.From that we obtain the range (215, 350) for c 0 , as can be read from the figure.To make the values of 2 easier to visualize, we also present, in Fig. 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 Fig. 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 jK 2 j $ 2 s =r $ Ã $ 0:6, in units of r 0 ).This gives us much confidence in the consistency of 220 240 260 280 300 320 340 for different values of the K 2 and c 0 parameters.The color of each point in the c 0 À K 2 plane represents the value of the ratio, lighter colors correspond to higher values and darker colors to lower values.We show values of the ratio up to 1.

(color online)
for different values of c 0 (left) or K 2 (right), the other parameter is fitted.our analysis.Note also that the best fit value for c 0 is smaller than the Pade ´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 Pade ´value c 0 ¼ 313 (see the solid black curve in Fig. 3).
In Fig. 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 Pade ´esti- Comparison of the singlet static energy with lattice data.We plot r o ðE 0 ðrÞ À E 0 ðr min Þ þ E latt 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 dotted-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).Impact of the variation of c 0 , K 2 , Ã MS and effect of higher-order terms in the static energy, 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Þ [38] and the effect of higher-order terms has been estimated by adding the term AE10C 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 Fig. 8.The red points are the lattice data of [35].mate of Ref. [33]).However, if we add to ã3;s the corresponding nonlogarithmic parts coming from US , then we obtain the nonlogarithmic piece of the static energy at three loops (which was denoted as ã3 in [7]).The three-loop coefficient of the static energy ã3 is independent of the scheme used to factorize the ultrasoft contribution, as opposed to ã3;s .At the practical level, the nonlogarithmic part of the ultrasoft contribution just shifts the values of c 0 that we quote above by À1%.Concerning the factorization scale and implementation of the RS scheme dependences, we have checked that varying R s by 30%, R o by 10% [according to the corresponding last known terms of the series in (30)], the 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 Pade ´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 with ã2 ¼ ã2s and ã3 in the range ð1:08; 1:17Þ Â 10 5 .

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 nonperturbative 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 : where the central value corresponds to the best fit of the N 3 LL curve.For those values, an excellent agreement with lattice data is achieved in the region where the weakcoupling calculation is reliable.

FIG. 4 (
FIG.4(color online).Illustration of the hierarchy of scales relevant for the static energy.The scales are: 1=r (dotted red),ðV o À V s Þ (dotted-dashed green), s ðV o À V s Þ (dashed black) and Ã MS ¼ 0:602=r 0 (solid blue).ðV o À V s Þ is taken at tree level, ðV o À V s Þ ¼ N c s =ð2rÞ, and the one-loop running of s ¼ s ð1=rÞ is used.

FIG. 5 (
FIG. 5 (color online).Impact of the variation of c 0 , K 2 , Ã MS and effect of higher-order terms in the N 3 LL expression for the static energy, respectively, represented as the green bands.c 0 has been varied by 35% around the Pade ´value c 0 ¼ 313.K 2 has been varied from À5 to 5. Ã MS has been varied according to Ã MS r 0 ¼ 0:602ð48Þ[38], and the effect of higher-order terms has been estimated by adding the term AE10C F 5 s =r to the N 3 LL curve.c 0 ¼ 313 and K 2 ¼ 0 have been used as central values.The red points are the lattice data of[35].

FIG. 8 (
FIG.8(color online).Comparison of the singlet static energy with lattice data.We plot r o ðE 0 ðrÞ À E 0 ðr min Þ þ E latt 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 dotted-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).

FIG. 9 (
FIG. 9 (color online).Impact of the variation of c 0 , K 2 , Ã MS and effect of higher-order terms in the static energy, 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Þ[38] and the effect of higher-order terms has been estimated by adding the term AE10C 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 Fig.8.The red points are the lattice data of[35].

ã 3 ;
s ¼ 1:11 þ0:06 À0:03 Â 10 5 ; are needed at order 2 s .Strictly speaking the equations above hold for c s;o ¼ V s;o þ Ã s;o , provided that V s;o ) Ã s;o 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.

TABLE I .
Values of the difference c 0;n f À c 0;n f À1 for the exact result and the Pade ´estimates.