The QCD Potential at $O(1/m)$

Within an effective field theory framework, we obtain an expression for the next-to-leading term in the $1/m$ expansion of the singlet $Q{\bar Q}$ QCD potential in terms of Wilson loops, which holds beyond perturbation theory. The ambiguities in the definition of the QCD potential beyond leading order in $1/m$ are discussed and a specific expression for the $1/m$ potential is given. We explicitly evaluate this expression at one loop and compare the outcome with the existing perturbative results. On general grounds we show that for quenched QED and fully Abelian-like models this expression exactly vanishes.


Introduction
After the discovery of the first heavy-quark bound states, the ψ and the Υ systems, it was soon realized that a non-relativistic picture seemed to hold for them. This is characterized by, at least, three scales: hard (the mass m, of the heavy quarks), soft (the relative momentum of the heavy-quark-antiquark |p| ∼ mv, v ≪ 1), and ultrasoft (US, the typical binding energy E ∼ mv 2 of the bound state system). It was also seen that, if one wanted to describe the whole spectrum of the ψ and the Υ systems, a perturbative evaluation of the potential was not sufficient. This triggered the investigation of these systems by all sorts of potential models (see [1] for some reviews), which are, in general, quite successful phenomenologically. Since then, a lot of effort has been devoted to obtaining the relevant potentials to be used in the Schrödinger equation of such models from QCD [2,3,4,5,6,7,8,9] by relating these potentials with some Wilson loops that could eventually be computed on the lattice or by using some vacuum models. The expression for the leading spin-independent potential, of O(1/m 0 ), has been known since long and corresponds to the static Wilson loop [2,3]. The expressions for the leading spin-dependent potentials in the 1/m expansion, of O(1/m 2 ), have been calculated in Refs. [4,5]. The 1/m corrections to these potentials have proved to be very difficult to obtain. To our knowledge the spin-independent case has been addressed only in Refs. [6,9]. The result of Ref. [6] does not reproduce the one-loop perturbative QCD potential (see discussion at the end of section 5.1) and, therefore, appears to be incomplete. In Ref. [9] the author does not succeed in obtaining suitable finite expressions. We conclude, therefore, that the question of the 1/m corrections to the QCD potential has not been settled yet and hence deserves further studies. In this work we will present an ab initio and systematic calculation of the QCD potential up to O(1/m). We will get a new expression of the 1/m potential that is finite, consistent with one-loop perturbative QCD and suitable to be evaluated by lattice simulations.
We will perform the calculation by integrating out in two steps the hard and the soft scales characterizing the heavy-quark-antiquark system. This is implemented by introducing suitable effective field theories. This approach allows us to express the heavy-quark-antiquark dynamics in terms of systematic and controlled expansions. It has proved to be a powerful computational tool in several different situations. For instance, the hard log corrections (∼ ln m) to the Eichten-Feinberg-Gromes potentials were computed in this way in Ref. [10] (see also [8]). Moreover, we believe that the effective field theories provide a suitable framework where eventually some long-standing conceptual questions will be clarified. In particular, the extent of validity of the naïve potential picture for the heavy quarkonium dynamics, assumed in potential models, could be affected by the consideration of extra degrees of freedom such as hybrids and pions.
The two QCD effective field theories that arise from integrating out the scales m and mv are called NRQCD and pNRQCD respectively. Non-relativistic QCD (NRQCD) was first introduced in Ref. [11]. It has an ultraviolet cut-off much smaller than the mass m and much larger than any other scale (in particular much larger than Λ QCD , which means that the matching from QCD to NRQCD can always be done perturbatively [12,13]). NRQCD has proved to be extremely successful in studyingQQ systems near threshold. The Lagrangian of NRQCD is organized in powers of 1/m, making in this way explicit the non-relativistic nature of the described systems. The maximum size of each term can be estimated by assigning the soft scale to any dimensionful object. In order to connect NRQCD with a potential picture the degrees of freedom with energies much larger than mv 2 have to be integrated out. Once this is done, one is left with a new QCD effective field theory called potential NRQCD (pNRQCD) [14,15] 1 . Strictly speaking, pNRQCD has two ultraviolet cut-offs, Λ 1 and Λ 2 . The former fulfils the relation mv 2 ≪ Λ 1 ≪ mv and is the cut-off of the energy of the quarks, and of the energy and the momentum of the gluons, whereas the latter fulfils mv ≪ Λ 2 ≪ m and is the cut-off of the relative momentum of the quark-antiquark system, p. This theory has been thoroughly studied, and its matching with NRQCD performed, in the situation mv ≫ Λ QCD . In this case the matching can be performed perturbatively. In this paper we will allow Λ QCD to be as large as mv and, therefore, we cannot rely on perturbation theory. Nevertheless, we will assume that the matching between NRQCD and pNRQCD can be performed order by order in the 1/m expansion. We will present, for the general situation Λ QCD < ∼ mv, the matching of NRQCD to pNRQCD at the next-to-leading order in the 1/m expansion in the singlet sector (to be defined later). This will prove to be equivalent to computing the 1/m potential. Formulas that are similar to those in the classical papers of Refs. [4,5,6] will be worked out. No further degrees of freedom with US energy besides the singlet will be considered. This means that non-potential effects will be neglected in this paper (in the perturbative situation this would be equivalent to working at zero order in the multipole expansion). A detailed study of the matching between NRQCD and pNRQCD in the situation Λ QCD < ∼ mv, including US effects, will be worked out elsewhere. The paper is organized in the following way. In section 2 we introduce NRQCD up to order 1/m and discuss its static limit. Moreover, we define what will be pNRQCD in the present context. In sections 3 and 4 we derive the 1/m corrections to the potential, by matching NRQCD to pNRQCD. The Green functions are worked out in the Wilson loop language. The 1/m potential can be written in a simple way as insertions of chromoelectric fields on a static Wilson loop. In section 5 we compute the potential perturbatively up to one loop. Quenched QED and Gaussian models of the QCD long-range dynamics are also discussed. Finally, section 6 is devoted to the conclusions and the Appendix to show how unitary transformations affect the form of the potential.

