S-wave heavy quarkonium spectrum with next-to-next-to-next-to-leading logarithmic accuracy

We obtain the Potential NRQCD Lagrangian relevant for S-wave states with next-next-to-next-to-leading logarithmic (NNNLL) accuracy. We compute the heavy quarkonium mass of spin-averaged $l= 0$ (angular momentum) states, with otherwise arbitrary quantum numbers, with NNNLL accuracy. These results are complete up to a missing contribution of the two-loop soft running.


Introduction
High order perturbative computations in heavy quarkonium require the use of effective field theories (EFTs), as they efficiently deal with the different scales of the system. One such EFT is Potential NRQCD (pNRQCD) [1,2] (for reviews see [3,4]). The key ingredient of the EFT is, obviously, its Lagrangian. At present the pNRQCD Lagrangian is known with next-to-next-to-next-to-leading order (NNNLO) accuracy [5] (for the nonequal mass case see [6]).
One of the major advantages of using EFTs is that it facilitates the systematic resummation of the large logarithms generated by the ratios of the different scales of the problem. For the case at hand we are talking of • the hard scale (m, the heavy quark mass), • the soft scale (mv, the inverse Bohr radius of the problem), • and the ultrasoft scale (mv 2 , the typical binding energy of the system).
At present, the pNRQCD Lagrangian is known with next-to-next-to-next-to-leading log (N 3 LL) precision as far as P-wave states is concerned [7]. For S-wave observables the present precision is NNLL [8]. The missing link to obtain the complete N 3 LL pNRQCD Lagrangian is the N 3 LL running of the delta(-like) potentials 1 . For the spin-dependent case, such precision for the running has already been achieved in [9,10]. Therefore, what is left is to obtain the N 3 LL result for the spin-independent delta potential. This is an extremely challenging computation. We undertake this task in this paper.
The new results we obtain in this paper are the following: • We compute the α/m 4 and the α 2 /m 3 spin-independent potentials. These potentials are finite. The expectation value of them produce energy shifts of order mα 6 , which contribute to the heavy quarkonium mass at N 3 LO. Nevertheless, since some expectation values are divergent, some of these energy shifts are logarithmic enhanced, i.e. of order O(mα 6 ln( ν mα )). Such corrections contribute to the heavy quarkonium mass at N 3 LL. This divergence, and the associated factorization scale ν, gets canceled by the corresponding divergence in the spin-independent delta potential. By incorporating the HQET Wilson coefficients with LL accuracy 2 in the α/m 4 and α 2 /m 3 spin-independent potentials, the divergent structure of their expectation value (tantamount to compute potential loops) determines the piece associated to these potentials of the renormalization group (RG) equation of the spin-independent delta potential with N 3 LL precision.
• We compute the (soft-)α 3 /m 2 contribution to the spin-independent delta-like potential proportional to [c (1) F ] 2 , [c (2) F ] 2 ,c (1)hl 1 andc (2)hl 1 . Unlike before, this potential is divergent. Therefore, for future use, we also give the renormalized expression. The divergent pieces produce corrections of O(mα 6 ln( ν mα )) (i.e. of order N 3 LL). From these divergences we generate the (soft) RG equation of the spin-independent delta potential and resum logarithms with N 3 LL precision. In order to reach this accuracy, we need the NLL running of the 1/m 2 HQET Wilson coefficients. For c F this is known [15,16] but not forc hl 1 (the associated missing term is of O(T f n f mα 6 ln(1/α)) and is expected to be quite small. Its computation will be carried out elsewhere). The possible mixing between the (soft-)α 3 /m 2 and the α 2 /m 3 spin-independent potential computed in this paper is also quantified.
The computation of the (soft-)α 3 /m 2 contribution to the spin-independent delta-like potential, proportional to other NRQCD Wilson coefficients, like [c (1) k ] 2 , [c (2) k ] 2 , and c (1) k c (2) k , will be performed in a separated paper. The associated contribution to the running is expected to be small in comparison with the total running of the heavy quarkonium potential. We will estimate its size using the result of the running of the already computed soft contribution.
• The N 3 LL ultrasoft running of the static, 1/m and 1/m 2 potential was originally computed in [17,18,19] (see also [20,21]). This is enough for P-wave analyses [7], where such corrections produce a N 3 LL shift to the energy. Nevertheless, it is not so for S-wave states, as already noted in [9,10] for the case of the hyperfine splitting. The reason is the generation of singular potentials through divergent ultrasoft loops. We revisit it in Sec. 5.2 and incorporate the missing contributions needed to have the complete ultrasoft-potential running that produces N 3 LL shifts to the energy.
• Finally, we compute the complete (potential) RG equation of the delta potential with N 3 LL accuracy (the first nonzero contribution). Solving this equation we obtain the complete N 3 LL running of the delta potential. This allows us to obtain the S-wave mass with N 3 LL accuracy. It is also one of the missing blocks to obtain the complete NNLL RG improved expression of the Wilson coefficient of the electromagnetic current. This, indeed, is what is needed to achieve NNLL precision for non-relativistic sum rules and t-t production near threshold. As the spin-dependent (and l = 0) contribution has already been computed in earlier papers [9,10,7], we only consider here energy averages of S-wave states where the spin-dependent contributions vanish, and only include terms relevant for the N 3 LL S-wave spin-average energy.
Throughout this paper we work in the MS renormalization scheme, where bare and renormalized coupling are related as (D = 4 + 2 ) where T F n f , n f is the number of dynamical (active) quarks and α = g 2 ν 2 /(4π). This definition is slightly different from the one used, for instance, in [22]. In the following we will only distinguish between the bare coupling g B and the MS renormalized coupling g when necessary. The running of α is governed by the β function defined through α(ν) has n f active light flavours and we define z = α(ν) α(ν h ) 1 β 0 1 − 1/(2π)α(ν h ) ln( ν ν h ). Note that with the precision achieved in this paper we need in some cases the two-loop running of the coupling when solving the RG equations.

