Potential NRQCD for unequal masses and the $B_c$ spectrum at NNNLO

We determine the $1/m$ and $1/m^2$ spin-independent heavy quarkonium potentials in the unequal mass case with $\mathcal O(\alpha^3)$ and $\mathcal O(\alpha^2)$ accuracy, respectively. We discuss in detail different methods to calculate the potentials, and show the equivalence among them. In particular we obtain, for the first time, the manifestly gauge invariant $1/m$ and $1/m^2$ potentials in terms of Wilson loops with next-to-leading order (NLO) precision. As an application of our results we derive the theoretical expression for the $B_c$ spectrum in the weak-coupling limit through next-to-next-to-next-to-leading order (N$^3$LO).


Introduction
In analogy to Hydrogen, the dynamical properties of quark-antiquark systems near threshold and with large quark masses (or Heavy Quarkonium for short) can be obtained by solving a properly generalized nonrelativistic (NR) Schrödinger equation. Whereas potential models have been used for years with reasonable phenomenological success, their connection with QCD has always been obscure, to say the least. On the other hand, the use of Effective Field Theories (EFT's), in particular of potential NRQCD (pNRQCD) [1,2] (for reviews see [3,4]), allows us to quantify this connection, and to derive the Schrödinger equation and its corrections from the underlying theory, in a model independent and efficient way. In the extreme weak-coupling limit we will consider in this paper, the EFT can be summarized schematically by This EFT makes the NR nature of the problem manifest and exploits the strong hierarchy of scales that govern the system: where v is the heavy-quark velocity in the center of mass frame.
A key ingredient in the EFT is the heavy quarkonium potential that appears in the Schrödinger equation. It consists of the static potential V (0) at leading order, i.e. O(m 0 ), and relativistic corrections, which are suppressed by inverse powers of the heavy quark masses. 1 The potential is obtained by matching NRQCD [5,6] to pNRQCD. There are several ways to carry out the matching in practice. The most common are i. On-shell matching, ii. Off-shell matching, iii. Wilson-loop matching.
In the on-shell matching one equates S-matrix elements of NRQCD and pNRQCD order by order in an expansion in the QCD coupling constant α and the velocity v (∼ α). 2 The S-matrix elements are defined for asymptotic external quark states satisfying the equations of motion (EOMs) of free quarks. This necessarily requires the incorporation of potential loops, i.e. loops with loop momenta (k 0 , k) ∼ (mv 2 , mv), in both calculations. The reason is that the free-quark on-shell condition produces an imperfect cancellation between potential loops in NRQCD and pNRQCD, and mixes different orders in the 1/m expansion. This obscures the mass dependence of the potential, as it invalidates a strict 1/m expansion for the determination of the potentials, i.e. in the on-shell matching computation the potentials at a given order in 1/m also receive contributions from matrix elements involving operators of higher order. On the other hand, the on-shell matching result for the potential is gauge invariant (to a fixed order in v), as are the S-matrix elements.
In the off-shell matching one equates off-shell Green functions computed in NRQCD with the corresponding off-shell Green functions in pNRQCD (still respecting global energymomentum conservation). In other words, we do not impose that the external quark fields fulfill the free EOMs. This allows us to perform the matching within a strict 1/m expansion, since potential loops in NRQCD and pNRQCD exactly cancel each other. Hence we can keep exact track of the mass dependence of the resulting matching condition for the potentials. The drawback is that the expression we get from the off-shell matching for the individual potentials may depend on the gauge. The total expression for the potential, though, still yields of course gauge invariant results for observables, in particular for the bound state energies, within the accuracy of the computation. In addition, the potentials may acquire some polynomial energy dependence, of which one should get rid by using field redefinitions, or, equivalently, the complete leading order EOMs (including the static potential) if working at lowest order.
In the Wilson-loop matching one equates NRQCD and pNRQCD gauge-invariant offshell Green functions, i.e. Wilson loops (with chromo-electric/magnetic insertions), directly in position space. Working in position space is not the major difference with respect to the previous matching schemes. (Obviously, by Fourier transforming the three-momentum, the on-and off-shell matching computations could also be done in position space.) The key point is that the time of the quark and antiquark fields are set equal. This is not a restriction, and is in fact the natural thing to do for the heavy quark-antiquark system near threshold. We also incorporate gluon strings between the quark and the antiquark fields such that the whole system is gauge invariant. The details of how this matching is performed can be found in Refs. [7,8]. In the static limit, it reduces to the original computation of the static potential by Wilson [9]. The advantage of this procedure is twofold: the matching can be done in a strict 1/m expansion (potential loops do not have to be considered), and closed expressions in terms of rectangular Wilson loops (with chromoelectric/magnetic field operator insertions) can be obtained for each potential. They are therefore explicitly gauge invariant. This makes this procedure quite appealing. In fact the static potential is typically computed this way. We will see that also the relativistic corrections can be efficiently computed using this method.
At present, the heavy quarkonium potential is known with N 3 LO precision (V ∼ mv 5 ) for the equal mass case in the on-shell matching scheme [10]. Nevertheless, there are several reasons why we would like to know the heavy quarkonium potential with N 3 LO precision for the unequal mass case, and also in other matching schemes. Let us highlight two of them: • The B c system The LHC provides a unique opportunity to study the properties of the B c bound states in great detail. In particular, the possibility to measure a good deal of the B c spectrum and decays is now a reality. Obviously, a major ingredient in such analyses is a detailed knowledge of the heavy quarkonium potential and spectrum in the short distance limit. In this paper we calculate both.
• The heavy quarkonium potential in terms of Wilson loops It is possible to give closed expressions for the potentials in terms of Wilson loops that can be generalized beyond perturbation theory. They are therefore suitable objects for the study of nonperturbative QCD dynamics by comparing different models with lattice simulations. (The Wilson loop representation of the potentials indeed allows for exact results in the case of QED, e.g. that the 1/m potential is zero to all orders [7].) For such analyses it is also important to control the short distance behavior of the potentials.
Another important motivation for this paper is to set the ground for higher order computations of the potentials, which we stress again are key ingredients in any observable related to heavy quarkonium we can think of (spectrum, decays, NR sum rules, tt production near threshold, ...). We would like to systematize their computation as much as possible, since, as one goes to higher orders, and as soon as ultrasoft effects start to play a role, the understanding of the relation of the computed potential to the EFT framework becomes compulsory.
In this respect, we believe that it is important to clarify the relation between the different matching schemes and to explore their advantages and disadvantages. The three matching methods mentioned above have been employed more or less independently over the years. The on-shell method has mostly been used to obtain the relativistic corrections to the heavy quarkonium potential [10][11][12][13][14][15]. Earlier, low-order computations, did not require the whole EFT machinery, and some recent computations have profited from the threshold expansion of scattering diagrams [16]. The off-shell method has mainly been used in QED [17,18] but also in some QCD computations [19]. The Wilson loop matching has been the less developed for weak-coupling calculations except for the very relevant case of the static potential [20][21][22][23][24], and the leading, O(α 2 ), contribution to the 1/m potential [7].
The results obtained with these methods are often different, which makes a comparison difficult. On top of that, there is the problem of how to renormalize the potentials in pNRQCD, i.e. how the ultrasoft divergences are subtracted from the bare potentials. There is much freedom here as well. One can perform the subtraction in momentum or position space. In the latter case one can define the subtraction for the potentials in D = 4 + 2 or in four dimensions. These different subtraction/renormalization prescriptions give rise to different expressions for the renormalized potentials (even if all of them only account for soft physics), but not for physical observables. We also note that, while computations using on-shell/off-shell Green functions are naturally done in momentum space, the Wilson loop calculations are naturally carried out in position space (as is the computation of the spectrum). We will discuss these issues in some detail. In particular we will put a special emphasis on matching schemes that admit a strict 1/m expansion of the potential in this paper.
In this work we will focus on the spin-independent potentials. The spin-dependent potentials are not affected by ultrasoft divergences, nor by field redefinitions, to the order required for the calculation of the heavy quarkonium mass with N 3 LO accuracy. Therefore, we will not consider them in detail in this paper and only use known results for the final determination of the B c spectrum. Nevertheless, we will present the spin-dependent results in a form compatible with our EFT computation.
Throughout this paper we will use the abbreviations FG for Feynman gauge and CG for Coulomb gauge.
The outline of the paper is as follows. In Sec. 2 we present the NRQCD and pNRQCD Lagrangians. We also discuss how the potentials are affected by field redefinitions. In Sec. 3 we determine the full D-dimensional result of the O α 2 /m 2 spin-independent potential for different schemes: on-shell, off-shell in CG and FG, and with Wilson loops. In Sec. 4, we obtain the O α 3 /m potential in the different schemes to O( ). In Sec. 5, we present the renormalized potentials. In Sec. 6, we confirm that our expressions comply with Poincaré invariance constraints. Finally, in Sec. 7 we compute the full NNNLO spectrum for unequal masses. We conclude in Sec. 8.
2 Preliminaries: NRQCD and general structure of the potential

NRQCD
The NRQCD Lagrangian is defined uniquely up to field redefinitions. In this paper we use the following NRQCD Lagrangian density for a quark of mass m 1 , an antiquark of mass m 2 (m 1 ∼ m 2 ∼ m Λ QCD ) and n f massless fermions to O(1/m 2 ) [5,6,25,26]: 3 3) 3 We also include the D 4 /(8 m 3 ) terms since they will be necessary later on.