NRQCD and pNRQCD
The Lagrangian of NRQCD up to order 1/m reads where ψ is the Pauli spinor field that annihilates the fermion, χ is the Pauli spinor field that creates the antifermion, iD 0 = i∂ 0 − gA 0 , and iD = i∇ + gA. The matching coefficients c (j) are not going to be relevant here. For simplicity, light fermions are not explicitly displayed. Their inclusion does not change the results of this paper, although it changes the expression of some intermediate formulas.
The Hamiltonian associated to the Lagrangian (1) is and the physical states are constrained to satisfy the Gauss law. We are interested in the one-quark-one-antiquark sector of the Fock space. It will prove convenient to study the static limit. The one-quark-one-antiquark sector of the Fock space can be spanned by where |n; x 1 , x 2 (0) is a gauge-invariant eigenstate (up to a phase), as a consequence of the Gauss law, of H (0) with energy E (0) n (x 1 , x 2 ); |n; x 1 , x 2 (0) encodes the gluonic content of the state, namely it is annihilated by χ † (x) and ψ(x) (∀x). It transforms as a 3 x 1 ⊗ 3 * x 2 under colour SU (3). The normalizations are taken as follows Notice that since H (0) does not contain any fermion field, |n; x 1 , x 2 (0) itself is also an eigenstate of H (0) with energy E (0) n (x 1 , x 2 ). We have made it explicit that the positions x 1 and x 2 of the quark and antiquark respectively are good quantum numbers for the static solution |n; x 1 , x 2 (0) , whereas n generically denotes the remaining quantum numbers, which are classified by the irreducible representations of the symmetry group D ∞h (substituting the parity generator by CP). The ground-state energy E (0) 0 (x 1 , x 2 ) is usually associated to the static potential (see [15] for important qualifications on this association) and the remaining energies E (0) n (x 1 , x 2 ), n = 0, are usually called gluonic excitations between static quarks. They can be computed on the lattice (see for instance [18]). Translational invariance implies that The gap between different states at fixed r will depend on the dimensionless parameter Λ QCD r. In a general situation, there will be a set of states {n us } such that E (0) nus (r) ∼ mv 2 for the typical r of the actual physical system. We denote these states as ultrasoft. The aim of pNRQCD is to describe the behaviour of the ultrasoft states (pNRQCD has been introduced in [14] and discussed in different kinematic situations in [15]). Therefore, all the physical degrees of freedom with energies larger than mv 2 will be integrated out from NRQCD in order to obtain pNRQCD. In the perturbative situation Λ QCD r ≪ 1, which has been studied in detail in [15], {n us } corresponds to a heavy-quark-antiquark state, in either a singlet or an octet configuration, plus gluons and light fermions, all of them with energies of O(mv 2 ). In a non-perturbative situation, which we will generically denote by Λ QCD r ∼ 1, it is not so clear what {n us } is. One can think of different possibilities. Each of them will give, in principle, different predictions and, therefore, it should be possible to experimentally discriminate between them. In particular, one could consider the situation where, because of a mass gap in QCD, the energy splitting between the ground state and the first gluonic excitation is larger than mv 2 , and, because of chiral symmetry breaking of QCD, Goldstone bosons (pions/kaons) appear. Hence, in this situation, {n us } would be the ultrasoft excitations about the static ground state (i.e. the solutions of the corresponding Schrödinger equation), which will be named the singlet, plus the Goldstone bosons. If one switches off the light fermions (pure gluodynamics), only the singlet survives and pNRQCD becomes totally equivalent to a quantum-mechanical Hamiltonian, thus providing us with a qualitative explanation of how potential models emerge from QCD. In addition, we shall see below how quantitative formulae can be provided in order to calculate the potentials in QCD.
In this paper, we will only study the pure singlet sector with no reference to further ultrasoft degrees of freedom. In this situation, pNRQCD only describes the ultrasoft excitations about the static ground state of NRQCD. In terms of static NRQCD eigenstates, this means that only |0; x 1 , x 2 (0) is kept as an explicit degree of freedom whereas |n; with n = 0 are integrated out 2 . This provides the only dynamical degree of freedom of the theory. It is described by means of a bilinear colour singlet field, S(x 1 , x 2 , t), which has the same quantum numbers and transformation properties under symmetries as the static ground state of NRQCD in the one-quark-one-antiquark sector. In the above situation, the Lagrangian of pNRQCD reads where h s is the Hamiltonian of the singlet (in fact, h s is only a function of r, p 1 , p 2 , which is analytic in the two last operators but typically contains non-analyticities in r), p 1 = −i∇ x 1 and p 2 = −i∇ x 2 . It has the following expansion up to order 1/m In this work we will present the matching between NRQCD and pNRQCD within an expansion in 1/m using the static limit solution as a starting point. Whereas this can be justified within a perturbative framework, in the non-perturbative case the validity of the 1/m expansion cannot be generally guaranteed 3 . We note that this assumption is implicit in all the attempts at deriving the non-perturbative potentials from QCD we are aware of. Furthermore, we would like to emphasize the following important point. The matching calculation can be done independently of what the specific counting in pNRQCD is. As argued in [13,17,15], when doing a matching calculation we are integrating out high energy degrees of freedom. Hence, the Wilson coefficients (potentials) that this integration produces are independent of what the low energy dynamics is. The 1/m expansion just provides a convenient way to organize the matching calculation. One should not conclude that the relative size of the different potentials in the pNRQCD Lagrangian is trivially dictated by the 1/m expansion.