NRQCD Lagrangian: 1/m 3 and beyond
Instrumental in the determination of the Wilson coefficients of the pNRQCD Lagrangian is the determination of the Wilson coefficients of the Lagrangian of the EFT named NRQCD [23,24]. We first need to assess which NRQCD operators we have to include in our analysis. We will include light fermions, which we will take to be massless.
The HQET 1/m 3 Lagrangian can be found in [25], and including light fermions, though in a different basis, in [26]. Here we use the basis and notation from [14], which also includes light fermions. In [14] one can find the resummed expressions of the Wilson coefficients with LL accuracy for the spin-independent operators. For the spin-dependent 1/m 3 operators, not relevant for this work, the LL running can be found in [27,28]. Note that there are no pure gluonic operators of dimension seven.
To obtain the complete 1/m 3 NRQCD Lagrangian, one also has to consider possible dimension-seven four heavy-fermion operators. There are no such operators, as mentioned in [29]. At O(1/m 4 ), we do not need the complete Lagrangian for the purposes of this paper. For the heavy-quark bilinear sector the complete set of operators was written for the case of QED in [30] and for QCD in [31] (in the last case without light fermions). Of those we can neglect most, we do not need the spin-dependent 1/m 4 operators, nor terms proportional to a single B, nor terms with two (either B or E) terms. The reason is that we only need 1/m 4 tree level potentials. Therefore, we can take all relevant operators from the QED case. Following the notation of [30], the possible relevant operators are and similarly for the antiquark. The dots stand for terms that one can trivially see that do not contribute to the S-wave spin-independent spectrum at NNNLL, either because involved the emission of two gluons or because they are spin-dependent. In principle we need three new coefficients. Nevertheless, we will see later that only c X1 contributes to the running of the spin-independent delta potential. Still, we will compute any tree level potential proportional to c X1 , c X2 and c X3 . The fact that we need c X1 , one of the Wilson coefficients of the 1/m 4 heavy quark bilinear Lagrangian, could make it necessary to consider the Wilson coefficients of the 1/m 4 heavy-light operators as well [light-light operators are subleading for the same reason they are at O(1/m 3 )], as they may enter through RG mixing. Fortunately, c X1 can be determined by reparameterization invariance, which gives us the following relation [30]: (where one should take Z = 1 for QCD). Note that it depends on c D , so indeed c (i) X1 is gauge dependent. Nevertheless, we will see later that it always combines with c M to produce gauge invariant combinations. This indeed is a nontrivial check of the computation. Note also that the above coefficient has an Abelian term, so it can be checked with QED computations.
Finally, we consider the heavy four-fermion sector of the 1/m 4 Lagrangian. They generate local or quasi-local potentials, which do not produce divergent potential loops. The same happens for the potentials generated by c X2 , c X3 . Therefore, in both cases, such potentials do not generate contributions to the heavy quarkonium mass at N 3 LL, and we can neglect them.
Out of this discussion, we conclude that we have the LL running of all necessary Wilson coefficients of the 1/m 4 NRQCD Lagrangian operators.