δL
(1) Here ψ is the NR fermion field represented by a Pauli spinor and χ c ≡ −iσ 2 χ * is the respective antifermion field also represented by a Pauli spinor. The matrix (T a ) T is the transpose of the SU (N c ) generator T a in the fundamental representation, and T a → (T a ) T in Eq. (2.8) only applies to the matrices contracted with the heavy quark color indexes. The components of the vector σ are the Pauli matrices. We define iD 0 = i∂ 0 − gA 0 , iD = i∇ + gA, E i = G i0 and B i = − ijk G jk /2, where ijk is the three-dimensional totally antisymmetric tensor 4 with 123 = 1 and (a × b) i ≡ ijk a j b k . For a list of the relevant Feynman rules derived from Eq. (2.1) we refer e.g. to Refs. [4,27]. The relevant NRQCD Wilson coefficients c i and d ij are collected in Appendix B. Throughout this paper we work in the MS renormalization scheme, where bare and renormalized coupling are related as and α = g 2 ν 2 /(4π). In the following we will only distinguish between the bare coupling g B and the MS renormalized coupling g when necessary.

pNRQCD: Potentials
Integrating out the soft modes in NRQCD we end up with the EFT pNRQCD. The most general pNRQCD Lagrangian compatible with the symmetries of QCD that can be con-structed with a singlet and an octet (quarkonium) field, as well as an ultrasoft gluon field to NLO in the multipole expansion has the form [1,2] and M = m 1 + m 2 . We adopt the color normalization for the singlet field S(r, R, t) and the octet field O a (r, R, t). Here and throughout this paper we denote the quark-antiquark distance vector by r, the center-of-mass position of the quark-antiquark system by R, and the time by t. Both, h s and the potential V s are operators acting on the Hilbert space of a heavy quark-antiquark system in the singlet configuration. 5 According to the precision we are aiming for, the potentials have been displayed up to terms of order 1/m 2 . 6 The static and the 1/m potentials are real-valued functions of r only. The 1/m 2 potentials have an imaginary part proportional to δ (3) (r), which we will drop in this analysis, and a real part 5 Therefore, in a more mathematical notation: h →ĥ, Vs(r, p) →Vs(r,p). We will however avoid this notation in order to facilitate the reading. 6 Actually, we also have to include the leading correction to the nonrelativistic dispersion relation for our calculation of the Bc spectrum: that may be decomposed as: where, Note that neither L 1 nor L 2 correspond to the orbital angular momentum of the particle or the antiparticle.
Due to invariance under charge conjugation plus m 1 ↔ m 2 interchange we have This allows us to write Invariance under charge conjugation plus m 1 ↔ m 2 also implies V (2,0) Our aim is to calculate the potentials. In order to do so we can neglect the center-ofmass momentum, i.e. we set P R = 0 in the following and thus L 1 ≡ r × p 1 = r × p ≡ L, L 2 ≡ r × p 2 = −r × p ≡ −L. As explained in the introduction we will not consider the spin-dependent potentials for most of the paper and focus on the spin-independent ones.

Potentials in momentum space
Unlike the position space potential V s , the momentum space potentialṼ s is a c-number, not an operator. It is defined as the matrix element (with P R = 0 from now on) For the static potential we have (k = p − p ) 3 ( ) can be found in Ref. [28]. For the one-loop resultD (0) 2 ( ) we have also done the computation in CG. Throughout this paper we will use the notation For the 1/m potential we follow the standard practice of making the prefactor 1/k explicit: . (2.32) In momentum space we choose the following basis for the 1/m 2 potentials: (1,1) The Wilson coefficientsD (1,1) In our convention the different coefficients of the Taylor expansion are dimensionless except forD (1,1) r,0 ( ). Implicit in the definitions above is the fact that the mass dependence of the potentials admits a Taylor expansion in powers of 1/m 1 and 1/m 2 (up to logarithms). This is so in the off-shell and Wilson-loop matching scheme but not in the on-shell scheme. An exception is againD (1,1) r,0 ( ), since it depends on the NRQCD four-fermion Wilson coefficients, which have a non-trivial mass dependence. 7 We will discuss these issues further in the following sections.

The L 2 operator and potentials in D dimensions
We work with dimensional regularization. Therefore, we need to define the potentials in D = 4+2 dimensions. In the previous section we have given D-dimensional expressions for the potentials in momentum space. In position space, for the spin-independent potentials, everything works as in four dimensions except for the L 2 operator. The definition of the operator L 2 in D dimensions is ambiguous. In this paper we choose the definition (2.41) The right-hand-side of the equation is equal to L 2 r 2 in four dimensions and commutes with pure functions of r in D dimensions, i.e. [f (r), L 2 r 2 ] = 0, as we would expect for an angular momentum operator.

Position versus momentum space
We now proceed to relate the potentials in position and momentum space. For the static and 1/m potentials the relation is straightforward. After Fourier transformation to position space Eq. (2.30) becomes is the d-dimensional Fourier transform of |k| −n .
For the 1/m potential we have To Fourier transform the 1/m 2 potentials some preparation is required. Given two generic functions of r, f (r) and g ij (r) = A(r)δ ij + B(r) r i r j r 2 , the following equalities hold: 8 for arbitrary functions f (r).
Furthermore, we can write 47) and analogously forD (1,1) off . The last equality is especially useful, because the first term has the structure of the Fourier transform of the left-hand-side of Eq. (2.46). It allows us to relateṼ with the potentials in position space: 9 and similarly for V (1,1) off . Note that the three potentials V L 2 , V p 2 and V r receive contributions from the Fourier transform ofṼ off . On the other hand,Ṽ p 2 andṼ r only directly contribute to V p 2 and V r . We stress thatṼ p 2 /r is not the Fourier transform of V p 2 /r .
In summary, we have the following relations: (2.55) 9 For the inverse Fourier transform the following relation is useful: (2.58) (1,1) off In the second equality of each expression we have expanded in powers of g 2 B . Again, F n (r) has been defined in Eq. (2.43) and expressions with F −2 (r) should be treated with care, as such operators are singular.
At each order in g 2 B it is possible to obtain closed expressions relating the 1/m 2 coefficients in momentum and position space. The position space expressions are, however, more complicated than for the static and 1/m potentials. The off-shell potentialṼ off obscures the relation between the momentum and position space potentials. Note also that V r can always be written as [p i , [p i , V r,bis (r)]], where V r,bis (r) has the same dimensions as V p 2 and V L 2 .

