Relativistic corrections to the static energy in terms of Wilson loops at weak coupling

We consider the ${\mathcal O}(1/m)$ and the spin-independent momentum-dependent ${\mathcal O}(1/m^2)$ quasi-static energies of heavy quarkonium (with unequal masses). They are defined nonperturbatively in terms of Wilson loops. We determine their short-distance behavior through ${\mathcal O}(\alpha^3)$ and ${\mathcal O}(\alpha^2)$, respectively. In particular, we calculate the ultrasoft contributions to the quasi-static energies, which requires the resummation of potential interactions. Our results can be directly compared to lattice simulations. In addition, we also compare the available lattice data with the expectations from effective string models for the long-distance behavior of the quasi-static energies.


Introduction
The near-threshold dynamics of quark-antiquark systems with large masses can be described by a properly generalized nonrelativistic (NR) Schrödinger equation. Such equation can be derived from first principles using effective field theories (EFT's) like potential NRQCD (pN-RQCD) [1,2] (for reviews see [3,4]). This EFT systematically exploits the nonrelativistic scale hierarchy m mv mv 2 . . . , where v is the heavy-quark velocity in the center of mass frame. In the case of unequal heavy quark masses, we will assume m 1 ∼ m 2 ∼ m.
In this paper we will mostly work in the weak-coupling regime, where the EFT can be schematically summarized as i∂ 0 − p 2 2mr − V s (r) φ(r) = 0 + interactions with other low-energy degrees of freedom Here, m r = m 1 m 2 /(m 1 +m 2 ), φ(r) represents the quarkonium wave function, p is the momentum operator conjugate to the distance vector r between the heavy quarks and V s is the color-singlet potential. We denote the color-octet potential by V o . The potential V s consists of the static potential V (0) and its relativistic corrections, which are suppressed by inverse powers of the heavy quark masses. For the general case of unequal masses we have where the spin-independent (SI) 1/m 2 potentials are with L 1 ≡ r × p 1 , L 2 ≡ r × p 2 , and p 1 , p 2 being the momenta of the heavy quarks in the center of mass frame.
In what follows we will not consider spin-dependent potentials.
The potential in Eq. (1) is obtained by matching pN-RQCD to NRQCD [5,6]. There are several ways to perform this matching in practice. In general, the resulting expressions for the relativistic corrections to the potential at different orders in the 1/m expansion, also referred to as (higher order) 'potentials', depend on the matching scheme. Nevertheless physical observables are of course unaffected. For a detailed discussion see Ref. [7].
In this paper we focus on the potentials obtained from Wilson-loop matching. 1 This scheme requires certain gauge invariant NRQCD and pNRQCD Green functions in position space to be equal at a given order in 1/m. On the NRQCD side these Green functions correspond to rectangular Wilson loops with chromo-electric/-magnetic field operator insertions. The potentials are Wilson coefficients of pNRQCD operators. Appropriate correlation functions of these operators are then matched onto the Wilson loops to determine the potentials. For details on the Wilsonloop matching we refer to Refs. [13,14]. In the static limit, i.e. at lowest order in 1/m, the potential reduces to the original Wilson-loop expression for the static potential V (0) [15]. The Wilson-loop expression for the 1/m potential was obtained in Ref. [13], for the 1/m 2 potential without light fermions in Ref. [14], and including light fermions in Ref. [7]. The corresponding expressions in terms of Wilson loops for the SI momentum-dependent 1/m 2 potentials, V p 2 and V L 2 , were first derived in Ref. [16]. The latter reference however missed the 1/m potential 2 , which we also study in this paper.
Only (off-shell) soft modes, i.e. those modes with energies and momenta of order 1/r, contribute to the potential in the Wilson loop scheme. This is achieved by Taylor expanding in powers of 1/m before integrating over the gauge and light-quark degrees of freedom. Thus, soft modes are properly accounted for, but ultrasoft modes still need to be eliminated. At weak-coupling this is achieved by evaluating the Wilson loops strictly in perturbation theory. Hence, the only physical scale left in the Wilsonloop calculation of the potential is 1/r. Still, the definition of the potential in terms of Wilson loops naturally suggests its generalization beyond perturbation theory. We will refer to this generalized potential as the quasi-static energy. 3 To distinguish it from the potential V we will use the symbol E for the quasi-static energy. The definition in terms of Wilson loops has the same form for both quantities, but E is evaluated nonperturbatively, for instance through lattice simulations. Varying the distance r allows us to move continuously from the extreme weak-coupling limit (r → 0) to the strong-coupling regime (r → ∞) of the quasi-static energies. They are therefore ideal observables to study the interplay between both regimes in a controlled manner. 1 In the on-shell matching scheme for the equal mass case they can be found in Ref. [8] with N 3 LO precision in the pN-RQCD power counting, except for the static potential, which is given in Refs. [9,10]. Previous partial results are given in Refs. [11,12]. The expressions in the unequal mass case have been derived in Ref. [7]. 2 as well as the SI momentum-independent 1/m 2 potential Vr, which we do not consider here. 3 In analogy to the potentials, we refer to the different (higher-order relativistic correction) terms in the 1/m expansion of the total quasi-static energy Es, c.f. Eq. (10), as the 'quasi-static energies'. At long distances, the quasi-static energies are suitable for the study of nonperturbative QCD dynamics, e.g. by comparing different models with lattice simulations. In this regime there is a direct connection to the strongcoupling version of pNRQCD, as (omitting possible ultrasoft modes of chiral origin and treating V (0) as small compared to the soft scale) the total quasi-static energy corresponds to the Hamiltonian operator that describes the heavy quarkonium bound state, see Refs. [13,14]. The determination of the quasi-static energy might therefore, for example, enable predictions of the higher excitations of bottomonium and charmonium. For a recent attempt in this direction see Ref. [17].
Even at short distances, the regime we will focus on in this paper, the behavior of the quasi-static energies cannot be obtained from a strictly perturbative computation in the strong coupling α ≡ α s alone. The reason is that the 1/r is generated dynamically by the exchange of an arbitrary number of Coulomb gluons. Including the ultrasoft effects requires a proper resummation of such potential interactions. Typically this is done with the help of EFTs, but it can also be done using a diagrammatic approach, as we will show. Thus, in addition to the perturbative soft part V s , the quasi-static energy E s also includes ultrasoft contributions, which are inaccessible in perturbation theory, even at weak coupling. Therefore it differs from the potential, which enters the Schrödinger equation that describes heavy quarkonium dynamics. Still, both objects are obtained by expanding in powers of 1/m before the integration over the relevant gauge or light-quark dynamical modes.
In Ref. [7], we determined the soft contributions to the quasi-static energies up to O(α 3 /m) and O(α 2 /m 2 ). For the complete results at short distances we also need the ultrasoft contributions. We will compute them in this paper for the 1/m and the SI momentum-dependent 1/m 2 quasi-static energies (E p 2 , E L 2 ). Adding soft and ultrasoft contributions we obtain the 1/m quasi-static energy with O(α 3 ) precision and E p 2 , E L 2 with O(α 2 ) precision. We also discuss the possible resummation of large logarithms for these quantities.
Our results can be directly confronted with the most recent lattice data from Refs. [18,19,20]. 4 One motivation for such a comparison is to perform quantitative checks of ultrasoft effects. The static energy E (0) starts to be sensitive to the ultrasoft scale only at O(α 4 ), i.e. N 3 LO [22,23]. 5 The high α suppression makes it difficult to study ultrasoft physics with this observable. For the quasi-static energies E (1,0) , E (0,1) and E (1,1) ultrasoft effects already contribute at NLO, for E (2,0) and E (0,2) even at leading non-vanishing order. They are therefore a good testing ground for both the ultrasoft physics and the reliability of the weak coupling prediction at a given distance r. This is important for several reasons. For example to give sup-port to perturbative analyses of the heavy quarkonium spectrum that include ultrasoft contributions at weak coupling [25,26,27]. It is also crucial for determinations of the strong coupling constant α from lattice simulations of the static energy [28], as they rely on the correct treatment of ultrasoft effects over the whole range of (short) distances they probe. Finally, we would also like to remark that in N = 4 supersymmetric Yang-Mills theories the static energy E (0) is also sensitive to ultrasoft effects starting at NLO, see [29,30,31,32]. Preliminary results from lattice simulations for this quantity have been reported in [33]. They qualitatively agree with the perturbative predictions (including the ultrasoft contributions) at small distances. A more detailed analysis would however be welcome.

The quasi-static energy and its relativistic corrections
We will use the following definitions for the Wilson-loop operators and their vacuum expectation values. The angular brackets . . . stand for the average value over the QCD action, W 2 for the rectangular static Wilson loop of dimensions r × T W , and where P is the path-ordering operator. Moreover, we define the connected Wilson loop average of operator insertions O 1 , O 2 in the Wilson (static quark) lines at times t 1 and We also define the short-hand notation where T W is the time length of the Wilson loop and T the upper limit of the time integrals in the definitions below. By taking first the T W → ∞ limit, the averages . . . become independent of T W and thus invariant under global time translations, see Ref. [14] for details.
The total color-singlet quasi-static energy is defined as the sum of the kinetic terms and the nonperturbative generalization of the potential in terms of Wilson loops. It reads where d = D − 1 = 3 + 2 ,r = r/r and E n is the chromoelectric field-strength operator to be inserted on the Wilson line associated with the static quark (n = 1) or on the Wilson line associated with the static antiquark (n = 2). The coupling g in the above expressions will be replaced by the bare coupling g B in the actual calculations carried out in the next section. We do not consider spin (S 1,2 )dependent or momentum (p 1,2 )-independent quasi static energies in this paper. The potential V s , which enters the Schrödinger equation that describes heavy quarkonium dynamics, is determined by expanding in powers of 1/m before integrating over any gauge and light-quark dynamical degrees of freedom. For the ultrasoft scales this implies ∆V p 2 /m. Concerning the ultrasoft contributions to the energy of a physical bound state of dynamical heavy quarks such a 1/m expansion would be incorrect, as according to the Virial theorem ∆V ∼ p 2 /m. This is the reason why ultrasoft modes are removed from the potential. They are retained in the quasi-static energy, which can therefore be (formally) regarded as the heavy quarkonium energy in the quasi-static limit: ∆V p 2 /m. The relation between the total (singlet) quasi-static energy E s and the potential in the Wilson loop scheme V s,W is where δE s,US denotes the ultrasoft contribution to the quasi-static energy. The analogous splitting into the corresponding (soft) potential and ultrasoft contribution is understood for each of the quasi-static energies E (...) ... in Eqs. (11)- (16).
In order to compute the potential V s,W (r) in perturbation theory, the ultrasoft contribution must be completely subtracted. The latter can be achieved in dimensional regularization by expanding in the ultrasoft scale ∆V before performing the loop integrations. As ∆V ∝ α this means we simply have to ensure a strictly perturbative calculation in α. Thus, only the soft scale 1/r is left and the loop integrals become homogeneous in that scale. The potentials then take the form of a power series in α ∝ g 2 (and, eventually, in ∆V , when working beyond the order we work at in this paper). In summary we have, We have presented the results for the potential V s,W up to O(α 3 /m) and O(α 2 /m 2 ) in Ref. [7].
In the weak-coupling regime ultrasoft degrees of freedom, with energies and momenta of order ∆V ≡ V Next, we explain how to compute δE s,US in detail.

Ultrasoft calculation
In Ref. [7] we obtained the NLO soft color-singlet contribution to Eqs. (12)- (16). In contrast, the ultrasoft contribution to the NLO quasi-static energies cannot be calculated within conventional perturbation theory. As the ultrasoft scale is only generated dynamically, it requires the resummation of LO potential (Coulomb) interactions. In the Wilson loop calculation the ultrasoft scale ∆V ∼ C A α/r is generated by the ratio of the resummed (exponential) LO expressions for the rectangular color-octet Wilson loop in the numerator and the color-singlet Wilson loop in the denominator in the definition of the quasistatic energies in Eqs. (12)-(16), cf. Eq. (7). This leads to a factor exp(−i∆V t), where t is the (positive) time interval spanned by a (ultrasoft) gluon exchange between the quarks, in the otherwise perturbative calculation. As will be explained below this exponential factor can be effectively treated as a position-space quark-antiquark pair propagator in our diagrammatic approach to the ultrasoft computation, similar to the singlet and octet (boundstate) propagators in pNRQCD [1,2].
To achieve a clean and systematic separation from the soft contribution within dimensional regularization, US . Left-right mirror graphs are understood. Dotted lines represent the A 0 , wiggled lines the A components of a gluon. The time axis in these diagrams extends from left to right. The crosses symbolize the chromo-electric operator insertions at times 0 and t on the upper Wilson line (solid). Vertical lines are associated with potential gluons. The ultrasoft gluon is depicted as diagonal, horizontal or curvy line. We adopt the NRQCD convention for the fermion flow (little arrows on the Wilson lines) explained in Figure 3. Diagrams (1b), (2c), (3e), (3f ), (3g) vanish in Coulomb gauge as they are proportional to scaleless integrals.
we follow the method of regions [34]. That is, we assign a definite energy and momentum scaling to the degrees of freedom involved in the ultrasoft calculation and consistently expand according to the power counting. In our case the relevant degrees of freedom are the ultrasoft gluon with four-momentum q µ ∼ ∆V and potential gluons with energy k 0 ∼ ∆V and three-momentum k ∼ 1/r ∆V . The Wilson lines themselves can be interpreted as static (anti)quarks interacting with both types of gluons. Modes with soft energies ∼ 1/r are insensitive to the ultrasoft scale and therefore decouple from the ultrasoft modes. They are not needed for the ultrasoft contribution to the NLO quasi-static energies. 6 In this calculation only one ultrasoft gluon and a (infinite) number of potential gluons are involved. The couplings of the ultrasoft gluon yield a factor α relative to the LO, whereas any additional potential gluon exchange does not change the power counting.
To compute the correlator E i a (t)E j b (0) c with two chromo-electric field operators (a, b ∈ {1, 2}) we use a (non-standard) diagrammatic approach. The relevant ('connected') diagrams for the ultrasoft part of the quasi-static energies in Eqs. (12)- (16) are displayed in Figures 1 and 2. Vertical lines represent potential gluons, which do not propagate in time (after expansion in the ultrasoft momentum). Potential gluon interactions that will be resummed are not displayed. We will come back to this later. The ultrasoft gluon is depicted as horizontal, diagonal or curvy line. Before we comment on the selection criteria for these diagrams, we describe how to construct and evaluate them.
The procedure is as follows. We first locate the vertices of the ultrasoft and potential gluons at specific times. As a consequence there is a priori no energy conservation at these vertices. The two vertices for the chromo-electric field operator insertions on the static quark (Wilson) lines are fixed at times 0 and t, respectively. The time coordinates of the other vertices have to be integrated over. The energies in the potential gluon propagators are expanded Fig. 3. Feynman rules for an E i (t) operator insertion (denoted by a cross in the diagram) on a static quark (Wilson) line. Dotted and wavy lines represent A0 and A gluons, respectively. All momenta (k, l, q, p) are incoming. Note that, in contrast to the usual convention for Wilson lines, we treat the antiquark (lower Wilson line) as a particle living in the antirepresentation of SU (3), i.e. the fermion flow (little arrows) of the antiquark is the same as for the quark. The corresponding Feynman rules for the antiquark are then obtained by replacing out and do not contribute at NLO. 7 The integration over these energies can therefore be carried out immediately yielding delta functions of time intervals. Many of the time integrations can now be performed trivially. The result is that vertices connected by potential gluons are all located at the same time. This has been indicated already in our graphs in Figures 1 and 2, where time proceeds from left to right.
The vertices for the chromo-electric operator insertions are the same as the ones given in Ref. [7] and shown in Figure 3. For the Wilson line interaction of an A 0 gluon, the quark and gluon propagators, and the gluon self-interactions we can use standard NRQCD Feynman rules at leading order in 1/m. In particular the static quark propagator in position space is nothing else than a θ-function of the time difference. We emphasize however that, according to the method of regions, we have to expand the potential gluon propagators in the energy and ultrasoft momenta.
The way we implement the resummation of potential gluon interactions, which actually generates the ultrasoft scale, in our diagrammatic calculation is probably best explained at a concrete example. Let us consider diagram (3d) in Figure 1. In Figure 4 we have depicted it with time labels for the vertices on the upper Wilson line and indicated a possible choice for the potential three-momentum flow. The momentum k is the Fourier conjugate to the distance r. After expansion and integration over the energies associated with the potential propagators the two vertices on the lower Wilson line are fixed at times −t and 0, respectively. We have highlighted the time interval, where the static quark-antiquark system is in a color octet state, because the ultrasoft gluon has been radiated off, by thick red static quark (Wilson) lines. During this time interval it is crucial to take into account the exchange of an infinite number of potential Coulomb gluons, between the two Wilson lines. The resummation of such Coulomb interactions yields the LO expression for the average of a rectangular Wilson loop with extensions t 2 − t 1 r in the color octet configuration, see e.g. Ref. [2]: The analogous expression with the octet static potential V (0) o (r) replaced by the singlet static potential V (0) (r) holds for the color singlet configuration. Outside the red marked time interval in Figure 4 this exponential however exactly cancels with the singlet Wilson loop average in the denominator of E i 1 (t)E j 1 (0) , whereas inside the octet time interval we can assign a resummed 'propagator', θ(t 2 −t 1 ) exp[−i∆V (t 2 −t 1 )], to the static quark pair.
In Feynman gauge we thus have where D = d + 1 = 4 + 2 , q denotes the ultrasoft fourmomentum, g B is the bare version of the coupling, and we have multiplied by 2 to add the left-right mirror graph exploiting time reversal symmetry. We only kept the leading order in the ultrasoft momentum expansion. Higher orders would come with higher powers of α and are beyond the order we are interested in. For the same reason, we only need the LO expression in D dimensions and we can replace the bare coupling g B by the renormalized g as A few practical comments on this example are in order: After performing the last two time integrations we have written Eq. (22) in a suggestive form in terms of momentum space building blocks. The analogous form can be directly obtained for all diagrams in Figures 1 and 2 by assigning the expression in square brackets to the momentum space quark-pair 'propagator' between the two crossed vertices and i/(−q 0 − ∆V ) to the other quark-pair 'propagators' between an uncrossed and a crossed vertex, and using momentum space Feynman rules for the vertices and gluon propagators otherwise.
In Coulomb gauge all diagrams with an A 0 ultrasoft gluon vanish, because the q integral is scaleless after the ultrasoft momentum expansion. To change from Feynman to Coulomb gauge in the diagrams with an A ultrasoft gluon we only have to replace δ ij → δ ij − q i q j /q 2 in the numerator of the A gluon propagator. It is easy to see that, due to the ultrasoft momentum expansion, this will result in nothing more than an additional overall factor (D −2)/(D −1) for the Coulomb gauge diagrams. In other words, due to gauge invariance, the sum of all diagrams with an ultrasoft A 0 gluon must be equal to 1/(D − 2) times the sum of all diagrams with an ultrasoft A gluon. Note that in order to check this explicitly, it is crucial to use the D-dimensional expression for ∆V in Eq. (24).
Also, from our calculation in Eqs. (22) and (23) and the Feynman rules in Figure 3 it should be clear that the diagrams (3d) in Figures 1 and 2 yield exactly the same result and in fact all diagrams in Figure 1 equal their counterpart in Figure 2. As we will see below this is related to Poincaré invariance. We therefore conveniently define the quantity with n = 1, 2 for the 1/m and 1/m 2 quasi-static energies, respectively. We can then write the ultrasoft contributions to the SI quasi-static energies in Eqs. (12)- (16) as δE (1,1) From these equations we can directly read off This verifies Poincaré invariance, see Ref. [35]. It requires that, given that there is no ultrasoft contribution at O(α 2 /m 0 ), the O(α 2 /m 2 ) ultrasoft part of the total energy of the quark-antiquark system must be proportional to 1/m 2 r = 1/m 2 1 + 1/m 2 2 + 2/(m 1 m 2 ). The corresponding relations for the soft potentials at the same order were checked in Ref. [7].
From the sum of all relevant ultrasoft one-loop diagrams we obtain the total bare result in both Coulomb and Feynman gauge with ∆V as given in Eq. (24). We emphasize that the connected diagrams in Figure 1 or Figure 2 are not the only diagrams contributing to Eq. (33). In our definition a 'connected' graph cannot be split into two (lower loop) Wilson loop diagrams by cutting the quark and the antiquark line at the same time, i.e. by a vertical cut through the diagram. We also have to take into account disconnected (2-Wilson-line reducible) one-loop diagrams as well as (products of) tree-level diagrams from the E i a (t) E j b (0) term in the definition of the connected Wilson loop average, Eq. (8). Their net effect is however only to remove the terms proportional to C 2 F and C 3 F from the connected diagram contributions as explained in detail in the next section. Keeping this in mind it is straightforward to check that the diagrams given in Figures 1 and 2 give the complete C F contribution we are interested in: Every additional diagram with one ultrasoft gluon exchange you can draw, is either disconnected or redundant, in the sense that it is one of the diagrams with additional Coulomb interactions, which we automatically take into account by our resummation. The individual contributions (∝ C F ) to W ij n from each diagram in Figures 1 and 2 are given in Appendix A in Coulomb and Feynman gauge respectively.
We can now use Eq. (33) together with the LO expression Eq. (24) for ∆V in the ultrasoft contributions Eqs. (26)- (29) and add the bare results for the corresponding potentials in Ref. [7] to obtain finite expressions for the quasi-static energies. With the bare coupling β 0 4π 1 (34) expressed in terms of the MS renormalized coupling α ≡ α(ν) ≡ g 2 ν 2 4π , we find after setting → 0.
2.2 Cancellation of C 2 F and C 3 F terms As already mentioned in the previous section, the Feynman diagrams we depict in Figure 1 and Figure 2 are the connected (2-Wilson-line irreducible) diagrams relevant for our ultrasoft calculation. It is clear that most of these diagrams are accompanied by a number of disconnected diagrams, which differ from the connected diagrams only by the order of the vertices on the quark lines. For each connected diagram in Figure 1 and Figure 2, let us define a diagram 'family', which also includes its disconnected relatives.
Regarding the associated color factor, within each family the connected diagram contains a maximally non-Abelian term, i.e. a term that is linear in the fundamental Casimir constant C F . The color factor of the disconnected diagrams consists of terms with at least two powers of C F . Summing all diagrams of a family and subtracting the corresponding E i a (t) E j b (0) term in the definition of the quasi-static energies according to Eq. (8) and Eq. (25), all terms ∝ C n F with n ≥ 2 cancel out and only the maximally non-Abelian term survives. This mechanism is very much reminiscent of the well-known non-Abelian exponentiation [36,37,38] of the static potential. For the case of the quasi-static potentials/energies we are lacking however a general proof. We have therefore checked the cancellation of the C 2 F and C 3 F terms present in the diagrams with one ultrasoft loop explicitly for each family. In the following we will explain how this works for two diagram families, whose connected graphs are displayed in Figure 1. From this it should be clear how to show the cancellation also for the other families. In particular, all statements in this section can be straightforwardly extended and are therefore equally valid also for the diagrams in Figure 2, which contribute to δE (1,1) U S . Let us first introduce the notion of a color-stripped diagram. We write e.g.
The graph on the RHS symbolizes the original Feynman diagram with removed (fundamental, anti-fundamental or adjoint) color generators. The red zigzag line represents the ultrasoft gluon.
On the level of such color-stripped diagrams we have This can be seen easily as follows. The only difference between the two graphs on the LHS is the lower quark propagator. Labeling the vertex of the ultrasoft gluon with the lower quark line by the time t this propagator is θ(t < t) and θ(t > t) in the two diagrams, respectively. After Coulomb resummation some of the θ-functions actually come from the octet bound state propagators and all diagrams with an ultrasoft gluon have a factor exp(−i∆V |t − t |) in common, where in this case t = 0. For the purpose of this section we can however effectively assign a θ-function to each static quark line in the color-stripped graphs. In the sum of the two one-loop diagrams the θfunctions drop out. As we assume t > 0 by definition, the upper 'quark propagator' θ(t > 0) can also be dropped. Furthermore, by shifting the potential three-momentum k we have (with q being the ultrasoft four-momentum, in the integrand of the two one-loop diagrams. Thus we are left with the product of the two tree-level diagrams on the RHS, which will eventually cancel with the corresponding term from E i a (t) E j b (0) ). Note that in this case the LO in the ultrasoft momentum expansion (|q| |k| ∼ 1/r) vanishes due to the antisymmetric qintegrand. The NLO in this expansion is nonzero and represents the contribution at the order we are interested in.
With these preliminaries it is straightforward to show the cancellation of the C 2 F and C 3 F terms for most of the other diagram families, and we can tackle already the most complicated case, namely the family 3e. Without the maximally non-Abelian C 2 A C F term the sum of the diagrams in this family can be written as Exploiting again the static quark line θ-functions, we have used the equalities and their analogs for the graphs, where the ultrasoft gluon is up-down mirrored, in Eq. (44) to write the C 2 F term as a sum of products of one-loop and tree-level diagrams with a single insertion of an E-field operator. Recall that in our notation the operator insertion vertex on the left is always located at time 0 and the one on the right at time t with t > 0. The two expressions on the RHS of Eq. (45) and Eq. (46) are therefore unequal. The C 2 F term of Eq. (44) will be directly cancelled by the corresponding C 2 F term from the E i 1 (t) E j 1 (0) part of the definition of the quasi-static energies.
The cancellation of the C 3 F term is a bit more subtle. First we note that by combining the θ-functions from the upper and lower static quark lines. Also, (after shifting a potential threemomentum) we directly have With the equalities in Eqs. (45)- (48) and their analogs for the mirror graphs we can write the C 3 F term of Eq. (44) in factorized form: Apart from the overall factor 1/2 this is what one might naively expect for the C 3 F term from E i 1 (t) E j 1 (0) at NLO in the ultrasoft loop expansion. On the first sight there therefore seems to be miscancellation of the C 3 F terms. Let us however take a closer look at Eq. (49). Similar to Eq. (47) the expressions in round brackets can be factorized further and we are left with Now, recall that in the definition of E 1 (t)E 1 (0) there is the singlet Wilson loop average in the denominator. At LO we simply have W 2 = exp(−iV s T W ), but at the NLO in the ultrasoft loop expansion we obtain as an additional contribution to the quasi-static energies. Thus the total C 3 F term from the E i 1 (t)E j 1 (0) vanishes.
It is easy to see that the same happens separately for the C 3 F term associated with E i 1 (t) E j 1 (0) . There we have two times Eq. (50) plus two times the RHS of Eq. (51) from the expansion of the two denominators in E i 1 (t) E j 1 (0) . This completes the proof of the cancellation of the C 2 F and C 3 F terms related to family 3e. Analogous proofs hold for families 3f and 3g. Adding the contributions to W ij n from diagram families 3e, 3f and 3g the LO terms in the ultrasoft momentum expansion, i.e. terms ∼ α 3 ∆V −1−n r −4 , cancel out, but a contribution ∼ α 3 ∆V 1−n r −2 ∼ α 4−n r n−3 , which is the order of our interest, survives. We note that the diagram families (2c, 3e, 3f , 3g) discussed above actually vanish in Coulomb gauge.
For the families with an ultrasoft A gluon (3d, 4a, 4b, 5a) relevant also in Coulomb gauge the proofs are, in fact, even simpler than the ones presented here and straightforward to do with the above techniques. We have thus shown 'non-Abelian exponentiation' of the 1/m and 1/m 2 quasi-static energies at one ultrasoft loop.

Resummation of large logs
Our finite order results for the quasi-static energies in Eqs. (35)-(39) depend on the factorization/renormalization scale ν. Setting ν ∼ 1/r large logarithms associated with the soft scale are resummed in the running coupling α(ν). For the complete resummation of large logarithms, also ultrasoft logarithms have to be resummed. In this section we discuss how this can be achieved.
The bare potential (in the Wilson loop scheme) can be written in terms of the renormalized potential and its counterterm as The counterterm δV s for the potentials in terms of Wilson loops can be determined at leading order from the (IR) 1/ terms of the soft contribution computed in Ref. [7]. It can also be obtained (at leading order) from our ultrasoft calculation in the present paper. Indeed the 1/ terms of both computations match, which guarantees that we obtain finite results in Eqs. (35)- (39). At leading order the structure of the counterterm is (cf. Eq. (55) below) where n and F depend on the specific potential V . This information by itself is not enough to obtain the leadinglog (LL) running for each potential in the Wilson loop scheme. Still, this problem can be bypassed at LL, since EFT arguments tell us that the structure of the LO overall counterterm should be (note that ν here refers to ν us ) where G is a function of the singlet/octet static potentials, and V A = 1 at this order.
Let us now link this discussion to previous determinations of the counterterm δV s in the direct context of EFTs. If the counterterm is determined in an EFT calculation in terms of the Wilson coefficients of the EFT it is possible to resum the ultrasoft logarithms into the potentials by solving the associated renormalization group equations (RGEs). Since our derivation of the quasi-static energies is based on the method of regions rather than on an EFT, we have a priori no RGEs at hand to resum such logarithms. Using the EFT, the counterterm and the associated RGE for the next-to-leading-log (NLL) ultrasoft running of V s were obtained in Ref. [39] 8 (the LL ultrasoft running of V s was obtained in Ref. [46]). These expressions are suitable for the potentials in the on-shell or Feynman/Coulomb gauge off-shell matching schemes. However, as explained in Ref. [7], the counterterm of Ref. [39] fails to cancel the divergences of the individual potentials in the Wilson loop matching scheme. Therefore we proposed a modified counterterm, that removes the divergences from the individual potentials at LO for that scheme [7]. Such a counterterm directly follows from the computation and the knowledge of the structure of the counterterm in Eq. (54), as we discussed above. We conjecture that the structure of this counterterm is preserved at NLO. Adding the NLO piece according to Ref. [39] gives the following expression for the counterterm in the Wilson loop matching scheme (W ): At LO this is correct, because δV s,W correctly reproduces the UV divergence of the ultrasoft contribution to all quasistatic energies we have calculated. At NLO and beyond it is only a conjecture, because in our approach we have no proper control on (potential) operator mixing effects in the RGE. We will nevertheless use the resulting expressions as an educated guess for the renormalized NLL potentials V s,W . We hope that an EFT analysis will eventually settle this issue.
The RGE then reads where Here and in the following the Wilson loop scheme for the potentials is understood and we drop the corresponding subscript W for brevity. We also set the pNRQCD Wilson coefficient V A = 1, which is consistent at LL [43] and NLL [44].
We write the running potential as with δV s,RG (r; ν, ν us = ν) = 0. We now only consider the 1/m and the 1/m 2 SI momentum-dependent potentials. Solving the RGEs and truncating at the appropriate order in α we obtain δV (1,0) RG (r; ν, ν us ) = 2 δV (2,0) δV (1,1) δV (1,1) where We computed the MS renormalized (Wilson loop) potentials at ν us = ν in Ref. [7]. For convenience we quote the results here: 9 Including the resummation of ultrasoft logarithms the total quasi-static energy then reads E s,RG (r) = V s,RG (r; ν us ) + δE US (r; ν us ) , and similarly for the individual quasi-static energies we discuss in this paper. The corresponding MS renormalized ultrasoft parts are Formally, these ultrasoft contributions are independent of ν, as in pNRQCD the soft scale has already been integrated out. In practice they do depend on ν once we allow for the running of α(ν). This however affects the result only beyond the order to which the ultrasoft computation has been carried out. For a meaningful evaluation of the (RG improved) quasi-static energies ν ≈ 1/r and ν us ≈ C A α/r must be chosen.

Comparison with lattice data
Monte Carlo lattice simulations of the 1/m and 1/m 2 SI momentum-dependent quasi-static energies have been performed over the years [21,19,20]. The data for the SI momentum-dependent 1/m 2 contributions is usually displayed in terms of the following linear combinations of the quasi-static energies: 10 Therefore, we will directly confront our results with these quantities in the following section. On the lattice side we will use the most recent results obtained in Refs. [19,20], also for the comparison with the 1/m quasi-static energy.

Short distances
Lattice data is very scarce at short distances. It is also worth mentioning that at short distances one can easily see finite discretization effects in the lattice data. On top of that, due to possible power-like divergences in the lattice simulation, we cannot compare with the lattice data in absolute values. We are thus forced to consider differences between quasi-static energies computed at different distances. We then choose to normalize our results to the lattice point at the shortest available distance for each observable. This leaves us with at most three independent lattice points (in fact only two for E b and E c ) for distances smaller than 0.5 r 0 ≈ 0.8 GeV, which is already at the borderline for our weak coupling results to be reliable. Therefore, the whole discussion in this section will stay at the qualitative level.
The analytic results we present in the plots of this section are based on Eqs. (35)- (39). For the resummation of the ultrasoft logarithms we will use the the results of Section 3. We will first define the different curves shown in these plots and discuss them in comparison with the available lattice data afterwards. The different colors of the lattice data points in our plots (with hardly visible error bars) indicate different β values (blue: β = 5.85, green: β = 6, brown: β = 6.2), see Refs. [19,20]. For E (1,0) , E b and E c Eqs. (35)- (39) give the first two (LO and NLO) nontrivial terms, whereas for E d and E e they only yield the first (NLO) nontrivial term in the weak coupling expansion. The scale dependence of the coupling constant is uncancelled beyond NLO. We use this fact to grasp the dependence of our predictions on higher order corrections. On the one hand we produce curves with a fixed, r-independent factorization scale ν. For these, we take as a reference scale the shortest distance of the available lattice points, and vary ν around this value. For E (1,0) , E d and E e , we take ν = x/0.153990 r −1 0 and for E b and E c we take ν = x/0.371627 r −1 0 . Here and in the following the parameter x always represents a number of order 1. On the other hand we also produce curves with ν = x/r. The obvious motivation of this is to resum logarithms of νr. 11 This is however not sufficient to resum all large logarithms, as there are also ultrasoft logarithms. In order to perform a complete resummation of logarithms, we have to use Eqs. (58)-(75) and count α(ν) ∼ α(ν us ), so that the resummation kernel F (ν, ν us ) ∼ 1. With this counting the O(α 2 ) term in the expression for the 1/m and the O(α) term in the expression for the 1/m 2 quasi-static energies, represents the results with LL accuracy. Our complete expressions for the quasi-static energies according to Eq. (70) represents an educated guess of the NLL expression and we use it as an estimate of the next order. The dependence on ν and ν us cancels in the theoretical expression at NLL but not beyond. Again, this residual dependence may give hints on the size of higher order terms. We will generate lines with ν = x/r and ν us = x us C A α s (ν)/r, with x ∼ x us ∼ 1. Note that the RG improved expressions 11 Setting ν = 1/r may introduce problems with renormalons (see the discussion in Ref. [47]), which we do not explore here. The first reason is that a complete analysis of the renormalon structure of these potentials is missing (besides the trivial fact that E We first consider the 1/m quasi-static energy at LO and NLO. We observe that the uncertainties are quite large. This is reflected in a strong scale dependence, which could be expected as we are evaluating our result at very low scales and already the LO is O(α 2 ). Even using the one-loop or two-loop expression for the running α produces sizable differences. Nevertheless, it is comforting that we can find reasonable values of ν such that the curves match the lattice data quite well (ν ≈ 3.4/0.264 r −1 0 ≈ 5.2 GeV for the NLO and ν ≈ 1.6/0.264 r −1 0 ≈ 2.4 GeV for the LO), even though the difference between LO and NLO is relatively large. To illustrate this difference we show curves for the LO and NLO results with fixed ν = 2 × 1/0.263821 r −1 0 ≈ 3 GeV in Figure 5. We also show the corresponding curves with ν = 3/r. They hardly differ from the choice ν = 2/r and lie below the lattice points, but are still consistent with them within the expected uncertainties. The effect of setting ν ≥ 2/r is making the series more convergent: The difference between LO and NLO becomes smaller compared to the fixed scale results as observed in Figure 5. However, ν ≤ 1/r is a too small scale and the perturbation series blows up (which may be due to renormalons). A systematic resummation of all large logarithms requires setting ν us ∼ C A α(ν)/r and ν ∼ 1/r. The additional resummation of ultrasoft logarithms further improves the convergence (LO and NLO curves are closer), but the curves also lie further below the lattice points, see Figure 5. Just like for the 'fixed order predictions' with r-dependent factorization scale, we have to choose relatively large values of ν to avoid unphysical behavior.
Next we consider E b . The situation here is far from optimal, as illustrated in Figure 6. We have one less lat- tice data point (the one at shortest distance) than for E (1,0) , which makes the region, where we can compare with the lattice data, even smaller. There are also visible finite size effects in the lattice data, when comparing the points for different values of β. This indicates that a global negative shift of our lines would be appropriate in the comparison with the lattice data. We do not implement this shift explicitly in Figure 6, but it is understood in the following discussion. With fixed ν the agreement is reasonable. For ν ≈ 1.5 × 1/0.371627 r −1 0 the NLO result agrees well with the data and so does the LO result, as there are only small differences between the two curves. Choosing ν ∼ 1/r deteriorates the agreement. As for the 1/m quasi-static energy, setting ν ≤ 1/r causes the perturbation series to blow up (maybe due to renormalons) as the scale becomes too small. The complete resummation of logarithms yields strongly scale-dependent results. If one takes ν us = C A α(ν)/r and ν = 4/r the agreement with the data is very good at NLL but the difference to the LL prediction is very large. Also, we observe a strong dependence of the result on x us (not displayed) in the range between one and two. In any case the uncertainties are large.
For E c , the situation is similar as for E b as far as the number and location of the lattice points are concerned, see Figure 7. Note that there are again sizable finite size effects in the lattice data. So, it would be reasonable to add a global shift upwards to all our lines when making the comparison with data. Working at fixed order with fixed scale the slope is smaller than predicted by the data, as can be seen in Figure 7  the situation, unless we choose x ≈ x us ≈ 1, where one can get good agreement with the data. Nevertheless, once more the scale dependence is large.
We now turn to E d and E e . They are linear combinations of E (2,0) L 2 (r) and E (2,0) p 2 (r). This makes them particularly interesting, as they vanish at O(α). Previous lattice determinations obtained zero slopes within uncertainties [21]. Nevertheless, in a dedicated series of works [19,20], the SI momentum-dependent quasi-static energies have been computed with increased precision on the lattice and a non-trivial behavior has been found at short distances. A non-zero slope is also predicted by our computation of the quasi-static energies. Therefore, it is very interesting to check whether they follow the same pattern, as we will do in the following.
For E d the agreement, if we choose a r-independent scale ν, is bad. The hardly visible slope of the corresponding curve in Figure 8 is by far not enough to match the lattice data. The asymptotic behavior is however correct (in sign) assuming that the ultrasoft log dominates, which is not the case for the value of ν used in Figure 8 though. If we set ν = 2/r the agreement is much better. (It gets even better for x ≈ 1). The resummation of ultrasoft logs for reasonable values of x and x us (optimally x ≈ 2 while 2 ≥ x us ≥ 1) further improves the agreement. However, the observed scale dependence varying x between one and two is relatively large.
The situation for E e is similar. Setting ν constant the agreement is bad, as can be seen in Figure 9. If we set ν = 1/r the agreement is much better (optimal for x ≈ 0.83). One can also find good agreement if we resum the ultrasoft logs for reasonable values of x and x us . Again, the uncertainties due to scale variations in the r range where lattice data is available are sizable and prevent a more quantitative analysis.
Overall, we find that for all studied observables our predictions can be made compatible with the available lattice data. With this statement we mean that for (more or less) reasonable values of the factorization scale we can get reasonable agreement with the lattice data, and/or that the expected size of higher order corrections is large enough to accommodate the difference to the lattice data. We cannot make quantitative statements though. Whereas adjusting the factorization scale or incorporating the resummation of logarithms can improve the agreement with the lattice data in some cases, we do not observe a general pattern for all potentials. Estimates of higher order effects from scale variations produce large effects, related to the fact that we have evaluated our predictions at quite low scales (large distances). Also, as already mentioned, the lattice data is rather scarce. Still, we find it quite remarkable that the asymptotic behavior (the sign of the slope) predicted by the ultrasoft logarithms ∝ ln(C A α) in E d and E e complies with the lattice data.

Long distances
It is also interesting to consider the long-distance behavior of the quasi-static energies. Their behavior cannot be predicted by QCD using analytic tools. Still, it can be simulated by computations in the discretized version of QCD, i.e. on the lattice. One can then compare these numerical results with model predictions. Here we focus on the effective string theory predictions for the 1/m and the SI momentum-dependent quasi-static energies obtained in Ref. [48]. For the SI momentum-dependent quasi-static energies, they are actually equivalent to the predictions obtained first with the minimal area law [16]. 12 The latter however misses the 1/m potential. The string model gives the following results at leading order in the 1/(rΛ QCD ) expansion: These expressions only depend on the string tension σ. 13 The string tension can be obtained from fits to the static potential at long distances. The number we use is taken from [50] and reads: σ = 1.3882r 2 0 . Therefore, the above expressions are in fact predictions. On the lattice side, to obtain absolute normalizations is problematic. Hence, just like at short distances we compare energy differences. We choose to normalize our prediction to the longest distance (rightmost) lattice point available, where the above models are expected to work best. The comparison with the lattice is shown in Figures 10, 11 and 12. There is qualitative agreement between the lattice and the effective string theory predictions. 14 On the other hand it is evident that the agreement is far from perfect and the slopes demand for corrections to the leading string behavior. Such corrections can in principle be computed in the same way as for the static potential, where they give rise to the Lüscher term [51]. It would be interesting to see if the analogous computations for the quasi-static energies (in case they happen to be free of new constants) improve the agreement with the lattice data. 12 Other models also share the same behavior in the longdistance limit, see the discussion in Ref. [49]. 13 The dependence on µ1 that appears in E (1,0) vanishes when comparing with lattice simulations as we will only compare differences of static energies evaluated at different r. 14 In Figure 11 it is obvious that there are still sizable finite size effects as there is no good scaling between the lattice data points for β = 5.85 (blue) and β = 6 (green).

Conclusions
In this paper we have determined the short-distance behavior of the nonperturbative generalizations of the 1/m potential and of the 1/m 2 SI momentum-dependent heavy quarkonium potentials, in terms of Wilson loops (as defined in Section 2.1), with O(α 3 ) and O(α 2 ) accuracy, respectively. We refer to these objects as quasi-static energies. Our results can be found in Eqs. (35)- (39). For our (NLO) computation it is crucial that the quasi-static energies start to be sensitive to the ultrasoft scale ∼ C A α/r. Therefore, their calculation requires a resummation of an infinite number of potential interactions to account for the ultrasoft effects.
We have computed these ultrasoft contributions to the quasi-static energies in a diagrammatic approach. This included Feynman graphs with one ultrasoft and up to three potential loops. We expect that a proper EFT setup will significantly simplify such computations. At the same time this would allow for a systematic RG analysis beyond LL precision. With our method we were able perform the correct LL resummation of ultrasoft logarithms, but can only provide an educated guess for the NLL results.
Our results should reproduce the short-distance limit of lattice simulations of the same quasi-static energies. Some of these were evaluated long ago in Ref. [21], and more recently in Refs. [18,19,20]. We have performed a preliminary comparison with the lattice data [18,19,20] at short distances in Section 4.1. As there is not enough lattice data at short distances, we were unable to perform a quantitative analysis. We had to compare at rather low scales. This causes large uncertainties from scale variations. Within these uncertainties we can find qualitative agreement with the data. It is quite remarkable that the sign of the slopes obtained by the lattice simulations agrees with the asymptotic behavior predicted by our weak coupling analysis. Whereas for three of the considered quasi-static energies, the behavior is dominated by the soft LO contribution, for two of them the sign of the slope is dictated by the ultrasoft contribution. This makes them particularly interesting. Those quasi-static energies vanish at O(α), but are non-zero at O(α 2 ). Earlier lattice simulations predicted an approximately zero slope for them in the short distance limit. The most recent simulations show a non-zero behavior consistent with our results, albeit within large uncertainties.
We also compare the lattice data with the expectations from effective string models at long distances. We find qualitative agreement, but there is obvious room for improvement. In the long term one could try to study (qualitatively) how the short and long distance predictions (for a specific model) combine.
The fact that the sum of the W X,ij n,CG gives again the gauge invariant result in Eq. (33) is a strong cross check of our computation.