pNRQCD Lagrangian
Integrating out the soft modes in NRQCD we end up with the EFT named pNRQCD. The most general pNRQCD Lagrangian compatible with the symmetries of QCD that can be constructed with a singlet and an octet (quarkonium) field, as well as an ultrasoft gluon field to NLO in the multipole expansion has the form [1,2] h o (r, p, P R , S 1 , , O], P R = −i∇ R for the singlet, P R = −iD R for the octet (where the covariant derivative is in the adjoint representation), p = −i∇ r , and M = m 1 + m 2 . We adopt the color normalization for the singlet field S(r, R, t) and the octet field O a (r, R, t). Here and throughout this paper we denote the quark-antiquark distance vector by r, the center-of-mass position of the quark-antiquark system by R, and the time by t.
Both h s and the potential V s are operators acting on the Hilbert space of a heavy quark-antiquark system in the singlet configuration. 3 V s (and V o ) can be Taylor expanded in powers of 1/m (up to logarithms). At low orders we have where, [17,18]. The N 3 LL result for the 1/m and 1/m 2 momentum dependent potential is also known in different matching schemes [19,32,7]: on-shell, off-shell (Coulomb, Feynman) and Wilson. In terms of the original definitions used in these papers they read (in four dimensions) The spin-dependent and momentum-dependent potentials are also known with N 3 LL precision [7]. We use the following definitions in this paper (again we refer to [7]): More delicate are V r , as their running is sensitive to potential loops, which are more efficiently computed in momentum space. Therefore, it is more convenient to work with the potential in momentum space, which is defined in the following way: Then the potential reads where the (Wilson) coefficientsD generically stand for the Fourier transform of the original Wilson coefficients in position space D. For them (and for αṼ ) we use the power counting LL/LO for the first nonvanishing correction, and so on.
is indeed known with the required N 3 LL accuracy [9,10] (one should be careful when comparing though, as there is a change in the basis of potentials used there, compared with the one we use here). In terms ofD where and we neglect higher order logarithms (as they are subleading). Finally we consider V r . In terms ofD Unlike all the other potentials, we do not know V (2) r with N 3 LL expression (though the N 2 LL expression is known [8]). This leads us to the main purpose of this paper: the computation of V r with N 3 LL accuracy. This is equivalent to obtaining the NLL expression ofD (2) d . This will require the use of the other Wilson coefficients to one order less: LL. Indeed in Eq. (18) we have already approximated the Fourier transform of V (2) L 2 by its N 2 LL expression (otherwise the momentum dependence is more complicated).
At LL the Wilson coefficients are equal in position and momentum space. We only explicitly display those that we will need later. For the static potential we would have at LL that α V = αṼ = α. For the rest, we show the results in the off-shell Coulomb (which are equal to the Feynman at this order) and on-shell matching schemes, except for D (2) LS i , which we do not need for S-wave: We now turn toD (2) d . ExpandingD (2) d (k, ν) in powers of ln k, we obtaiñ where we have made explicit the dependence on the different factorization scales. So far we have not made explicit the dependence on ν h ∼ m. Nevertheless, it will play an important role later, when solving the RG equations. Therefore, in the following, we use the notationD can be written in several ways: as a sum of the LL (D (2),LL d (ν h ; ν)) term and the NLL (δD (2),N LL d (ν h ; ν)) correction, or as the sum of the initial condition (D d (ν h )) at the hard scale and the running contribution (δD This Wilson coefficient may depend on the matching scheme. Here we mainly consider the off-shell Coulomb gauge matching scheme. Still, for later discussion, we also give expressions in the on-shell matching scheme (see [6] for more details). The LL running is known [8]: is a gauge invariant combination of NRQCD Wilson coefficients, for which its LL running can be found in [8]. In order to visualize the relative importance of the NLL corrections compared with the LL term, we plot the later in Fig From the LL result (using the ν s independence of the potential at LO) one obtains k d dkD This term contributes to the N 3 LL energy shift of the spectrum. 4 Unlike in the other plots, we use here the two-loop running for α. The effect is small.
Since we know the NLO expression ofD (2) d , we can determine the initial matching condition. It reads c D , and the four-fermion Wilson coefficients d ss and d vs , were computed at one loop in [25] and [33] respectively, where one can find the explicit expressions.
At the order we are working δD can be split into pieces. Thus, the NLL approximation for the Wilson coefficient is given by the sum where the second line is zero when 142 for n f =4, 3, and 0 respectively. We nicely observe that these numbers generate small corrections to the leading order results. At present the NLL running is only known for the ultrasoft term [19]: where a 1 = 31/9C A − 20T F n f /9. We show the size of this correction in Fig. 2. Note that the ultrasoft contribution to the delta potential vanishes in the large N c limit (it is 1/N 2 c suppressed). Nevertheless, it quickly becomes big at relatively small scales because the overall coefficient is large and the ultrasoft scale quickly becomes small. Finally, note also that part of the ultrasoft correction (proportional to ln k) is included in Eq. (32). The missing terms to obtain the complete NLL running ofD (2) d are then δD and δD (2),N LL d,p (ν h ; ν). For δD (2),N LL d,s (ν h ; ν) we need the two-loop soft computation ofD and the associated soft RG equation, which we partially obtain in Secs. 4.3 and 5.1, respectively. We also discuss the mixing with higher order 1/m potentials in Sec. 4.4. For δD (2),N LL d,p (ν h ; ν) we need to determine and solve the potential RG equation. This requires first the matching between NRQCD and pNRQCD to higher orders in 1/m, which we do in Secs. 4.1 and 4.2, an extra (ultrasoft associated) running, which we obtain in Sec. 5.2, and obtaining the potential RG equation, which we do in Sec. 5.3.

NRQCD-pNRQCD matching, spin-independent
In this section we compute the potentials for which their expectation values produce corrections to the spectrum of O(mα 6 ). This means the O(α/m 4 ), O(α 2 /m 3 ) and O(α 3 /m 2 ) potentials. Of them we mostly care about those that produce logarithmic enhanced contributions to the spectrum. Therefore, in particular, we do not need to consider the p 6 /m 5 correction to the kinetic term, since it does not give an ultraviolet divergent correction. The O(α/m 4 ) and O(α 2 /m 3 ) potentials are finite. Some of them can be traced back from the QED computation. We mainly compare with [34] (but one could also look into [35] for the equal mass case). Logarithmic enhanced corrections are produced by the divergences generated when inserting these potentials in potential loops. On the other hand the logarithmic enhanced contribution to the spectrum due to the O(α 3 /m 2 ) is not generated by potential loops but by the divergent structure of the potential itself, which we then refer to as soft running. This case will be discussed separately in Sec. 5.1.
The spin-dependent case was computed in [9,10]. Explicit expressions for the potentials can be found in the Appendix of [36]. They produced corrections to the hyperfine splitting (but not to the fine splittings, as shown in [7]).

O(α/m 4 ) potential
From a tree level computation (see the first diagram in Fig. 3) we obtain the complete (spin-independent) α/m 4 potentials in momentum space: In this result we have already used the (full) equations of motion replacing [37] Such k 2 0 terms are generated by Taylor expanding in powers of the energy k 0 the denominator of the transverse gluon propagator.
Not all terms in Eq. (37) contribute to the NLL running of the delta potential. The ones that are local (or pseudo-local) do not contribute, as they do not produce potential loop divergences, since the expectation values of these potentials are proportional to |ψ(0)| 2 and/or (analytic) derivatives of it (kind of ∇ 2 |ψ(0)| 2 ), which are finite. This happens for instance for the potentials proportional to c 2 D , c X2 and c X3 . It is also this fact that allows us to neglect 1/m 4 potentials generated by dimension eight four-heavy fermion operators of the NRQCD Lagrangian.
As we have incorporated the LL running of the HQET Wilson coefficients, these potentials are already RG improved.
Note that with trivial modifications these potentials are also valid for QED.

O(α 2 /m 3 ) potential
We now compute the complete set of the O(α 2 /m 3 ) spin-independent potentials. We show the relevant topologies that contribute to the α 2 /m 3 potential in Fig. 3. By properly The (d) type diagrams in Fig. 3 The (f) type diagrams in Fig. 3 The rest of topologies ((g), (h), (i), (j)) do not contribute. Note that those topologies include, in particular, the one-loop diagrams proportional to c hl i or d hl i , as they may produce ∼ α 2 /m 3 potentials. We find that such contributions vanish.
As we have incorporated the LL running of the HQET Wilson coefficients, these potentials are already RG improved.
Note that with trivial modifications these potentials are also valid for QED.
In this section we perform a partial computation of the O(α 3 /m 2 ) soft contribution to the V r potential. The contributions we compute here are those proportional to the HQET Wilson coefficientsc Using the notation of [6],  (12 + 17)) + 45) + 15) csc 2 (π ) Γ 2 + 5 2 + 24π 5/2 ( ( ( (4 ( + 12) + 127) + 130) + 65) + 15) csc(π ) csc(2π ) (4 ( + 2) + 3)Γ 2 + 5 2 + 12π 4 (2 − 1) sec 2 (π ) Γ 2 ( + 1) With obvious changes the same result is obtained forD (0,2) r,3 . It is worth emphasizing that this expression vanishes in pure QED. A non trivial check of this result is that c D and c hl 1 appear in the gauge invariant combinationc hl 1 = c D + c hl 1 . Another nontrivial check is that the counterterm is independent of k and that the 1/ 2 terms comply with the constraints from RG. This computation has been done in the Feynman gauge (with a general gauge parameter ξ) in the kinematic configuration p = k and p = 0. We also set the external energy to zero. Not setting it to zero produces subleading corrections (we recall that the one-loop computation of this contribution has no energy dependence [6]). The result is shown to be independent of the gauge fixing parameter ξ.
For future computations, it is useful to explain the convention we have taken for the D-dimensional spin matrices. For the c (i) F vertex we typically take a covariant notation ∼ σ µν (see for instance [38]) and project to the particle to single out the spin-independent part: ∼ T r[(I + γ 0 )/2(· · · · · · )(I + γ 0 )/2]. At one loop this procedure gives the same result than using Pauli matrices with the conventions used in [6].
Though not directly relevant for this work, we also give the MS renormalized expression of the bare potential computed above. It will be of relevance for future computations of the spectrum (and decays) at N 4 LO. The result reads (α = α(ν)) where we have also included the O(α 2 ) term. Note that this contribution does not mix with V L . Therefore, it really corresponds to the contributions proportional to c Finally, note that the missing part of the soft term should carefully be computed in a way consistent with the scheme we have used for the rest of the computation, in particular of the α 2 /m 3 potential, as a strong mixing (if using field redefinitions) of the terms proportional to c 2 k is expected.