The Wilson loop matching
In this section we will work out the matching up to order 1/m, in a language that is close to the traditional approach of Eichten and Feinberg [4,9]. The matching between NRQCD and pNRQCD is done by enforcing suitable Green functions to be equal in the two theories at the desired order in 1/m. In this section we shall consider space-time Green functions, which are more conventional in non-perturbative studies.
Let us consider an interpolating state with a non-vanishing overlap with the ground state, where φ may in principle be everything that makes the above state overlap with the ground state |0 (0) . We will use here the following popular choice 4 : where P is a path-ordering operator. We also define φ(y, x; t = 0) ≡ φ(y, x). It is going to be useful to introduce a n (x 1 , x 2 ), defined by The states |n (0) being defined up to a phase, it is convenient to fix them in such a way that all the coefficients a n are real. The identification of the singlet from the state (7) depends on the interpolating field used in NRQCD. This dependence is reflected in different normalization factors Z. The matching condition reads As in h s , the normalization factor Z only depends on r, p 1 and p 2 , and is given in the form of an expansion in 1/m as where we have made use of the fact that the NRQCD Lagrangian (1) as well as the pNRQCD Lagrangian (5) are invariant under a CP plus m 1 ↔ m 2 transformation. We will match the Green function G NRQCD defined by with the corresponding Green function in pNRQCD (here and in the rest of the paper, if not explicitly stated, dependence on x 1,2 , p 1,2 is understood)