Field redefinitions
The bases of potentials, Eqs. (2.44), (2.54)-(2.59), in position space, and (2.32)-(2.34) in momentum space, are ambiguous. There is a large freedom to reshuffle (parts of) some potentials into others using unitary transformations of the pNRQCD fields S and O, which leave the spectrum unchanged. It turns out that we can even eliminate the 1/m potential or, alternatively, the off-shell 1/m 2 potentialṼ off , completely by such field redefinitions. In fact, the latter is achieved in the on-shell matching scheme, which provides us with a minimal basis of operators by construction, as it systematically uses the free EOMs and sets p 2 = p 2 . The drawback is that, as it relies on the free EOMs, the determination of the potentials has to be corrected order by order in α, through potential loops. Still, once a minimal basis is fixed, there is no ambiguity left and each potential is well defined on its own. This can also be seen by looking at the energy shifts, or corrections to the S-matrix, produced by each individual potential in a minimal basis. This also implies that unitary transformations that keep the Hamiltonian in a given minimal basis cannot move terms between the potentials.
In this work, however, we want to keepṼ off , in order to enable a strict 1/m expansion and to maintain the Poincaré invariance relations, see Sec. 6. We are also not particularly interested in completely eliminating the 1/m potential, as it naturally appears in the Wilson loop matching, as well as in the off-shell/on-shell matching schemes.
Instead, the goal of this section is to determine the field redefinitions that translate the results of different matching schemes into each other. This will eventually allow us to combine our calculation of the 1/m 2 potential with the result of the 1/m potential in the on-shell matching scheme for the equal mass case computed in Ref. [15] to obtain the 1/m potential in the unequal mass case in Sec. 4. Following Ref. [7] we proceed as follows. The Hamiltonian has the form where · · · stands for O(1/m) (or higher order) potentials that we are not interested in eliminating.
The unitary transformation More explicitly, under the condition {W, p} m r (which is necessary in order to maintain the standard form of the leading terms in the Hamiltonian, i.e. a kinetic term plus a velocity independent potential) h s reads (2.62) By choosing we completely eliminate δV 1 mr from h s . Moreover, since δV 1 ∼ α 2 (as there is no tree-level 1/m potential), for the precision of the calculations in this paper, we can neglect some terms in Eq. (2.62): (2.64) Therefore, eliminating δV 1 /m r is equivalent to introducing an extra 1/m 2 potential: and Eq. (2.46), we obtain where, without loss of generality, Hence, we find in momentum space (see Eq.(2.46) and following equations) where we have defined This has the important consequence that through O(α 2 ) the coefficientsD r andD p 2 remain invariant under the field redefinitions discussed above. One can also check that the O(α/m 3 ) potential is invariant under the field redefinition (2.64). We will make use of these results in the following. At higher orders in α the neglected terms in Eqs. (2.62) and (2.64) may give an extra contribution toD r . On the other hand, note thatD p 2 is unaffected by the field redefinition, Eq. (2.61), at any order in the α expansion.
Finally, we stress that, since the unitary transformation used in this section can move us into a minimal basis, and, since the static and the α/m 3 potential remain invariant under such transformation, the result for these two potentials is independent of the specific matching scheme used to determine them.
3 Determination of the O(α 2 /m 2 ) potential for unequal masses The spin-independent potential at O(α 2 /m 2 ) for unequal masses is so far unknown. We fill this gap in this section by explicitly calculating it in different matching schemes and with full dependence. Our results directly fix the bare coefficientsD in each case. With little effort and using the equations in Sec. 2.2.3, one can then obtain the expressions for the bare coefficients D of the potential in position space. Note that all 1/m 2 position space potentials depend on the matching procedure (albeit some of them weakly, in the sense that the matching scheme dependence vanishes when → 0), as do g off and h off in Eqs. (2.54)-(2.59). Therefore, instead of presenting explicit expressions, we give the position space results only in terms of the momentum space coefficients in Sec. C.

Matching with Green functions
In the off-shell matching we equate four-point off-shell Green functions computed in NRQCD with the analogous four-point off-shell Green functions in pNRQCD. In this way we will determine the 1/m 2 pNRQCD potential at O(α 2 ). We do not require the quarks to fulfill the free EOMs, i.e. the only restriction on the external momenta is total energy-momentum conservation. This allows us to perform the matching in a strict 1/m expansion, since NRQCD and pNRQCD potential loops cancel each other exactly. Hence, we can directly equate soft NRQCD diagrams (computed with static quarks) with the bare potentials in pNRQCD at a given order in 1/m.
By contrast, in the on-shell matching S-matrix elements of NRQCD and pNRQCD are equated order by order in an expansion in α and v (∼ α). These S-matrix elements are computed with asymptotic quarks satisfying the free EOM. This necessarily requires the incorporation of potential loops in both calculations, since the on-shell condition causes an imperfect cancellation between potential loops in NRQCD and pNRQCD. The latter mix different orders in the 1/m expansion, i.e. potential loops involving a potential at a given order can contribute to the matching of a potential at lower orders. See, for instance, Ref. [11] for an illustrative example.

Off-shell matching: Coulomb gauge
In Refs. [17,18] the off-shell matching between NRQED and pNRQED has been studied in detail with O(mα 5 ) precision in CG. The FG matching has also been discussed in Ref. [18] with O(mα 4 ) precision.
We now perform the matching for the case of QCD. We focus on the relativistic 1/m 2 corrections to the potential. The tree-level matching is analogous to the one in QED up to the straightforward incorporation of color factors: The gauge-dependent off-shell coefficientsD are given here in CG and labeled accordingly. Now we consider the one-loop corrections. In Appendix E we present the result of the (sum of the) relevant diagrams in CG as well as in FG. It has always been assumed that the evaluation of Feynman diagrams in the CG can be quite cumbersome, especially for non-Abelian gauge theories. We find that this is not the case, at least for the computation we perform in this paper. More details on the computation will be shown in Ref. [29].
The diagrams depend on the energies of the four external quarks E i , see Appendix E. This dependence can be eliminated in the potentials using field redefinitions. Their implementation can however be cumbersome. Fortunately, for our purposes it is not necessary. As discussed in Appendix E, at the order we are working, we should use the complete EOMs, which include the static potential. Effectively, though, we can neglect the static potential, as the difference contributes to the 1/m potential, which we will determine in an independent way, anyhow. This is equivalent to using the free EOMs, but still keeping p 2 = p 2 for the incoming and outgoing quark momenta (in the center of mass frame), unlike in the on-shell matching. Finally, we obtain the following (bare) CG results where a 1 and β 0 are defined in Sec. A. Note that, strictly speaking, there are subleading contributions in powers of α encoded in the NRQCD Wilson coefficients.