Equations of motion
Some of the potentials we have obtained in Sec. 4.2 are energy dependent. If we want to eliminate such energy dependence, and write an energy independent potential, this could be achieved by using field redefinitions. At the order we are working it is enough to use the full equation of motion (at leading order), which includes the Coulomb potential. Let us see how it works. We first consider Eq. (39). It depends on the total energy of the heavy quarkonium and does not contribute to the running of the delta potential. We next consider Eq. (51), which is the only energy dependent potential proportional to c (i)2 F . Such potential is generated by the following interaction Lagrangian (56) For this Lagrangian one can use the equations of motion (V C (x) = −C F α/|x|): and similarly for the other fields. We then obtain where the dots stand for the analogous contribution for the antiparticle. The first term in Eq. (58) yields the potential we had obtained after using the free on-shell equations of motion in Eq. (51). It reads The second term is a six-fermion field term. After contracting two of them, a new α 3 /m 2 potential is generated (here we only care about the divergent part). It reads It is worth mentioning that this contribution has a different color structure as those (purely soft) computed in Sec. 4.3, and that they are π 2 enhanced compared to those also. Therefore, one could expect them to be more important than the strict pure-soft contribution.
Remarkably enough, we will see later that the contributions from Eq. (59) and Eq. (60) to the running of the delta potential cancel each other in the equal mass case (but not for different masses). This was to be expected, since in the equal mass case, the potential can be written in terms of the total energy of the heavy quarkonium, which does not produce divergences that should be absorbed in the delta potential.
It is worth mentioning that this exhausts all possible c (i)2 F structures that can be generated. To be sure of this statement, we have to check that the result does not depend on the gauge. Therefore, we have redone the diagrams proportional to c 1loop andṼ (f, 9) 1loop ) in the Feynman gauge and found the same result.
The other potentials that are dependent on the energy are proportional to c 2 k . As before, these contributions will mix with the α 3 /m 2 pure-soft contribution proportional to c 2 k , which we have not computed anyhow. Therefore, in this paper, we only include the explicit contribution generated using the free equations of motion and postpone the incorporation of the other contribution to have the full result. The contributions we explicitly include in this paper then read:

