Renormalization group improvement of the NRQCD Lagrangian and heavy quarkonium spectrum

We complete the leading-log renormalization group scaling of the NRQCD Lagrangian at $O(1/m^2)$. The next-to-next-to-leading-log renormalization group scaling of the potential NRQCD Lagrangian (as far as the singlet is concerned) is also obtained in the situation $m\alpha_s \gg \Lambda_{QCD}$. As a by-product, we obtain the heavy quarkonium spectrum with the same accuracy in the situation $m\alpha_s^2 \simg \Lambda_{QCD}$. When $\Lambda_{QCD} \ll m\alpha_s^2$, this is equivalent to obtain the whole set of $O(m\alpha_s^{(n+4)} \ln^n \alpha_s)$ terms in the heavy quarkonium spectrum. The implications of our results in the non-perturbative situation $m\alpha_s \sim \Lambda_{QCD}$ are also mentioned.


Introduction
Heavy quark-antiquark systems near threshold are characterized by the small relative velocity v of the heavy quarks in their center of mass frame. This small parameter produces a hierarchy of widely separated scales: m (hard), mv (soft), mv 2 (ultrasoft), ... . The factorization between them is efficiently achieved by using effective field theories, where one can organize the calculation as various perturbative expansions on the ratio of the different scales effectively producing an expansion in v. The terms in these series get multiplied by parametrically large logs: ln v, which can also be understood as the ratio of the different scales appearing in the physical system. Again, effective field theories are very efficient in the resummation of these large logs once a renormalization group (RG) analysis of them has been performed. This will be the aim of this paper for the cases of NRQCD [1] and potential NRQCD (pNRQCD) [2,3] 1 .
First, we will obtain the RG improved matching coefficients of the NRQCD Lagrangian at one loop and up to O(1/m 2 ). Since, by construction, the matching coefficients of HQET are equal to the analogous ones of NRQCD, these can already be obtained from the literature [4,5]. Therefore, only the four-heavy fermion matching coefficients need to be computed to obtain the complete leading-log (LL) RG improvement of the NRQCD Lagrangian at O(1/m 2 ). We will perform such a calculation in this paper. For the spin-dependent fourheavy fermion matching coefficients there already exists a computation in Ref. [6]. We differ with their results. Our evaluation is relevant in the study of the situation mα s ∼ Λ QCD . On the one hand they could be used to improve the accuracy of phenomenological studies or lattice simulations of NRQCD. On the other hand, this situation has been studied within an effective field theory framework in [7] where the matching between NRQCD and a Schrödinger like formulation has been achieved in a controlled fashion. In particular, it is possible to write the potentials as Wilson loops multiplied by the matching coefficients inherited from NRQCD. The obvious application is that the matching coefficients here computed are the ones that multiply the Wilson loops in the non-perturbative potentials. This is especially relevant now that the complete expression for the potential at O(1/m 2 ) in terms of Wilson loops is available [7,8,6]. In particular, it would be welcome to have an updated evaluation of the lattice analysis of the heavy quarkonium spectrum made in Ref. [9] taking into account the complete O(1/m 2 ) potential as well as the now known complete set of LL NRQCD matching coefficients.
In the situation when mα s ≫ Λ QCD , the matching between NRQCD and pNRQCD, i.e. the computation of the potentials, can be done perturbatively. In this case ultrasoft gluons as well as the quark-antiquark in an octet configuration do exist at the matching scale between NRQCD and potential NRQCD producing further divergences. By taking into account these divergences as well as the divergences computed before we have obtained the next-to-next-to-leading-log (NNLL) RG improved pNRQCD Lagrangian (as far as the singlet is concerned).
If we are in the situation Λ 3 QCD /(mα s ) 2 ≪ mα 2 s , the leading solution of the spectrum corresponds to a Coulomb-type bound state and the non-perturbative effects are corrections. In this situation, by using the previous result of the NNLL RG improved pNRQCD Lagrangian, we are able to obtain the heavy quarkonium spectrum with the same accuracy.
If instead, we are in the situation Λ QCD ≪ mα 2 s , from our previous result, we are able to obtain the whole set of O(mα (n+4) s ln n α s ) terms in the heavy quarkonium spectrum. There already exists an evaluation [10] within the vNRQCD framework [11,12,13] of the RG improved Heavy Quarkonium mass when Λ QCD ≪ mα 2 s . We agree for the spin-dependent terms (since we do for the spin-dependent potentials computed in Ref. [12]) but differ for the spin-independent ones.