Off-shell matching: Feynman gauge
The matching in FG involves considerably more (soft) NRQCD diagrams. In particular, diagrams with only A 0 gluon exchanges now give a nonzero contribution. As a consequence, the dependence on the (off-shell) external quark energies is more complicated. The complete expression for the sum of all one-loop diagrams can be found in Appendix E. Yet, after using the free EOMs (which is sufficient at the order we are working at) we find that the coefficientsD r andD p 2 agree with their CG results. This is indeed what we expected, as these potentials remain the same in the on-shell limit. The differences to CG therefore manifest themselves only in theD off coefficients. At tree level in FG (and at one loop in the CG) an energy-dependent term ∝ k 2 0 = (E 1 − E 1 ) 2 occurs. In principle, the redefinition of the quark energies in terms of threemomenta is ambiguous. In this special case, however, there is a preferred prescription (see Ref. [18]) to transform away the energy dependence, namely Eq. (E.7). It is the only way to preserve the 1/(m 1 m 2 ) structure and at the same time leave the 1/m potential unchanged, see Appendix E for details. Adopting this prescription we arrive at the same result as in CG:D With the energy replacement rules given in Appendix E we obtain at one loop

On-shell matching
Finally, we determine the potential in the on-shell matching scheme. In this scheme we haveD off,on−shell = 0 by construction. At the order we are working at, this means It turns out that for the other potentials a dedicated on-shell matching computation is not necessary. A priori, we must take into account potential loops, which are not needed in the off-shell computation, in addition to the soft NRQCD loops. The discussion on field redefinitions in Sec. 2.3 however shows that the transformation from an off-shell to the on-shell scheme leaves the coefficientsD p 2 andD r , as well as the O(α/m 3 ) potential, unchanged at the order we are working at. Hence, potential loops can neither contribute toD p 2 ,2 andD r,2 , 10 nor to the O(α/m 3 ) potential. Therefore, these coefficients are equal irrespectively of computing them on-or off-shell, and in the latter case they are independent of the gauge, as we have seen. This is in fact the reason why we have not labeled them according to the matching procedure.
For equal masses and in the on-shell matching scheme the potential has been computed in Refs. [10][11][12][13][14]. In particular, we compared our results with the ones of Refs. [10] and [14]. The complete dependence for the equal mass case can be found in Ref. [27]. We agree with their results. The novel results of the present section are the potentials for unequal masses (keeping track of the NRQCD Wilson coefficients).
As another cross check we have calculated the O(α 2 /m 2 ) potential for unequal masses from soft on-shell scattering amplitudes in vNRQCD [30] using the Feynman rules given in Ref. [14]. We found complete agreement with our momentum space results in the on-shell matching scheme.
The momentum space results of sections 3.1.1 -3.1.3 can straightforwardly be transformed to position space using Eqs. (2.54)-(2.59). In Appendix C we give the corresponding expressions for all schemes in terms of the respective momentum space coefficients.