5D
(2) d NLL running We now compute the NLL soft and potential running ofD  . In practice, such computation can be understood as getting the NLL soft running of d ss + C Fdvs (see Eq. (29) or Eq. (30)). It reads

Soft running
The O(α 3 ) stands for terms proportional to NRQCD Wilson coefficients different from c . This equation is meant to represent the pure soft running of the NRQCD Wilson coefficients. It does not give the full running ofD (2) d , as one should also include the potential and ultrasoft running. We fix the initial matching condition to zero, since we only need the initial matching condition of the total potential, which can be determined in the final step, when combining the different contributions.
The strict NLL contribution to the solution of this equation reads (the LL is already included in Eq. (29)) We do not aim in this paper to give a full fledged phenomenological analysis. Still, we  To this contribution one should also add the contributions generated by the new α 3 /m 2 potentials that appear after using the full equations of motion. Of those we only This generates a new contribution to the soft RG equation: Its solution reads We then show the size of this new contribution in Fig. 5.
Finally, let us note that the c 2 k terms can also mix with α 2 /m 3 potentials through field redefinitions, see the discussion in the Appendix. Therefore, this contribution could be different for other matching schemes.

Ultrasoft running
To obtain the complete potential RG equation, we also need an extra potential divergence that is generated by ultrasoft divergences. This term was already computed in [36], and applied to the spin-dependent case. Here, we give the full term, which contributes to both, the spin-dependent and spin-independent term. It is generated by the following diagram which produces the following ultrasoft RG equation or alternatively (but equivalent at this order) Using that the LL running of c F is independent of the masses (we take the initial matching condition to be ν h for both heavy quarks), its solution reads or where (we use the same notation as in [36]) V S 2 ,1/r 3 is singular and will contribute to the potential running ofD (2) d . We now have all the necessary preliminary ingredients to obtain the complete potential RG equation. The next step is to compute all potential loops that produce ultraviolet divergences that get absorbed inD (2) d and are at most of O(α 3 ). Since the delta-like potential is of O(1/m 2 ) we must construct potential loop diagrams of O(α n /m 2 ) with n ≤ 3 describing the interaction between the two heavy quarks in the bound state through several potentials. The first non-vanishing contribution to the potential running is indeed of O(α 3 /m 2 ). To construct such potential loop diagrams, we must consider the power of α and m of each potential and take into account that each propagator adds an extra power   of the mass m in the numerator. We summarize all kind of diagrams that contribute to the NLL potential running ofD (2) d , in Figs. 6, 7, 8 and 9. The ultraviolet divergences arising in such diagrams must be absorbed in the 1/m 2 potentials. However, after the computation, we observe that all divergences are only absorbed by the delta-like potential. It is important to mention that the iteration of two or more spin-dependent potentials can give a contribution toD (2) d , associated to a spin-independent potential. The relevant diagrams are shown in Figs. [6][7][8][9], where V C is the tree level, O(α), Coulomb potential, V α r /m s is the O(α r /m s ) potential and V 1/m 3 corresponds to the first relativistic correction to the kinetic energy, and it is proportional to c 4 .