NRQCD
NRQCD has an ultraviolet (UV) cut-off ν N R = {ν p , ν s } satisfying mv ≪ ν N R ≪ m. At this stage ν p ∼ ν s . ν p is the UV cut-off of the relative three-momentum of the heavy quark and antiquark. ν s is the UV cut-off of the three-momentum of the gluons and light quarks. This does not seem to give problems at the order we are working at but one should be eventually careful upon the possible gauge dependence of this splitting.
Indeed, the above cutoffs plus the matter content of the theory given below correspond to our definition of NRQCD. Within the threshold expansion framework [14] this corresponds to integrate out the hard modes of QCD in order to obtain NRQCD. Unfortunately, NRQCD already contains non-physical degrees of freedom for the phase space region it is aimed to describe, which implies that the terms in the Lagrangian do not have a unique size nor, therefore, power counting (to avoid this problem is one of the motivations for the construction of pNRQCD, which will be done in the next section). Nevertheless, this poses no problem for the NRQCD running considered in this section.
The NRQCD Lagrangian including light fermions reads at O(1/m 2 ) (up to field redefinitions) [1,4,5] L where ψ is the Pauli spinor that annihilates the fermion, χ is the Pauli spinor that creates the anti-fermion, iD 0 = i∂ 0 − gA 0 , iD = i∇ + gA, analogously for L χ and We have also included the D 4 /m 3 term above since it will be necessary in the evaluation of the heavy quarkonium mass once the power counting is established. Moreover, we will consider that the kinetic term matching coefficients are protected by reparameterization invariance (c k = c 4 = 1) [15], however, we will often keep them explicit for tracking purposes.
The NRQCD matching coefficients are functions of ν N R = {ν p , ν s }. Somewhat by definition, the matching coefficients of the bilinear in the heavy quark fields and of the pure gluonic terms are just functions of ν s , i.e. c = c(ν s , m) ≡ c(ν s ). In any case, it will explicitely come out from the calculation. The complete LL running of these matching coefficients in the above basis (2-4) have been calculated by Bauer and Manohar [5] in the (background) Feynman gauge 2 (some partial previous results already exist in the literature [4]). Therefore, in order to complete the RG running of the NRQCD Lagrangian we only need to compute the four-heavy-quark operators with LL accuracy. This will be our aim in the following.
A procedure to get the ν s dependence of the NRQCD matching coefficients is by using HQET-like rules in NRQCD (by this we mean to perform the perturbative expansion in 1/m prior to the computation of the Feynman integrals). In fact, in our case, at one loop, all the dependence of the matching coefficients is only due to ν s since no ν p dependence appears at one-loop, i.e. d(ν p , ν s , m) ≡ d(ν p , ν s ) ≃ d(ν s ). This will be discussed below within pNRQCD.
Formally, we can write the NRQCD Lagrangian as an expansion in 1/m in the following way: where the above fields and parameters should be understood as bare and the renormalization group equations of the renormalized matching coefficients read The RG equations have a triangular structure (the standard structure one can see, for instance, in HQET RG evaluations, ie. for the Lagrangian (1) setting the heavy antiquark field to zero): where the different B's can be power-expanded in λ 0 (λ 0 corresponds to the marginal operators (renormalizable interactions)). For NRQCD we have λ 0 = α s and λ 1 = {c k , c F }, At this stage, we would like to stress that we are working in a non-minimal basis of operators for the NRQCD Lagrangian. Consequently, the values of (some of) the matching coefficients are ambiguous (only some combinations with physical meaning are unambiguous). In particular, some of the matching coefficients could depend upon the gauge in which the calculation has been performed. Therefore, it is important to perform the matching calculation in the same gauge (at least for those operators which could suffer the ambiguity). We will further discuss this point latter on.
The RG equations for the {c} in the Feynman gauge can be read from Bauer and Manohar results [5]. Because of latter comparison, we only explicitely write the equation for where, These equations have been obtained by explicit computation in the Feynman gauge by Signer [16] within the threshold formalism [14]. We have obtained them by using the results of Ref. [17], which were performed in the Feynman gauge, plus doing the explicit calculation of the terms that depend on c D in NRQCD in the Feynman gauge. This proves to be enough since the dependence on c 2 k of Eq. (11) can be inferred from the results of Ref. [17] once the dependence on c D is known (since the spin-dependent terms will depend on c 2 F ). Both 3 The renormalization group evolution of the one-heavy quark sector has also been done in a minimal basis in Ref. [5] by elliminating the operator multiplying c hl 1 . In that case the expression of c D obtained in Ref. [5] is indeed gauge independent. 4 For the record, we also display the non-equal mass case equations with the definitions d/m 2 → d/(m 1 m 2 ).
The equations for d ss and d sv remain equal, for d vv one has to change c 2 F and the equation for d vs reads calculations agree. Note that it was needed not to have ν p dependence at one-loop in order the argument to go through.
As we have mentioned we are not working in a minimal basis. This shows up in the ambiguity of the value of the matching coefficients of some operators. At the practical level, this means that they will depend on the specific basis of operators we have taken for the NRQCD Lagrangian and on the procedure used (in particular on the gauge). Therefore, if working in a non-minimal basis, one should be careful and do the matching using the same gauge for all the operators (or at least for those that are potentially ambiguous).
For illustration, let us consider the case without light fermions. In this case, c D and d vs are ambiguous but not an specific combination (see Eq. (34)). In particular, c D could be absorbed by other matching coefficients by using a field redefinition [5]. We can check these statements by doing the calculation in the Coulomb gauge. In this case we obtain (no non-trivial change is now required in the non-equal mass case for the RG equations of the four-heavy fermion matching coefficients) One can see that, as far as the combination that appears in Eq. (34) is concerned, the physical result is unchanged. In the following we will use the Feynman gauge results for the NRQCD matching coefficients.
With the above results we have completed the RG equations of the NRQCD Lagrangian at one loop at O(1/m 2 ). In order to solve these equations, we need the (tree-level) matching conditions of the matching coefficients at some matching scale. We choose as the matching scale m. The {c(m)} can be read, for instance, from [5]. The tree-level matching conditions for the four-heavy fermion operators read We can then obtain the solution of the RG equations. We only explicitely display those which are new or will be necessary later on (we define z = αs(νs) The {c} matching coefficients can be found or deduced from the results in Ref. [5]. Finally note that it is very important to know the basis of operators one has been working in NRQCD as well as in which gauge the calculation has been performed. In practice this means that one should make sure that c D and d vs have been computed in the same way in order to obtain the correct result.

pNRQCD
The above results are also a necessary step towards the RG improvement of pNRQCD when mα s ≫ Λ QCD , which we consider in the following.
By integrating out some soft degrees in NRQCD one ends up in pNRQCD. This latter theory is defined by the cut-off ν pN R = {ν p , ν us }, where ν p is the cut-off of the relative threemomentum of the heavy quarks and is such that mv ≪ ν p ≪ m and ν us is the cut-off of the three-momentum of the gluons and light quarks with mv 2 ≪ ν us ≪ mv. In principle, we do not rule out the option of correlating ν p with ν us in order to efficiently perform the renormalization group improvement at higher orders [11]. Nevertheless, at the order we are working with, we not need to specify any relation between ν p and ν us since the dependence on ν p would be a subleading effect. Therefore, in this paper, we will treat them as independent.
The pNRQCD Lagrangian reads as follows: where we have explicitly written only the terms relevant to the analysis at the NNLL of the singlet sector; S and O are the singlet and octet field respectively. All the gauge fields in Eq. (15) are functions of the center-of-mass coordinate and the time t only. For a more extensive discussion we refer the reader to Refs. [3,18].

Potentials
We now display the structure of the matching potentials V (0) s , V (0) o , V (1) s and V (2) s , which are the relevant ones to our analysis.
1) Order 1/m 0 . From dimensional analysis, V (0) s (r) can only have the following structure: and similarly for the static octet potential: 2) Order 1/m. From dimensional analysis and time reversal, V (1) s can only have the following structure: 3) Order 1/m 2 . At the accuracy we aim, V (2) s has the structure where S 12 (r) ≡ 3r · σ 1r · σ 2 − σ 1 · σ 2 and S = σ 1 /2 + σ 2 /2. The coefficients, α Vs , D s , ... contain some ln r dependence once higher order corrections to their leading (non-vanishing) values are taken into account. In particular, we will have expressions like δ (3) (r) ln n r. This is not a well-defined distribution and should be understood as the Fourier transform of ln n 1/k. Nevertheless, in order to use the same notation for all the matching coefficients, and since it will be sufficient for the purposes of this paper, resum the leading logs, we will use the expression δ (3) (r) ln n r, although it should always be understood in the sense given above.