The quasi-static energy and general formulas
An alternative determination of the potentials is the direct matching of NRQCD and pNRQCD gauge-invariant Green functions in position space. One key point is that the time of the quark and antiquark are now set equal. This is not a restriction. Instead, it is rather natural to describe quark-antiquark bound states by fields that depend on a single time coordinate. Another difference to the off-shell matching scheme is the insertion of gluon strings (Wilson lines) between the static quark and antiquark in order to form a Wilson loop, so that the whole system is gauge invariant.
The details of the Wilson-loop matching procedure are given in Refs. [7,8]. In these references the emphasis was put on the matching in the nonperturbative scenario without ultrasoft degrees of freedom. Two alternative methods were worked out in detail. One is the direct matching between NRQCD and pNRQCD Wilson loops, and the other one is a generalized "quantum-mechanical" matching, which gives the spectral decomposition of the potentials, allowing them to be rewritten in terms of Wilson loops. Either way, the matching can be done in a strict 1/m expansion (potential loops do not have to be considered at all) and closed expressions in terms of Wilson loops can be obtained for each potential, which are then manifestly gauge invariant. This allows for a nonperturbative definition of the potential E s , to which we will refer to as the "quasi-static" energy in the following. Formally we write We use "E" to make the distinction to the potentials "V " explicit. The latter are, by definition, the potentials of the Schrödinger equation. In the strong-coupling regime (and provided there are no ultrasoft degrees of freedom), the "quasi-static" energy replaces the potential in the Schrödinger equation describing the nonperturbative heavy quarkonium bound state. Once ultrasoft effects are included (as e.g. in our calculation of the B c spectrum) this is not true anymore. That is why we distinguish explicitly between E and V . We will elaborate on this in Sec. 3.2.2 and in a forthcoming paper. We use the following definitions for the Wilson-loop operators. The angular brackets . . . denote the average value over the Yang-Mills action, W is the rectangular static Wilson loop of dimensions r × T W : and . . . 20) and similarly with extra operator insertions (see Ref. [8] for more details). At leading order in the 1/m expansion, we get nothing but the static energy already found by Wilson many years ago [9] The complete expression of the 1/m and 1/m 2 potentials in the quenched approximation (no light quarks) in terms of Wilson loops has been determined in Refs. [7,8] (partial results for the 1/m 2 potential can be found in Refs. [31][32][33][34]). For these we define the shorthand notation where T W is the time length of the Wilson loop and T is the time length appearing in the time integrals shown below. By performing the limit T W → ∞ first, the averages . . . Refs. [7,8], with the exception that we rewrite some of them so that they remain valid in D dimensions. For the spin-independent potentials we have where in the second-to-last line the light-quark operators are located in the heavy-quark Wilson line (i.e. at the position x 1 ). The last term contains the 1/m 2 operators in the NRQCD Lagrangian that only involve light degrees of freedom. Note also that the other Wilson-loop expectation values should be computed with dynamical light quarks. Equation (3.28) generalizes the result of Ref. [8] to the case with light fermions (as usual neglecting ultrasoft effects). Note that, although, formally, the first, the second-to-last, and the last lines of Eq. (3.28) depend on the time where the operators are inserted on the heavyquark lines, this is not so after performing the T W → ∞ limit, due to time translation invariance. Finally, the last term we need in Eq. (3.18) is 11 Let us further elaborate on the expressions for E (2,0) r and E (1,1) r . The first term of 11 The first term of this equation corrects a sign error in the first term of Eqs. (48) and (54) in Ref. [8].
Note that its spectral decomposition in Eq. (23) of that reference is correct though.
(2,0) r admits the alternative representation It is also possible to use the Gauss law to simplify Eq. (3.30). This was done in Ref. [8] for the case without light fermions. Including them we find i γ 0 T a q i everywhere. This makes their dependence on the light fermions more explicit. We obtain where the last six lines are due to light fermions, and the light quark operators are located on the heavy quark Wilson line (i.e. at the position x 1 ) except for the last operator, and In summary, the results of this subsection are the generalization of the results of Ref. [8] for Once we focus on the weak-coupling regime, ultrasoft degrees of freedom certainly contribute to the quasi-static energies. They do so with energies/momenta of order ∆V ≡ V Nevertheless, for a consistent description of the weaklycoupled quark-antiquark system, the potentials in the Schrödinger equation (i.e. in the pNRQCD Lagrangian) should only include contributions associated with the soft modes. Taylor expanding in powers of 1/m before integrating over the gauge or light-quark dynamical variables effectively sets the potential loops to zero. However, this does not eliminate the ultrasoft contributions from the potential expressed in terms of Wilson loops. Actually, as far as the ultrasoft modes are concerned, the 1/m expansion can be formally understood as exploiting the hierarchy ∆V p 2 /m, which is the limit implicit in the discussion of Sec. 3.2.1. 12 In order to obtain the potentials in perturbation theory, the ultrasoft contribution has to be subtracted. This can be easily achieved by expanding in the ultrasoft scale before performing the loop integration. Thus, only the soft scale appears in the integrals, which 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 are interested in). In summary we have, cf. Eqs. (3.18) and (2.11), where we have put the subscript W to indicate the Wilson-loop matching scheme. For the static potential V (0) , the Wilson loop definition is given in Eq. (3.21). Its perturbative evaluation in powers of α can be transformed into a calculation in momentum space, where the energies of the external quark and antiquark are set to zero for T W → ∞, since the time-dependent part of the external quark propagator, θ(T W − t), can be approximated by 1. See also Ref. [28] for a detailed discussion. In addition, for a certain class of gauges (including FG and CG), one usually neglects the exchange of asymptotic gluons from the boundaries of the Wilson loop at ±T W /2 for T W → ∞, see the discussion in Refs. [20,28]. In this setup, the Wilson-loop matching for the static potential is equivalent to a standard diagrammatic S-matrix calculation with off-shell static quarks, i.e. with zero (kinetic) quark energies, but nonzero external three-momenta. This is indeed equivalent to the off-shell matching computation at leading order in the 1/m, E 1 and E 2 expansion.
On the other hand, at lowest order in 1/m no kinetic propagator insertions are involved in a soft NRQCD S-matrix calculation, as they would inevitably come with factors of 1/m. 12 Obviously, this is not the kinematic situation we face in the bound state, where ∆V ∼ p 2 /m. It therefore does actually not matter for the latter calculation, whether the external quarks are on-or off-shell. Furthermore, potential loop contributions to the static potential in the on-shell matching scheme must vanish, because there are no field redefinitions (compatible with the symmetries of QCD) that could possibly remove them by modifying a higher order potential, cf. Sec. 2.3. Hence, we conclude, that the static potential is the same in any of the matching schemes discussed in this paper.
The Wilson-loop calculation for the higher-order potentials cannot be related to a purely momentum space S-matrix calculation due to the insertions of gluonic/light-quark operators that are integrated over time. Nevertheless, we will see that we can also compute the higher-order potentials in the Wilson-loop matching scheme efficiently based on Feynman diagrams. It is worth emphasizing that the expressions for the potentials in terms of Wilson loops encapsulate all effects at the soft scale in a compact way, and they are correct to any finite order in perturbation theory. In particular, compared to the standard calculation of the static potential, only a few extra Feynman rules for the operator insertions have to be introduced (see Appendix F) once the exchange of asymptotic gluons from the boundaries of the Wilson loop at ±T W /2 for T W → ∞ is neglected. This is to be contrasted with the matching of Green functions, where higher-order kinetic insertions on the propagators must be taken into account, both for on-shell and off-shell matching, which can be quite tedious at higher orders. When matching on-shell, in addition, potential loops must be considered. Let us now compute V (2,0) We use this case in order to illustrate how we perform the Wilson loop calculations.
At O(α) we only have contributions to V (1,1) The diagrams needed are drawn in Fig. 1. In CG only the second diagram contributes and using the Feynman rules derived in Appendix F, the detailed calculation reads where the projector P ij (k) = δ ij − k i k j k 2 . In FG both diagrams in Fig. 1 contribute, but we still obtain the same result, as expected due to gauge invariance of the Wilson loop. At this order the result coincides with the result obtained using off-shell matching. 13 Thereforẽ where, D = d + 1. Here we chose the energy l 0 to flow along the arrow between the crossed vertices in Fig. 2 and the momentum q to flow counter-clockwise in the loop. The (integrand of) the one-loop amplitude M ij can be obtained by applying standard static Wilson-loop Feynman rules together with the additional rules for the E i operator insertions as given in Appendix F. Note that we have pulled out a factor 1/(l 0 + i0), corresponding to the upper static quark propagator from the amplitude's integrand, in order to render M l 0 -independent. We emphasize that, as in the calculation of soft on/off-shell Green functions, we must neglect the (ill-defined) contribution of pinch singularities. The latter are related to iterations of lower order potentials and are not part of the soft regime. In fact, the pinch 13 In CG it coincides exactly, in FG only after using the EOMs as discussed above. singular terms are explicitly removed in the definition of connected Wilson loops according to Eq. (3.20).
In CG only the first two diagrams of Fig. 2 contribute (the first gives a divergent contribution). Using our Wilson-loop Feynman rules in Appendix F we find the (unintegrated) amplitude Plugging this in Eq. (3.38) gives We have also checked that we get the same result performing the calculation in FG, where all four diagrams contribute. Note that V L 2 ,W is carried out along the same lines. The diagrams contributing at O(α 2 ) are displayed in Fig. 3. In order to have a cross check we compute again in both, CG and FG, and indeed obtain the same result: The above calculations of the V L 2 ,W potentials to O(α 2 ) are actually all we need to fix also the other spin-independent position space potentials V p 2 ,W and V r,W with O(α 2 ) precision. The reason is that we can use Eqs. (2.54) and (2.57) to determine g off,W and then (by inverse Fourier transformation)D off,W (k) in momentum space. We find where and As a check, we can also directly compute V which yields the same result as Eq. (3.45). On top of that, we have also checked that we obtain the same result in FG. We have also computed V (1,1) p 2 (r) directly in CG and FG finding agreement with Eq. (3.48). This is an even stronger check, because the calculation is more difficult, as it involves more diagrams.
From the above analysis we can also determine (the soft part of) some Wilson loops that contribute to V r,W and can be treated separately, because they are multiplied by different NRQCD Wilson coefficients. Let us first focus on V We can also directly compute the Wilson loop and check this result. The relevant diagrams are the same as in Fig. 2. An analogous calculation for the V L 2 ,W potentials yields which is equal to Eq. (3.52). Using this result and Eq. (3.30) we obtain from the comparison to Eq. (3.32). Note that this results from a nontrivial cancellation of non-Abelian contributions so that only light-quark effects survive. This is precisely what should happen according to Eq. (3.32). We can also confirm confirm Eq. (3.54) by direct inspection of the c D +c hl 1 term of V (2,0) r,W (but now written in terms of light-quark operators), which, thus, provides us with an independent check.
Finally, by comparing the terms proportional to c 2 F we find With this we have already exhausted all contributions to V (2,0) r,W . Therefore, we conclude that all the remaining terms are O(α 3 ), i.e., Unfortunately a similar analysis for V

