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

We consider the O(1/m)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {O}}(1/m)$$\end{document} and the spin-independent momentum-dependent O(1/m2)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {O}}(1/m^2)$$\end{document} 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 O(α3)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal{O}(\alpha ^3)$$\end{document} and O(α2)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal{O}(\alpha ^2)$$\end{document}, 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 (EFTs) like potential NRQCD (pNRQCD) [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 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 pNRQCD to NRQCD [5][6][7]. 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. [8].
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 Wilson-loop matching we refer to Refs. [14][15][16]. 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) [17]. The Wilson-loop expression for the 1/m potential was obtained in Ref. [14], for the 1/m 2 potential without light fermions in Refs. [15,16], and including light fermions in Ref. [8]. The corresponding expressions in terms of Wilson loops for the SI momentumdependent 1/m 2 potentials, V p 2 and V L 2 , were first derived in Ref. [18]. 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 lightquark 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 Wilson-loop 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- 1 In the on-shell matching scheme for the equal mass case they can be found in Ref. [9] with N 3 LO precision in the pNRQCD power counting, except for the static potential, which is given in Refs. [10,11]. Previous partial results are given in Refs. [12,13]. The expressions in the unequal mass case have been derived in Ref. [8]. 2 As well as the SI momentum-independent 1/m 2 potential V r , 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 quasistatic energy E s , cf. Eq. (10), as the 'quasi-static energies'. static energies. They are therefore ideal observables to study the interplay between both regimes in a controlled manner.
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 strong-coupling 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. [14][15][16]. 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. [19].
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 so-called 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 quasistatic 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. [8], 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. [20][21][22]. 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 [24,25]. 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 nonvanishing 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 support to perturbative analyses of the heavy quarkonium spectrum that include ultrasoft contributions at weak coupling [27][28][29]. It is also crucial for determinations of the strong-coupling constant α from lattice simulations of the static energy [30], 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 [31][32][33][34]. Preliminary results from lattice simulations for this quantity have been reported in [35]. 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 for the rectangular static Wilson loop of dimensions r × T W , 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 Refs. [15,16] 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 P R is the center of mass momentum operator, p = −i∇ r , M = m 1 + m 2 , the reduced mass m r = m 1 m 2 /M, and is the static energy [17]. The formal definitions for the higher-order quasi-static energies to O(1/m) can be found in Ref. [14] and to O(1/m 2 ) in Refs. [8,15,16]. Here we recall the ones relevant for the present work. At O(1/m) we have The SI momentum-dependent 1/m 2 potentials are parametrized in complete analogy to Eqs. (3)-(5) and given by 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 quasistatic 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. [8].
In the weak-coupling regime ultrasoft degrees of freedom, with energies and momenta of order ΔV ≡ V (0) o − V (0) ∼ C A α/r , contribute to the quasi-static energies through δ E s,US . These contributions to the Wilson loops correspond to the contributions from the lower (ultrasoft) energy regions in an asymptotic expansion of the loop integrations according to the method of regions [36]. This implies that (some of) the loop integrals are sensitive to the ultrasoft scale. In analogy to Eq. (18) we write Next, we explain how to compute δ E s,US in detail.