Calculation in NRQCD
We expand G NRQCD order by order in 1/m Integrating out the fermion fields one gets Analogous formulas hold here and in the following for G (0,1) NRQCD . For simplicity we will not display them. The angular bracket . . . stands for the average value over the Yang-Mills action, W ✷ is the rectangular static Wilson loop the angular bracket . . . ✷ ≡ . . . W ✷ stands for the average over the gauge fields in the presence of the static Wilson loop ( 1 ✷ = W ✷ ). For further convenience we also define . . . ✷ ≡ . . . W ✷ / W ✷ . We use the convention that, if not specified, fields act on the first quark line, e.g. t), and so on. We have used time reversal: B(t) ✷ = − B(−t) ✷ to eliminate the spin-dependent term in Eq. (14). Therefore, we can already state that no spin-dependent potential appears at O(1/m) [4].
It is useful to introduce at this point two identities involving covariant derivatives and Schwinger lines: In the first Ref. [4], both identities i) and ii) were derived. Identity ii) corrects their Eq. (4.7c). As a by-product of the above identities, we get where Let us note that time-reversal symmetry gives Using the above relations, Eq. (14) can be worked out giving

Calculation in pNRQCD
Let us consider the Green function defined in Eq. (12). Expanding it up to order 1/m, we obtain Inserting Eqs. (6) and (10) into Eq. (12), we obtain where ∇ = ∇ r . As in the NRQCD case, we do not display the analogous formulas for G (0,1) pNRQCD . In order to keep Eqs. (17) and (18) simpler, we have already used the fact that Z 0 and Z 1,p can be chosen as real. This follows from the matching to G (1,0) NRQCD (16) once any constant phase in Z 0 is conventionally set to zero.

Matching
At O(1/m 0 ) we match Eq. (13) with Eq. (17). We obtain At O(1/m) we match Eq. (16) with Eq. (18). We obtain A similar expression with neither end-point string contributions nor normalization factors has been obtained in [9]. The right-hand side of Eqs. (19) and (20) can be shown to be finite in the T → ∞ limit. This is not obvious from the point of view of a pure Wilson-loop calculation in the spirit of Ref. [4], as is apparent from the difficulties met by the author of [9] (the normalization factors are crucial in order to get a finite expression). For this reason, and because such kind of arguments may become relevant to future analyses of Wilson-loop correlators, we mention here the relevant steps of the proof. a) Inserting the identity operator |n (0) (0) n| into the Wilson-loop average, the latter may be written as Doing the same for the average E(t) ✷ one obtains d) Inserting the identity operator into the correlator E(t) · E(t ′ ) ✷ , it may be written as for t > t ′ . With the above points a)-d) it is easy to show that the right-hand side of Eqs. (19) and (20) is finite in the large-T limit. Their explicit expressions read Analogously, we can obtain matching expressions for the normalization factors. At O(1/m 0 ), the result is Z 0 = |a 0 | 2 . At O(1/m) a more complicated equality is obtained, which involves a combination of Z 1 and Z 1,p . The fact that we have two equations and three independent functions: Z 1 , Z 1,p and V 1 at O(1/m) makes the determination of V 1 ambiguous. This ambiguity is intrinsic in the sense that any value of Z 1,p and V 1 that fulfils Eq. (22) will lead to the same physics (as far as one consistently works at higher orders in 1/m), and has to do with the fact that a quantum-mechanical Hamiltonian is defined up to time-independent unitary transformations (see Appendix). As a consequence, the ambiguity in the definition of V 1 also gives rise to ambiguities in the definition of the potentials of order O(1/m 2 ) and higher. In the next section we will fix Z 1,p (and then V 1 ) by imposing an extra matching condition, which will prove to be particularly convenient.