Determination of the O(α 3 /m) potential for unequal masses
The O(α 2 /m) potential for the unequal mass scheme was computed first in the on-shell matching in Ref. [11]. The D-dimensional expression in the same matching scheme, but for the equal mass case, can be found in Ref. [35]. The O(α 3 /m) potential for equal masses was obtained in Ref. [15] using on-shell matching and the O( ) piece can be found in Ref. [36]. Overall, the equal mass result (to the highest order in presently known) reads Note that, unlike the expressions for the potentials in the previous sections, we have written the potential in Eq. (4.1) in terms of the MS renormalized coupling g 2 evaluated at the scale ν (see Eq. (2.10)), 14 because this allows for an easier comparison with the results of Ref. [15].
It is the aim of this section to obtain the expression for the 1/m potential in the unequal mass case for the on-shell, off-shell (CG and FG) and the Wilson-loop matching schemes described in Sec. 3. We will rely on the 1/m 2 results obtained in Sec. 3, as well as on the results of Ref. [15]. A key point in our derivation will be the use of the field redefinitions discussed in Sec. 2.3.
Based on these field redefinitions we have argued in Sec. 3 that through O(α 2 ) the potential coefficientsD p 2 andD r are the same in all three matching schemes. We have checked this prediction explicitly for on-shell and off-shell matching. We have also determined all other 1/m 2 potentials at O(α 2 ). Our results of Sec. 3 thus represent the complete O(α 2 /m 2 ) potential in the Wilson-loop, off-shell and on-shell scheme.
The scheme differences can be compactly expressed in momentum space: . We emphasize that this is possible because we know the complete offshell 1/m 2 potential. In addition, our 1/m 2 results contain the full information on the m i dependence of the potentials. Therefore, we are also able to determine the O(α 3 /m) potential for unequal masses in the on-shell matching scheme, as we will see below.
We start with the results in the Wilson-loop scheme, where the appropriate field redefinition gives m δṼ with . (4.7) We can now use Eq. (4.5) to determineṼ (1,0) W . We find the following momentum space coefficients according to Eq. (2.32): Note that these coefficients refer to the expansion of the 1/m potential in powers of g 2 B . After Fourier transformation to position space we obtain , it is straightforward to identify the field redefinitions that relate the potentials obtained in the Wilson-loop and the CG/FG off-shell matching schemes. We emphasize that the differences δṼ (2) W − δṼ (2) CG and δṼ (2) W − δṼ (2) FG are precisely of the form of Eq. (2.69) with a mass-independentg(k). Hence, according to the field redefinition in Eq. (2.64), the corresponding differences in the 1/m potential are proportional to 1/m r = 1/m 1 + 1/m 2 . This explicitly verifies that the strict 1/m expansion also holds for the CG and FG off-shell schemes.
We now give expressions for the 1/m potentials in the latter schemes. The CG/FG coefficients in Eq. (2.32) read (4.13) In the CG computation it is easy to see that there are no C 2 F T F n f contributions to the 1/m potential by inspection of the possible diagrams at O(α 3 ).
Furthermore, as stated above, we can determine the NLO 1/m potential in the on-shell scheme for unequal masses. Now, however, the field redefinition relating it to the off-shell potentials induces a non-trivial dependence on the masses, because 2D (2,0) off,X = D (1,1) off,X in Eq. (4.4). We obtain  We remark that in the first equality we keep the complete dependence of the terms proportional to the color factors C 2 F C A and C 2 F T F n f . This is an outcome of our calculation. In the second equality we expand to O( ).
Finally, note that, unlike for the off-shell and Wilson-loop potentials, it does not make sense to define V

Renormalized potentials
So far we have obtained the bare potentials for different matching procedures. The different results can be related by unitary field redefinitions. Therefore, the physical spectrum of the quark-antiquark system will be the same irrespectively of the matching scheme used to determine the potentials. 15 In order to produce physical results one always has to add the ultrasoft contribution to the respective observable. The ultrasoft calculation relevant for the determination of the B c spectrum yields the following contribution to the (singlet) heavy quarkonium self-energy (in the quasi-static limit) [17,22,37]: . In general, ultrasoft contributions will depend on the basis of potentials used, but, up to the order we work at here, it only depends on the static octet potential, which is not affected by the field redefinition in Eq. (2.63).
The (ultraviolet) divergences of Eq. (5.1) that are associated with the pole of the heavy quarkonium propagator (i.e. those independent of h s − E) should cancel the divergences of the bare potential V s . We collect the latter in δV s : so that V MS s produces finite physical results. This does not necessarily mean that V MS s is finite in the four-dimensional limit, as the cancellation of divergences should only occur in physical quantities and not necessarily for each individual potential separately.
Let us elaborate on this point. We take Eq. (5.1) and move one factor of (h s − E) to the left, one to the right, and the remaining one in is moved such that one obtains a (h s − E)-free divergence that is cancelled by the counterterm (note that V L 2 does not appear in this expression): This expression was used in Ref. [4]. It is however not unique. If we take (h s −E) 3 and move one factor of (h s −E) to the left, one to the right, and the remaining one is split in half and symmetrically moved to the left and right in Eq. (5.1), we obtain there is no cancellation of the divergences: We cannot get finite four-dimensional expressions for the potentials. Formally this is not a problem, because the uncanceled divergences vanish in the calculation of the spectrum, but we are then forced to compute intermediate results in D dimensions. In practice, it is therefore convenient to find finite renormalized expressions that allow us to work in four dimensions. This is achieved by subtracting δV (W ) s from V s,W and δV (GF) s from the bare potentials in the CG/FG off-shell and on-shell schemes. Finally, using Eq. (2.10), 16 See also the discussion in Ref. [38].
that V A = 1 with leading logarithmic accuracy [39], and the relation between the bare and renormalized expressions of the NRQCD Wilson coefficients presented in Sec. 2.1, we obtain the renormalized potentials for the different matching prescriptions.
In order to simplify the notation we drop the index MS of the NRQCD Wilson coefficients in the expressions of the renormalized potentials we give below. Note also that the divergences of the bare NRQCD Wilson coefficient d sv we use in this paper (computed in FG), do not cancel the divergences of V . On the other hand, had we computed the NRQCD Wilson coefficients in CG, we would find no mixing between these potentials for the cancellation of divergences. See the renormalization group equations in Ref. [40] for the latter case.
We now list the final expressions for the renormalized potentials obtained in the different matching schemes in position space. In the off-shell CG scheme they read In the off-shell FG scheme, we have (5.14) The renormalized potentials obtained from the Wilson-loop prescription are Finally, we present the renormalized potentials in the on-shell scheme: We remark again that in Eqs. Any of the above sets of potentials produces the same spectrum. We also stress that our renormalization procedure does not just subtract the 1/ poles, but also adds some finite pieces and an dependence to the renormalized potentials. We do this in such a way that the ultrasoft bound state calculation is simplified, see Sec. 7.2.
The above renormalized potentials can be transformed to momentum space. We display the resulting expressions in Appendix D.