RG equations
Formally, we can write the pNRQCD Lagrangian as an expansion in 1/r(= 1/r, p) and 1/m in the following way: where the above fields and parameters should be understood as bare. As for the renormalized quantities, we define V as the potentials andṼ as the (almost) dimensionless constants in it. The latter are in charge of absorbing the divergences of the effective field theory. Therefore, they will depend on ν p and ν us . Note that the dependence on ν s of the NRQCD matching coefficients has to cancel inṼ since the new effective theory does not have any ultraviolet cutoff dependent on ν s . This discussion completely fixes the procedure to obtain the RG equations of the potentials: by studying the UV behavior of pNRQCD it is possible to obtain the scale dependence of the potentials on ν p and ν us and the independence on ν s trivially sets the ν s scale (in-)dependence of the potentials. Being more specific, the potentials have the following structure: V (d(ν p , ν s , m), c(ν s , m), ν s , ν us , r) =Ṽ (ν p , m, ν us , r) ≡Ṽ (ν p , ν us ) .
This produces the following RG equations: The first equation just reflects the independence of the potential on ν s . At the practical level, with the accuracy we are working, it is equivalent to set ν s = 1/r up to factors of order one. The second equation tells us that the dependence on ν p is inherited from the (four-heavy fermion) NRQCD matching coefficients. One of our aims will be to obtain the heavy quarkonium spectrum with NNLL accuracy when Λ 3 QCD /(mα s ) 2 ≪ mα 2 s . In that situation the leading order solution corresponds to a Coulomb-type bound state and, leaving aside non-perturbative corrections, a perturbative expansion is licit. In order to achieve this goal we will need the RG improvement of the pNRQCD Lagrangian for the singlet sector with the same accuracy. Being more precise, we will need ν d dνṼ as well as (due to mixing) within an strict expansion in α s .
Let us consider first the dependence on ν p . Eq. (22) tells us that the dependence on ν p appears due to the four-heavy fermion matching coefficients. These first appear at O(1/m 2 ). Therefore, we only need to obtain the equation ν p d dνpṼ s . In fact, at least at lowest non-vanishing order, only the delta potentials are dependent on ν p . Of this type is precisely the leading dependence on ν p of pNRQCD. It appears through (the iteration of) enough singular potentials when performing standard quantum mechanics perturbation theory. Explicit inspection shows that these kind of effects first appear at O(mα 6 s ) in a perturbative computation of the mass (this argument is based on the singularity of the (iteration of the) potentials plus taking into account their leading non-vanishing power in the α s expansion). Therefore, we obtain By using Eq. (22), this is equivalent to ν p d dνp d = 0 + O(α 3 s ), which was needed previously in NRQCD. In principle, this result could also be proved by explicit inspection on the possible diagrams at the quark-gluon level that could give divergences proportional to ν p . The final conclusion is that we can neglect any dependence on ν p at the order we are working at in the potentials, i.e.Ṽ (ν p , ν us ) ≃Ṽ (ν us ). Therefore, we only have to compute the ν us scale dependence.
The ν us -scale dependence could be taken from the computation in [19,18,20,3] (see also [21]) by keeping track of the dependence of the result on the differentṼ . Let us note that we only need to do a one-loop computation in order to achieve the necessary accuracy (plus the already known two-loop singlet static potential). This should be compared with the one-, two-and three-loop calculations that seem to be necessary if the calculation is performed at the quark-gluon level as in Ref. [12,13,10].
Formally, the renormalization group equations of the renormalized matching coefficients due to the ν us -dependence read From a practical point of view one can organize the RG equations within an expansion in 1/m.
At O(1/m 0 ), the analysis corresponds to the study of the static limit of pNRQCD, which has already been carried out in Ref. [20]. We repeat the basic points here for ease of reference. SinceṼ −1 = 0, there are relevant operators (super-renormalizable terms) in the Lagrangian and the US RG equations lose the triangular structure that we enjoyed for the RG equations of ν s . Still, ifṼ −1 ≪ 1, a perturbative calculation of the renormalization group equations can be achieved as a double expansion inṼ −1 andṼ 0 (for a similar discussion in the context of scalar λφ n -like theories see [22]), where the latter corresponds to the marginal operators (renormalizable interactions)). At short distances (1/r ≫ Λ QCD ), the static limit of pNRQCD lives in this situation. Specifically, we haveṼ −1 = {α Vs , α Vo }, that fulfills V −1 ∼ α s (r) ≪ 1;Ṽ 0 = α s (ν us ) andṼ 1 = {V A , V B } ∼ 1. Therefore, we can calculate the anomalous dimensions order by order in α s (ν us ). In addition, we also have an expansion iñ V −1 . Moreover, the specific form of the pNRQCD Lagrangian severely constrains the RG equations general structure. The result obtained in Ref. [20] reads At higher orders in 1/m, we only need to consider the singlet potentials. The same considerations than for the static limit apply here as far as the non-triangularity of the RG equations is concerned. At O(1/m, 1/m 2 ), we have the following matching coefficients: V (1) 1,s , D 2,s , D d,s , D S 2 ,s , D LS,s , D S 12 ,s }, and we obtain and zero for the other matching coefficients (in particular for the spin-dependent potentials).
In a more formal way, Eq. (28) has the following structure: In fact, in general, we have the structure ( and one has to pick up the leading contributions from all the possible terms. Eqs. (21), (25), (27) and (28) provide the complete set of RG equations at the desired order. By using Eqs. (21) and (25), we obtaiñ V =Ṽ (d(1/r), c(1/r), ν s = 1/r, ν us , r) .
We now need the initial condition in order to solve the US RG equations, i.e. the matching conditions. We fix the initial point at ν us = 1/r. In summary, we need to know the singlet static potential with O(α 3 s ) accuracy, the 1/m potential with O(α 2 s ) accuracy, the 1/m 2 potentials and the singlet octet potential with O(α s ) accuracy and V A with O(1) accuracy at ν us = 1/r. They read 1,s (r −1 ) = α s (r −1 ), D 2,s (r −1 ) = α s (r −1 ), D where β 1 = 34/3C 2 A −4C f T F n f −20/3C A T F n f and the values of a 1 and a 2 have been computed in Ref. [23]. We now have all the necessary ingredients to solve the RG equations.
Eqs. (27) and (28) give rise to subleading effects within an strict expansion in α s . Therefore, we can approximate them to (if not displayed the RG equation remains equal) 1,s = 8 3 We can finally obtain the RG improved potentials for the singlet: 1,s (ν us ) = D 1,s (r −1 ) + 2,s (ν us ) = D 2,s (r −1 ), LS,s (ν us ) = D This completes the RG evaluation of the pNRQCD Lagrangian at NNLL (as far as the singlet is concerned).

Heavy quarkonium spectrum
In the situation Λ 3 QCD /(mα s ) 2 ≪ mα 2 s , the heavy quarkonium behaves as a Coulomb-type bound state and the perturbative corrections can be computed in a systematic way. From the potential-like terms, we obtain the following correction to the NNLO energy expression (the derivation of this result would go along similar lines to those in Ref. [18]): where E n = −mC 2 f α 2 s /(4n 2 ), the scale ν s in z and in the NRQCD matching coefficients has been fixed to the soft scale . α s is also understood at the soft scale ν s = 2a −1 n unless the scale is specified, and Eq. (35) gives all the O(mα 4 s (α s ln) n ) terms for n ≥ 1 of the heavy quarkonium mass, where ln stands either for ln(α s ), arising from the hard scale, or for ln(mα s /ν us ), arising from the ultrasoft scale. After adding to Eq. (35) the NNLO [24] result with the normalization point at the same soft scale, ν s = 2a −1 n , that we have used here, the complete (perturbative) NNLL heavy quarkonium mass is obtained (note that for the LO result the three-loop running α s has to be used).
The ν us dependence of Eq. (35) cancels against contributions from US energies. At the next-to-leading order in the multipole expansion, the contribution from these scales reads where H o = c k p 2 m + V (0) o and ν us is the UV cut-off of pNRQCD. Then, the total correction to the energy reads δE n,l,j = δ pot E n,l,j (ν us ) + δ US E n,l (ν us ) .
Different possibilities appear depending on the relative size of Λ QCD with respect to the US scale mα 2 s . If we consider that Λ QCD ∼ mα 2 s , the gluonic correlator in Eq. (37) cannot be computed using perturbation theory. Therefore, in a model independent approach, one can leave it as a free parameter and fix it with experiment at some scale ν us (since the running of Eq. (37) with ν us is known, one can then obtain its value at another scale).
If we consider that mα 2 s ≫ Λ QCD , Eq. (37) can be computed perturbatively. Since mα 2 s is the next relevant scale, the effective role of Eq. (37) will be to replace ν us by mα 2 s (up to finite pieces that we are systematically neglecting) in Eq. (35). Then Eq. (38) reduces to Eq. (35) with ν us ∼ mα 2 s . In particular, we take ν us = −E n . As expected, Eq. (35) with ν us = −E n reproduces the already known O(mα 5 s ln α s ) correction [18] (see also [25,10]). Since in this situation one is assuming that Λ QCD mα 2 s ≪ 1, one can expand on this parameter. Therefore, non-perturbative corrections can be parameterized by local condensates. The leading and next-to-leading non-perturbative corrections have been computed in the literature [26,27].
There already exists an evaluation [10] within the vNRQCD framework [11,12,13] of the RG improved Heavy Quarkonium mass when Λ QCD ≪ mα 2 s . We agree for the spindependent terms (since we agree with the spin-dependent potentials computed in Ref. [12],