The quantum-mechanical matching
In the previous section we have done the matching comparing Green functions in NRQCD and in pNRQCD. In this section the comparison is between states and matrix elements in NRQCD and in pNRQCD. The calculation, in terms of states, will be closer in philosophy to the usual quantum mechanics calculations in perturbation theory [20] (see also [21]). Moreover, the whole procedure will share some similarities to the adiabatic approximation and the Born-Oppenheimer method as used in atomic physics calculations [22]. The underlying assumption is that the difference of energies among states labelled with different n is much larger than the difference of energies among states labelled with the same n. In our case, since we only aim at correctly reproducing the ground state spectrum, we only need that the splitting between the ground state and the first gluonic excitation be larger than the typical splitting of the states of the ground state (taken of O(mv 2 ) by definition). This is nothing but the condition we have assumed throughout the present paper.
H is not diagonal in the basis of H (0) (|n (0) ) with respect to the n labelling. We consider instead a basis of states, labelled as |n; such that the Hamiltonian H is diagonal with respect to these states where E n (x 1 , x 2 , p 1 , p 2 ) is an analytic function in p 1 , p 2 , and such that they are normalized as m; x 1 , x 2 |n; y 1 , Conditions (24) and (25) give A set of states |n and an operator E n that satisfy Eqs. (25) and (26) can be obtained from the static solutions |n (0) and E (0) n to any desired order of accuracy by working out formulas analogous to the ones used in standard quantum mechanics perturbation theory [20]. If we write |n and E n as an expansion in 1/m, we obtain at O(1/m) for n = 0 (when not specified, states and energies are calculated in x 1 and x 2 ): These equations may be derived from the identities which follows from symmetry considerations and which follows from explicit calculation. Analogous formulas hold for the antiparticle contribution.

Matching
The aim of pNRQCD is to describe the behaviour of |0 . The integration of high excitations is trivial using the basis (23) since, in this case, they are decoupled from |0 . Therefore, the matching of NRQCD to pNRQCD is basically to rename things in a way such that pNRQCD reproduces the matrix elements of NRQCD for the ground state. The matching conditions in this formalism read as follows: Using Eq. (31) and Eq. (9) we obtain the normalization factor, given the interpolating field, It is explicit now that h s and, hence, the potential, fixed by Eq. (31), depend neither on the normalization factors nor on the specific shape of the end-point strings. Moreover, Eqs.
(31) and (32) provide matching conditions at any finite order in 1/m. From Eqs. (31) and (32) one could draw the conclusion that the ambiguity in the computation of V 1 , discussed in the previous section, has disappeared. This is not the case. Actually, Eqs. (27) and (28) with Eqs. (29) and (30) only give one of the possible solutions of Eqs. (26) and (25). Indeed, these equations do not completely fix the state |n . In standard quantum mechanics the state is fixed up to an arbitrary constant phase. In our case, since we only diagonalize in the n space, this phase becomes a unitary operator 5 ; given a state |n that fulfils Eqs. (26) and (25), this means that the states with O † n = O n , still fulfil Eq. (25) and the Hamiltonian is still diagonal in n with This freedom reflects the fact that h s is defined by the matching up to an arbitrary unitary transformation or, which is the same, up to an arbitrary unitary field redefinition: and so does the Z. In order to make calculations easier, we have taken advantage of this freedom by fixing the relative phase between |0 and |0 (0) following the standard choice of quantum mechanics perturbation theory. With this choice, we obtain Eqs. (29) and (30); as we will see below, this will allow us to obtain a compact expression for V 1 in terms of Wilson loops. Furthermore, the procedure can be easily generalized at any finite order in 1/m. From Eqs. (29) and (32) we can perform the matching at O(1/m). We obtain From Eqs. (30) and (31) we get, up to O(1/m): At first sight, this result seems of limited practical utility, since it depends on the exact and complete solution of the bound state at order 1/m 0 . Fortunately this is not the case. In fact, Eq. (20) suggests a very simple form for the 1/m potential, Using the same technique as outlined in a)-d) of Sec. 3.3, it can be proved that Eq. (38) is finite and equal to Eq. (37). Equation (38) is the main result of the present work. It gives the leading 1/m correction to the heavy-quark potential expressed in terms of chromoelectric field insertions in a static Wilson loop, thus avoiding the explicit computation of the normalization factors. The potential appears, therefore, in the same form as the potential calculated in Refs. [4,5,6] and is suitable for lattice evaluations similar to those performed in [8]. Moreover, the above expression expanded in α s gives the full perturbative series of the 1/m potential in the regime where this expansion makes sense (i.e. mv ≫ Λ QCD ). In the next section we will calculate from it the leading non-vanishing perturbative contribution to the 1/m potential, which will coincide with the results of Refs. [23,24]. For the calculation of some higher order perturbative corrections to the V 1 potential we refer to [25]. Finally, it is interesting to note that Eq. (38) also holds if the correlators are evaluated on a Wilson loop with infinitely large time extension (i.e. if we make the limit T → ∞ in the Wilson loops for t and t ′ finite, and then evaluate the integrals). This observation greatly simplifies the perturbative calculation and may also do so for non-perturbative ones.