Potential running
It is interesting to discuss in more detail which, of the novel α 2 /m 3 potentials computed in Secs. 4.2 and 4.4 (we remind that here we use the potentials after using the (free) equations of motion, i.e. the expressions in Sec. 4.4 for the energy dependent potentials), contribute to the running ofD (2) d . The potentials in Eqs. (39)(40) do not contribute to the running ofD (2) d . Equation (39) does not because it is proportional to a total derivative, whereas Eq. (40) does not because of the following argument: the only possible potential loop that can be constructed with an O(α 2 /m 3 ) potential is the iteration of it with a Coulomb potential. As a consequence, the α 2 /m 3 potential is always applied to an external momentum. When the high loop momentum limit is taken in the integral in order to find the ultraviolet pole, all these external momenta vanish and all the terms become proportional to |k| 1+2 . After doing so and summing all the terms the overall coefficient is zero, explaining the fact that they do not contribute. This argument also applies toṼ (e,1) andṼ (f,i) (with i = 1 to 6). On the other handṼ (e,2) andṼ (f,7/8/9) do contribute to the running. Note thatṼ (f, 8) andṼ (f,9) were originally dependent on the energy.
Diagrams with V 1/m 3 in the extremes of a potential loop, i.e. acting over a external momentum have not been drawn because they do not produce any ultraviolet divergence. Similarly, diagrams with V 1/m 5 do not produce ultraviolet divergences. One can then easily convince himself that there are no diagrams with five potential loops or more that can contribute to the O(α 3 ) anomalous dimension ofD d . Therefore, the above discussion exhausts all possible contributions to the O(α 3 ) anomalous dimension ofD d , and the potential RG equation finally reads ν dD The first five lines are generated by potential loops with α 2 /m, α/m 2 and the p 4 /m 3 correction to the kinetic energy (besides the iteration of the Coulomb potential, accounted for by α V ). The 6th line is the term generated by the potential computed in Sec. 5.2.
The last four lines are generated by potential loops with the α 2 /m 3 and α/m 4 potentials (besides the iteration of the Coulomb potential, accounted for by α V ). Note that, for simplicity, we have already used c (i) k = 1 [39] in the terms that do not have NRQCD Wilson coefficients in the above expression. A part of this equation was already computed in [40]. Also, several of these terms (for QED) can be checked with the computations in [34].
It is interesting to see that there is a matching scheme dependence of the individual α 2 /m 3 and α/m 4 potentials that cancels out in the sum. In the above expression the coefficients c A 2 , c D , c M , c X1 appear (note that the last two coefficients are dependent on c D due to reparameterization invariance). They are gauge dependent quantities. Such gauge dependence should vanish in the final result. Indeed it does. This is actually a strong check of the computation. In Eq. (73) we can approximate α V = α (everything is needed with LL accuracy). Then we can show that everything can be written in terms ofc A 2 , which is gauge independent (it is an observable in the low energy limit of the Compton scattering, see the discussion in [14]), and the explicit dependence in c D , c M , c X1 , c A2 disappears. The resulting expression reads ν dD 1 + 12D .
From this result one may think that c A3 and c A4 contribute to the Abelian case. Nevertheless, the LO matching condition is zero for these Wilson coefficients, and all the running vanishes in the Abelian limit. Therefore, there is no contradiction with the pure QED case.
In order to solve Eq. (74), we need to introduce the D s, the Wilson coefficients of the potentials. The necessary expressions can be found in Sec. 3. Note that in those expressions we have already correlated the ultrasoft factorization scale ν us with ν and ν h using ν us = ν 2 /ν h . We also do so in Eq. (72) (where we also set 1/r = ν, consistent with the precision of our computation). This correlation of scales was first introduced and motivated in [41]. For n f = 3 or 4 it is not possible to get an analytic result for the solution of the RG equation, more specifically for the coefficients multiplying the different z functions (note that this comes back to the fact that the polarizability Wilson coefficients c A1 , c A2 , ..., cannot be computed analytically). On top of that the resulting expressions are too long. Therefore, we only explicitly show the analytic result with n f = 0. It reads Finally, in Fig. 10, we give the numerical evaluation δD for different values of n f . The contribution is sizable.