Ultrasoft calculation
In Ref. [8] 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 quasi-static 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 (bound-state) propagators in pNRQCD [1,2].
To achieve a clean and systematic separation from the soft contribution within dimensional regularization, we follow the method of regions [36]. 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 fourmomentum 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 Figs. 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 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 Figs. 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. [8] and shown in Fig. 3. For the Wilson-line interaction of an A 0 gluon, the 7 Strictly speaking, there are higher-order terms in the potential energy expansion contributing in the Feynman gauge diagrams (3e), (3 f ) and (3g), but they cancel in the sum. 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  Fig. 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 Fig. 4 this exponential, however, exactly cancels with the singlet Wilson-loop average in the denom- 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 g B = g + O(g 3 ).
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 Figs. 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 δ i j → δ i j − 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 Fig. 3 it should be clear that the diagrams (3d) in Figs. 1 and 2 yield exactly the same result and in fact all diagrams in Fig. 1 equal their counterpart in Fig. 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 From these equations we can directly read off δ E (1,1) This verifies Poincaré invariance; see Ref. [37]. 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 quarkantiquark 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. [8].
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 Figs. 1 or 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 Figs. 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 i j n from each diagram in Figs. 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. [8] to obtain finite expressions for the quasi-static energies. With the bare coupling expressed in terms of the MS renormalized coupling α ≡ α(ν) ≡ g 2 ν 2 4π , we find 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 term in the definition of the quasi-static energies according to Eqs. (8) and (25), all terms ∝ C n F with n ≥ 2 cancel and only the maximally non-Abelian term survives. This mechanism is very much reminiscent of the well-known non-Abelian exponentiation [38][39][40] 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 Fig. 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 Fig. 2, which contribute to δ E (1,1) U S . Let us first introduce the notion of a color-stripped diagram. We write e.g.
. Figure 1 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- (42) 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 q-integrand. 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 Eqs. (45) and (46) are therefore unequal. The C 2 F term of Eq. (44) will be directly canceled 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 = exp(−i V 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 3 f and 3g. Adding the contributions to W i j n from diagram families 3e, 3 f and 3g the LO terms in the ultrasoft momentum expansion, i.e. terms ∼ α 3 ΔV −1−n r −4 , cancel, 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, 3 f , 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. [8]. 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 leading-log (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. [41] 8 (the LL ultrasoft running of V s was obtained in Ref. [48]). These expressions are suitable for the potentials in the on-shell or Feynman/Coulomb gauge off-shell matching schemes. However, as explained in Ref. [8], the counterterm of Ref. [41] fails to cancel the divergences of the individual potentials in the Wilson-loop matching scheme. Therefore we proposed a modified counterterm, which removes the divergences from the individual potentials at LO for that scheme [8]. 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. [41] 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 [45] and NLL [46]. 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 (2,0) δV (2,0) δV (1,1) δV (1,1) where We computed the MS renormalized (Wilson-loop) potentials at ν us = ν in Ref. [8]. For convenience we quote the results here: 9 Including the resummation of ultrasoft logarithms the total quasi-static energy then reads and similarly for the individual quasi-static energies we discuss in this paper. The corresponding MS renormalized ultrasoft parts are δ E (2,0) δ E (2,0) δ E (1,1) 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) quasistatic 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][22][23]. 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. [21,22], 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 10 In the literature often the notation V b , V c , etc. is used. We stick, however, to the notation E b , E c , etc. to distinguish the quasi-static energies studied here from the potentials that appear in the Schrödinger equation.
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 results of Sect. 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. [21,22]. 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 uncanceled beyond NLO. We use this fact to grasp the dependence of our predictions on higherorder 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 11 Setting ν = 1/r may introduce problems with renormalons (see the discussion in Ref. [49]), 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 (2,0) p 2 depends on the leading pole mass renormalon). Moreover, we are at low orders in the weakcoupling expansion and indications of renormalons, even if existing, cannot be clearly distinguished from other uncertainties. 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 Fig. 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 Fig. 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  Next we consider E b . The situation here is far from optimal, as illustrated in Fig. 6. We have one less lattice 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 Fig. 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 from 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 to E b as far as the number and location of the lattice points are concerned; see Fig. 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 Fig. 7 for x = 2. Setting ν ∼ 1/r does not help much, until we set ν ≈ 1/r , where the NLO result gets quite close to the lattice data (unlike in the previous cases here we can set ν quite low). However, it still shows a huge difference from the LO result. Similarly, the resummation of ultrasoft logarithms does not improve 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 [23]. Nevertheless, in a dedicated series of works [21,22], the SI momentum-dependent quasi-static energies have been computed with increased precision on the lattice and a nontrivial behavior has been found at short distances. A nonzero 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 Fig. 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 Fig. 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 Fig. 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 from 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  Fig. 11 Comparison of the lattice data with the predictions from effective string theory at long distances for E b and E c 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 simu- Fig. 12 Comparison of the lattice data with the predictions from effective string theory at long distances for E d and E e lated 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 momentumdependent quasi-static energies obtained in Ref. [50]. For the SI momentum-dependent quasi-static energies, they are actually equivalent to the predictions obtained first with the minimal area law [18]. 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 poten-tial at long distances. The number we use is taken from [52] 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 Figs. 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 [53]. 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.

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 heavyquarkonium potentials, in terms of Wilson loops (as defined in Sect. 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 we 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. [23], and more recently in Refs. [20][21][22]. We have performed a preliminary comparison with the lattice data [20][21][22] at short distances in Sect. 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 nonzero at O(α 2 ). Earlier lattice simulations predicted an approximately zero slope for them in the short-distance limit. The most recent simulations show a nonzero 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.