Applications
In this section we will consider Eq. (38) in perturbative QCD and in quenched QED. In the first case we will evaluate the contribution to the potential at the leading non-vanishing order in α s . In the second we will show that the potential V 1 vanishes exactly. This may be relevant to several Gaussian models used in the phenomenology of the non-perturbative dynamics of the strong interaction. We will shortly comment on this.

Perturbative QCD
We perform the calculation in the Coulomb gauge. Since Eq. (38) is a gauge-independent quantity, the result will hold in any gauge. At order α s , the right-hand side of Eq. (38) only gives a self-energy type of contribution, since the chromoelectric fields are inserted on the same quark line. The first non-vanishing contribution to the potential (i.e. depending on r) appears at order α 2 s . Three types of diagrams arise: i) Disconnected graphs. An example is shown by graph a) of Fig. 1. These graphs cancel in the difference of the right-hand side of Eq. (38). ii) Connected graphs with end-point strings contributions. Non-disconnected graphs involving gluons attached to the end-point strings vanish in Eq. (38). As noticed in the previous section the T → ∞ limit can also be performed in the correlators before doing the time integrals in Eq. (38). This turns out to be quite useful here, making the irrelevance of these graphs to the potential manifest. iii) Connected graphs with no end-point strings contributions. The only non-vanishing graph contributing to Eq. (38) is graph b) of Fig. 1. The analogous graph involving the triplegluon vertex, but with two transverse gluons attached to the chromoelectric fields on the quark line, can be shown to vanish in the limit of Eq. (38) by explicit calculation. Therefore, considering all relevant contributions (which reduce to graph b) of Fig. 1) we get from Eq. (38) at order α 2 Equation (39) coincides with the result of Refs. [23,24] (which may also be obtained using the rules of [17]) and with the non-Abelian part of the perturbative calculation of the 1/m potential done in [26]. In Ref. [26] also an Abelian contribution to the 1/m potential is found. We have stressed throughout the paper that there is no unique way to define the 1/m potential and, therefore, different matching procedures are, in general, expected to give different expressions for it. One may suspect that this is due to the fact that in Refs. [23,24] the Coulomb gauge is used whereas in [26] the Feynman gauge is used. However, this has more to do with the freedom we have in quantum mechanics to change the form of the Hamiltonian by a unitary transformation without changing the physics, rather than with the gauge fixing dependence itself. What typically happens in perturbative calculations, which match non gauge invariant Green functions, is that depending on the gauge one uses, one gets a different form of the potential. The different forms are equivalent and can be obtained from each other by unitary transformations. More precisely, from the discussion in the Appendix it follows that V pert 1 /m, at O(α 2 s ), can be rewritten in terms of an O(α s /m 2 ) potential 6 . However, whereas the Abelian piece found in [26] can be obtained from the general expressions for the O(1/m 2 ) potential [6], this is not so for the non-Abelian piece in Eq. (39), since there is not non-Abelian contribution at tree level. Therefore, we conclude that Eq. (38) represents a genuine new potential not considered in the past in the nonperturbative parametrizations of the QCD potential in terms of Wilson loops. Finally, let us mention that, to our knowledge, V pert 1 has never been computed using Wilson loops before.

Quenched QED
Let us consider quenched QED. In this situation the Wilson-loop average is exactly known (it reduces to Gaussian integrals, see for instance [27]) and can be written as where f ij (t, r) = f ji (t, r) = E i (t, r)E j (0, 0) . As a direct consequence (for a derivation, see for instance [28]) we obtain Therefore, since the term f ii (t − t ′ , 0) is a self-energy-type contribution, the potential contribution of V 1 (as defined in Eq. (38)) vanishes exactly in quenched QED. The same result is obtained in Ref. [17]. Several QCD vacuum models [1,7,29] seem to approximate the Wilson-loop long-range non-perturbative dynamics with expressions analogous to Eq. (40). The above considerations made for quenched QED may be relevant to them. As a matter of fact, 1/m corrections do not seem to show up there. Therefore, it is tempting to generalize the exact result of quenched QED to Gaussian models of the QCD vacuum. Nevertheless some words of caution are needed. A model consists in an approximation of QCD, which is supposed to coincide with a relevant limit on the QCD dynamics. This limit is unknown by definition. Therefore, how to implement it is, in the practice of several models, not well established. While it is clear that starting from a Gaussian expression for the Wilson loop a relation like Eq. (41) is going to hold, it is not guaranteed that 1/m potentials will not show up if the Gaussian approximation of the Wilson loop is implemented at some intermediate step of the calculation 7 .