Poincaré invariance constraints
Poincaré invariance (of full QCD) poses constraints on the form of the heavy quark potential. In the context of our computation the following two relations can be derived Note that they do not involve the NRQCD Wilson coefficients. These relations were originally found in Ref. [33] by explicit calculation of the potentials in terms of Wilson loops. In the context of pNRQCD, and explicitly using the Poincaré algebra, they were deduced in Refs. [41,42]. We have checked that our results fulfill these equalities: We have explicitly verified that Eqs. (6.1) and (6. On the other hand, we stress that the Poincaré invariance constraints cannot be applied to the results in the on-shell matching scheme. The reason is that Eqs. (6.1) and (6.2) are derived assuming a certain mass scaling of the potentials. This assumption does not hold for the potentials obtained by on-shell matching, as the latter mixes different orders in the 1/m expansion.
Finally, it is easy to see that the above Poincaré invariance relations are not affected by our field redefinition in Eq. (2.64), as the latter produces shifts of the form δV

The B c mass to N 3 LO
We are now in the position to compute the spectrum of a heavy quarkonium bound state made of two heavy quarks with different masses with N 3 LO accuracy in the weak coupling limit. We have derived explicit expressions for the (spin-independent) relativistic corrections to the potential. For ease of reference, we quote the known expressions for the renormalized static potential and the spin-dependent potentials in Sec. 7.1. They are not affected by possible ambiguities due to field redefinitions of the kind discussed in Sec. 2.3 to the order we are working at. The static potential is, however, affected by ultrasoft divergences, which we renormalize following the discussion of Sec. 5. In Sec. 7.2 we quote the energy shift produced by the ultrasoft contribution. In Sec. 7.3 we quote the energy shifts associated with the static potential. In Sec. 7.4 we compute the energy shifts associated with the relativistic corrections to the potential, and in Sec. 7.5 we present our final expression for the heavy quarkonium mass.
Eqs. (7.3) and (7.6)-(7.8) correct misprints in Eqs. (70) and (71) of Ref. [4] (when setting the masses equal). For the spin-dependent potentials, in this paper, we can work with the four dimensional expressions for L · S i , S 12 and S 2 . Even though the (soft) matching calculation for these spin-dependent potentials exhibits ultraviolet divergences, they do not require renormalization in pNRQCD. The divergences exactly cancel the ones of the NRQCD Wilson coefficients, so that the overall spin-dependent potential in pNRQCD is finite (to the order of interest), cf. Eqs. (7.3)-(7.8).
The spin-dependent potentials are unambiguous (at least to the order we are working at). They were originally computed in Ref. [11] at NNLO, in Ref. [43] for the N 3 LO hyperfine splitting, and in Ref. [12] the complete expression for unequal masses was obtained.

The ultrasoft energy correction
Combining the results given in Refs. [18,19,37] we find for the ultrasoft contribution to the energy: and Numerical determinations of these non-Abelian Bethe logarithms were obtained for low values of n in Ref. [37] for l = 0 and in Ref. [44] also for l = 0.

Energy correction associated with the static potential
This contribution we extract from the results of Ref. [44]. We (partially) adopt their notation in the following. For the ground state and first excitations the contribution was computed in Refs. [45][46][47]. It follows from standard (time-independent) quantum mechanical perturbation theory up to third order and reads δE(n, l, s, j) where and where 20) and X LS , X LS i , D S and S 12 have been defined in Appendix G.
By default we will use the on-shell potential for the computation, as it will ease the comparison with other results, in particular those of Ref. [44]. We split the computation of the N 3 LO correction to the bound state energy into a spin-dependent and a spin-independent part. The spin-dependent contribution can be organized as follows: Using the expectation values given in Appendix G for single and Appendix H for double potential operator insertions we find The terms ξ SD i are given in Appendix I. For the spin-independent part of the energy we proceed in the same way. In this case the energy shift can be written as (7.25) Using again the expectation values in Appendix G for single and Appendix H for double potential insertions we obtain The terms ξ i and ξ SI i are given in Appendix I.

The O mα 5 spectrum for unequal masses
Summarizing the previous results, we can present the complete expression for the energy levels of a heavy quark-antiquark bound state with unequal quark masses and N 3 LO accuracy: where c i = c c i + c nc i . We have checked that for the ground state the result agrees with the NNLO B c energy given in Ref. [48]. For arbitrary quantum numbers the NNLO result can be found in Ref. [49] (though in a basis different from ours), and in Ref. [50] for the equal mass case. We have also checked that our results agree with the N 3 LO energy in the equal mass case, which was obtained in Ref. [45] for the ground state, in Refs. [46,47] for S-wave states, and in Ref. [44] for general quantum numbers. We also agree with the numerical results given in Ref. [51].
All relevant definitions for the functions and parameters in the previous formulae can be found in Appendix I.

Conclusions
In this paper we have computed the O(α 2 /m 2 ) contribution to the heavy quarkonium spinindependent potential in the unequal mass case. We have obtained the bare expressions (in D = 4 + 2 dimensions) for different matching schemes in momentum and position space. We have performed all our calculations in Coulomb and Feynman gauge. Perturbative evaluations of loop diagrams in Coulomb gauge have always been thought to be complicated and difficult to handle, especially for non-Abelian theories. On the other hand one typically has to compute less diagrams in that gauge. For the one-loop calculations (using dimensional regularization) carried out in this paper, Coulomb gauge has proven to be a competitive method.
In momentum space, the results are encoded in the coefficientsD 2 . The coefficients D These results, obtained from different matching procedures, can be related by field redefinitions. We have identified the field redefinitions that relate the O(α 2 /m 2 ) heavy quarkonium potentials in the different matching schemes. These field redefinitions are valid in D dimensions and can be applied to the bare potentials.
Our calculation yields an independent determination of the bare O(α 3 /m) potential proportional to the color factors C 2 F C A and C 2 F T F n f for unequal masses and for the different matching schemes considered in this paper. For the equal-mass on-shell case it agrees with the results of Ref. [15] up to O( ), but we remark that we also predict the complete dependence of these terms. Using the equal-mass on-shell result of Ref. [15] together with our new O(α 2 /m 2 ) potentials we have determined the other terms of the O(α 3 /m) potential (proportional to C F C 2 A and C F C A T F n f ) for unequal masses and the three different matching schemes to O( ).
For the 1/m potential in terms of Wilson loops we summarize our results in Eq. (4.10), and the corresponding momentum space coefficientsD can be found in Eqs. (4.11)-(4.13). In Eq. (4.14) we present the position-space expression for the unequal-mass 1/m potential in the on-shell scheme (note the non-trivial mass dependence). In the latter case it is actually meaningless to define the coefficientsD (1,0) , as they would depend on the heavy quark masses.
We remark that, in the Wilson-loops scheme, the terms of the O(α 3 /m) potential proportional to the color factors C 2 F C A and C 2 F T F n f vanish. For the Coulomb/Feynman gauge off-shell matching the C 2 F T F n f term is zero, whereas in the on-shell scheme all possible color structures contribute. This suggests that using Wilson loops might be the optimal setup to determine the 1/m potential.
In summary, we have obtained the bare heavy quarkonium potential for unequal masses with the required precision to compute the B c mass with N 3 LO accuracy. We have determined the renormalized potentials in the different matching schemes in Sec. 5 and discussed their dependence on the specific ultrasoft subtraction scheme. We have seen that the relativistic potentials obtained in the Wilson-loop and off-shell matching schemes (both the renormalized and bare expressions) satisfy certain constraints due to Poincaré invariance, unlike those obtained in the on-shell matching scheme.
We have performed the computation of the B c mass with N 3 LO accuracy for arbitrary quantum numbers in Sec. 7. The final theoretical expression is given in Eq. (7.27). Note that, even though the expressions have been obtained in the weak-coupling limit, one can easily obtain expressions valid for m r α 2 ∼ Λ QCD by subtracting the perturbative expression of the ultrasoft contribution, Eq. (7.10), and adding the corresponding expression in that regime (which then includes nonperturbative effects). A phenomenological analysis will be carried out elsewhere.
Other important results of our computation are the NLO expressions for the soft contribution of the 1/m and spin-independent (and velocity-dependent) 1/m 2 "quasi-static" energies in the short-distance limit. These "quasi-static" energies represent non-perturbative definitions of the heavy quarkonium potentials. At this order, the "quasi-static" energies start to be sensitive to ultrasoft effects. Therefore, our results are, in fact, factorization scale dependent. To obtain "physical" results that can be compared with Monte Carlo lattice simulations, like those performed in Refs. [52][53][54], the ultrasoft contributions to the relevant Wilson loops must be computed and added to the results of this paper. This calculation will be carried out in a forthcoming publication.
The analysis of this paper allows us to grasp the advantages and inconveniences of each matching scheme for perturbative evaluations of the potential. As a matter of fact, we find that all methods appear to be feasible in practice. In particular we found that perturbative computations using Wilson loops are not only feasible but may even have some advantages: The potentials in terms of Wilson loops encapsulate in a compact way all the information related to the soft scale, they are correct to any finite order in perturbation theory, and neither kinetic operator insertions nor potential loops have to be considered in the computation (which otherwise can be quite cumbersome at higher orders). We emphasize that, in the case of pure QED, it is possible to obtain closed expressions for some potentials, so that only a few orders in perturbation theory contribute. This implies some all-orders non-renormalization theorems (for the QED part) and, thus, constrains also the ultrasoft contributions.

A Constants and useful Formulae
(A.7)