Potential running, spin-dependent delta potential
Even though not relevant for this paper, we profit to present the potential RG equation ofD (2) S 2 in the basis we use in this paper, which is different from the basis used in [36]. The final solution is nevertheless the same.
This equation has slightly changed with respect to Eq. (36) in [36] because of the change in the basis of potentials. In particular the term proportional to D S 2 changes to compensate the fact that D (2) d is also different so that the result is the same.

N 3 LL heavy quarkonium mass
For the organization of the computation and presentation of the results we closely follow the notation of [7]. In particular we split the total RG improved potential in the following way: where . We then split the total energy into the N 3 LO result and the new contribution associated to the resummation of logarithms. The Swave spectrum at N 3 LO was obtained in Ref. [42] for the ground state, in Refs. [43,44] for S-wave states, and in Refs. [45,46] for general quantum numbers but for the equal mass case. The result for the nonequal mass case was obtained in Ref. [6].
From the RG improved potential one obtains the N i LL shift in the energy levels where the explicit expression for E N i LO (ν) can be found in Ref. [6], and in a different spin basis in Appendix B of [7]. The LO and NLO energy levels are unaffected by the RG improvement, i.e.
We now determine the variations with respect to the NNLO and N 3 LO results. We are here interested in the corrections associated to the resummation of logarithms. In order to obtain the spectrum at NNLL and N 3 LL we need to add the following energy shift to the NNLO and N 3 LO spectrum: which was computed in Ref. [8], and δE nl,RG was computed for l = 0 in Ref. [7], and for l = 0, s = 1 in Refs. [9,10].
To have the complete result for S-wave states, one needs to compute (and add) the new term for l = 0 and δV N i LL r is the delta-related potential contribution to δV N i LL s . The new term generated fromD (2) d then reads where δD (2)N LL d is defined in Eq. (35). The first three lines are generated by the term proportional to δ (3) (r). The last two lines are the contribution to the S-wave energy (l = 0) from the last term of Eq. (21) (the contribution to the P-wave energy, proportional to 1−δ l0 term, is already included in Ref. [7]. Therefore, we do not include it in the expression above). To this contribution we have explicitly subtracted the fixed order contribution already included in the N 3 LO result.
By adding δE new nl,RG N 3 LL to the results computed in these references 5 one obtains the complete result.