Conclusions
We have obtained an exact and non-perturbative expression for V 1 , the QQ potential at order 1/m. This expression is shown in Eq. (38). The perturbative calculation of it up to order α 2 s agrees with [23,24] and is consistent with the one-loop contribution of [26]. However, the existence of a non-perturbative contribution at order 1/m was, to our knowledge, not considered before in the literature, be it in phenomenological applications [1] or in attempts at obtaining the non-perturbative potentials from QCD [4,5,6]. We note that via a unitary transformation the 1/m terms can, in certain circumstances, be reshuffled in 1/m 2 (and higher) terms in the potential (see Appendix). When these circumstances apply (i.e. when V 1 ≪ m 2 v 2 ), our results do not have immediate consequences on phenomenological models where the full set of 1/m 2 potentials are considered and input as ansätze functions. However, when they do not (i.e. when V 1 ∼ m 2 v 2 ), our results imply that a 1/m potential should be included in those models. In either case, they are extremely important for the attempts at obtaining the potentials from lattice QCD, since, upon doing the above-mentioned reshuffling, the 1/m potential given by Eq. (38) generates 1/m 2 terms different from the ones calculated so far [6].
The most promising way to have an estimate of the non-perturbative behaviour of the correlators appearing in Eq. (38) is by a lattice simulation. In practice this can already be done from the lattice data of Ref. [8], since the correlators we get are of the same type. Such lattice data could be of interest for at least two reasons. If the potential found on the lattice happens to be small in the long range, it would support the Abelian dominance picture. In the short range it should show the interplay of the perturbative (known) region with the non-perturbative one.
When there are no additional US degrees of freedom, the adopted non-perturbative procedure can be easily generalized to the evaluation of higher-order terms in the inverse mass expansion [30]. The inclusion of US degrees of freedom in the non-perturbative regime will be discussed elsewhere.

Appendix A
In this appendix we show that there exists a unitary transformation that reshuffles momentum-independent 1/m terms in 1/m 2 momentum-dependent and momentum-independent terms. Let us consider the Hamiltonian By choosing all the 1/m potentials in the Hamiltonian (43) disappear. The price to pay is the appearance of new terms at order 1/m 2 (and higher). Of course, the leading size of these new terms is the same as the original V 1 /m. In particular, since p, ∇, ∼ mv and V 0 ∼ mv 2 , we get W ∼ V 1 /(m 2 v 3 ) and {p i , {p j , (∇ i W j )}}/m 2 ∼ V 1 /m. The size of the remaining 1/m 2 induced terms is V 2 1 /(m 3 v 2 ) and hence, it depends on the size of V 1 /m, which is a priori unknown. On general grounds the maximum size of V 1 is given by the largest available scale, namely mv, and hence at most V 1 ∼ m 2 v 2 . Reasoning in the same way, the 1/m 2 potentials (calculated via the quantum-mechanical matching in [30]) are not bigger than mv 3 . However, from the condition {W, p} ≪ m, it follows that the reshuffling of the 1/m potential to 1/m 2 potentials may be done in the way above only if V 1 ≪ m 2 v 2 . As a consequence, all the terms of O(V 2 1 /(m 3 v 2 )) have a size much smaller than mv 2 . More specifically, if, for instance, V 1 ∼ m 2 v 3 , then the terms of O(V 2 1 /(m 3 v 2 )) are of order mv 4 and hence suppressed by a factor v with respect to the original 1/m potential (as well as with respect to the 1/m 2 potentials obtained from the quantum-mechanical matching). We note that in perturbation theory (v ∼ α s ) V 1 ∼ m 2 v 4 and, therefore the terms of O(V 2 1 /(m 3 v 2 )) are suppressed by a factor v 2 .