B NRQCD Wilson coefficients
We have c (i) F − 1 due to reparameterization invariance [25]. In Ref. [55], c F was computed with NLO accuracy. The other NLO Wilson coefficients to O(1/m 2 ) were computed for the one and zero heavy-quark sector in Ref. [25] and for the two heavy-quark sector in Ref. [56], both in FG. Here we only list the Wilson coefficients that are directly relevant for our analysis. 17 Their bare expressions read as this equation corrects Eq. (218) in Ref. [4].
The color constants T F , C A , C F and the QCD β-function coefficients (β i ) are given in Appendix A. In the MS scheme the respective renormalized Wilson coefficients of the single quark sector are (for m j = m i ) 18 The four-quark Wilson coefficients for unequal masses are given by (note that for the equal mass case the annihilation contribution should be included, see Ref. [56] for the specific expressions): (B.12) 18 The term in cD proportional to TF does not appear in the result quoted in Ref. [25]. It is generated by the field redefinition that eliminates the operator GD 2 G from the NRQCD Lagrangian, see the discussion in Ref. [8].
At the order we are working, we can set c (i)hl 1 = 0. However, if we are interested in the resummation of large logarithms, we must keep c (i)hl 1 due to its non-trivial RG evolution. For future purposes we will therefore retain the contribution proportional to c (i)hl 1 in the potential and only set it to zero in the final determination of the heavy quarkonium mass with N 3 LO accuracy.
Since the basis of operators is not minimal, there are ambiguities in the values of some Wilson coefficients. In particular the expressions of d vs and c D depend on the gauge used to determine them (not only the finite pieces but also the logarithmic divergences, see the discussion in Ref. [40]). The expression we give here is the FG result. Also, as we have already mentioned, there is an ambiguity on how the Pauli matrices are treated in D dimensions, affecting the coefficient d vv . Here, we choose the prescription used in Ref. [56]. This will also affect the soft computation of the potentials.  where X is "CG"/"FG" for the CG/FG off-shell matching scheme, "W " for the Wilson-loop scheme, and "on-shell" for the on-shell scheme.

D Results for the renormalized potentials in momentum space
Upon Fourier transformation of Eqs. (5.5)-(5.26) and according to the definitions in Sec. 2.2.1 we obtain the following expressions for the coefficients of the renormalized potentials in momentum space: The coefficientsD off andD (1,0) depend on the matching scheme. For the cases considered in this paper we find Finally, in the on-shell scheme (where obviouslyD (2,0),MS off,on-shell (k) = 0), we havẽ Note that the above expressions are not what one obtains by just subtracting the 1/ poles in momentum space (which is what it is usually named MS scheme). The latter renormalization scheme would complicate the (ultrasoft part of the) bound state calculation for the spectrum. For our purposes, it is more convenient to do the subtraction in position space, and the prescription we have proposed here is particularly useful, because it avoids spurious logarithms of k 2 . We will refer to it as MS scheme in this paper. In this way we can efficiently carry out the bound state computations in four dimensions.
In CG the result reads In order to obtain energy-independent potentials we have to express the energies E i , E i and k 0 in Σ (2) in terms of three-momenta. We achieve this by redefining the heavy quark fields in the pNRQCD Lagrangian before projecting onto the quark-antiquark system, i.e. where the potentials are four-fermion operators (see, for instance, Eq. (42) in Ref. [4]).
For an example of such a field redefinition see Eq. (B13) of Ref. [57]. At the order we are working, this can be approximated by applying the full heavy quark EOMs x 2 ) + · · · , ∂ t ψ † (t, x 1 ) = −i ∇ 2 2m 1 ψ † (t, x 1 ) + i d d x 2 ψ † (t, x 1 )V (0) (|x 1 − x 2 |)χ † c χ c (t, x 2 ) + · · · , (E.4) (and analogously for the antiquark) in the matrix elements. In addition to the free EOMs, they include the static potential. Higher order terms in the coupling constant g and/or in 1/m produce subleading effects. Therefore, we neglect them in Eq. (E.4) for the following discussion.
Eq. (E.4) mixes different orders in 1/m (and sectors with different number of heavy quarks), but still maintains the strict 1/m expansion in the off-shell scheme. As we do not compute the 1/m potential explicitly in this work, instead of using Eq. (E.4), we can make the following replacement in the potentials (p 1 = p, p 1 = p ): and neglect the O(αm 0 ) contributions. In other words, for the matching of the O(α 2 /m 2 ) potentials, we can simply use the free EOMs in Σ (2) , while keeping p 2 − p 2 = 0, as neglecting the O(αm 0 ) contributions sets to zero terms that would contribute to the 1/m potential, which we extract by other means, see Sec. 4. On the other hand, replacing the heavy quark energies by means of the EOMs introduces a potential ambiguity in the determination of the 1/m 2 potentials, since each energy E i , E i can be written in terms of the others by the equality E 1 + E 2 = E 1 + E 2 (energy conservation). This can lead to different results for the individual 1/m 2 potentials. Consider e.g. a term proportional to (E.6) The first line in Eq. (E.6) is zero due to energy conservation. Therefore, we are free to add it to Σ (2) . Transforming the energies E i using Eq. (E.5) generates nonzero contributions to the off-shell terms in V (2,0) , V (0,2) and V (1,1) , as indicated by the second line in Eq. (E.6). However, using Eq. (E.4) in Eq. (E.6) also generates additional contributions to the 1/m potentials such that physical observables, like the heavy quarkonium mass, remain unaffected by the apparent ambiguities. For the k 2 0 -dependent terms we use the prescription As shown in Appendix B of Ref. [18], this transformation does not affect the 1/m potential, because it is based on the continuity equation, which does not contain potential terms. Eqs. (E.5) and (E.7) are relevant for both the FG and the CG result. In addition, we have chosen the prescriptions n f T F k 0 p 2 − p 2 (m 1 − m 2 ) → −n f T F p 2 − p 2 2 , C A k 0 p 2 − p 2 (m 1 − m 2 ) → C A 3m 2 1 + 2m 1 m 2 + 3m 2 2 4m 1 m 2 p 2 − p 2 2 , (E.8) for the energy dependent terms in Σ (2) FG in order to obtain the concrete off-shell matching results in Eqs. (3.15)- (3.16). These prescriptions are motivated by simplicity, and the fact that the resulting off-shell potentials are finite and, therefore, do not require renormalization, see Sec. 5.
Finally, note that the on-shell 1/m 2 potentials resulting from the calculation above are gauge-invariant and independent of the conventions for the field redefinitions we performed.

F Feynman rules for the matching with Wilson loops
In this appendix, we present a set of Feynman rules that can be used to calculate the soft contribution to Wilson-loop expectation values, such as V (1,0) and V L/p 2 ,W in Eqs. Besides the standard set of (static) Feynman rules of NRQCD at leading order in 1/m, the only additional Feynman rules needed to calculate V (1,0) and V L/p 2 ,W are the ones for the chromo-electric field operator E i (t) insertions given in Fig. 5. Because of the explicit factors of t in Eqs. (3.25) and (3.27) we are forced to retain the t dependence in the Feynman rules for E i (t). As a consequence there is no energy conservation at these vertices: Unlike for the pure momentum space Feynman rules for the static potential, the Feynman rules in Fig. 5 do not include an implicit energy conserving (Dirac) δ-function. Three-momentum conservation is however understood as usual. + 1 m 1 − m 2 (5m 2 − 2m 1 ) ln m 1 2m r − (5m 1 − 2m 2 ) ln m 2 2m r , (I.14) where L H = ln n C F α + S 1 (n + l). The following expressions are needed in Eq.