Conclusions
In this paper we have computed the α/m 4 , and the α 2 /m 3 , spin-independent potentials (in the Coulomb gauge), and an extra ultrasoft correction that contributes to the S-wave spin-average NNNLL spectrum. We have also obtained the potential RG equation of the delta potential with NLL accuracy (the first nonzero contribution). Combined with the previous results we solve this equation and obtain the complete (potential and ultrasoft) NLL running of the delta potential.
We have also computed the bare and renormalized (soft-)α 3 /m 2 contribution to the spin-independent delta-like potential proportional to [c Combining all these results with the results in [7] and Refs. [9,10] allows us to obtain the S-wave mass with N 3 LL accuracy. The missing terms to obtain the full results are to have the NLL running ofc hl 1 (the associated missing term is of O(T f n f mα 6 ln α) and is expected to be quite small. Its computation will be carried out elsewhere), and a piece of the soft running of the delta potential. This computation will be performed in a separate paper. The magnitude of this contribution is estimated to be smaller compared with the potential running computed in this paper. It is also expected to be smaller than the complete running of the heavy quarkonium potential. Nevertheless, a detailed phenomenological analysis is postponed to future publications. Finally, we remark that significant parts of the computations above are necessary building blocks for a future N 4 LO evaluation of the heavy quarkonium spectrum.
At the order we are working in this paper such differences produce the following difference between the RG equation forD which does not vanish. This difference can be understood through field redefinitions. The field redefinition that moves from the off-shell Coulomb to the on-shell scheme was already discussed in [47,6]. In the second reference the discussion was focused on effects to the spectrum up to O(mα 5 ). We now need to see (the logarithmically enhanced) differences of O(mα 6 ). They can be traced back by using the following Hamiltonian in the Coulomb (Feynman) gauge: where h (0) ∼ mv 2 is the leading order Hamiltonian: and h (2) CG ∼ mv 4 is the relativistic correction, with the explicit potentials: h CG correctly reproduces the O(mα 4 ) spectrum (for the purpose of the comparison we can neglect the O(α) corrections to the static potential: C F α Vs r −C F α 1 r ). This will be enough for our purposes. We now consider the field redefinition that transforms h CG into h OS , the on-shell Hamiltonian: W can be determined from the equation: Since the only possible tensor structure of W is W = W (r 2 )r i the above equation can be written as: V (1) We then obtain and h CG and h ON obviously produce the same spectrum. Therefore, δh cannot produce energy shifts, and any change in the RG equations has to be compensated among different terms. Let us see how it works. h (2) ON produces the differences reported in Eq. (88). Such differences should be eliminated by h |k| 1+2 + 4(p · p )|k| −1+2 + 2(p · k)(p · k) Note that the term proportional to |k| 1+2 in the first line gives a contribution to the potential RG equation through potential loops. It is equivalent to generating a new 1/m 3 potential. The other two terms in the first line do not contribute to the potential RG equation. Looking at the 2nd line, it is also interesting to see that there is a kind of soft contribution, which nevertheless, has ultrasoft in the on-shell scheme. The second term in the second line can also be interpreted as a pure soft contribution. This brings the interesting observation that even if the potential RG equation can be written in a matching scheme independent way, the implicit scheme dependence of the potentials allows for a mixing with the soft computation (at least in the on-shell scheme). Finally, for the dots in the 3rd line we refer to extra contributions toh (4) ON , generated by the field redefinitions, which nevertheless do not contribute to the running.