Inclusive Decays of Heavy Quarkonium to Light Particles

We derive the imaginary part of the potential NRQCD Hamiltonian up to order 1/m^4, when the typical momentum transfer between the heavy quarks is of the order of Lambda_{QCD} or greater, and the binding energy E much smaller than Lambda_{QCD}. We use this result to calculate the inclusive decay widths into light hadrons, photons and lepton pairs, up to O(mv^3 x (Lambda_{QCD}^2/m^2,E/m)) and O(mv^5) times a short-distance coefficient, for S- and P-wave heavy quarkonium states, respectively. We achieve a large reduction in the number of unknown non-perturbative parameters and, therefore, we obtain new model-independent QCD predictions. All the NRQCD matrix elements relevant to that order are expressed in terms of the wave functions at the origin and six universal non-perturbative parameters. The wave-function dependence factorizes and drops out in the ratio of hadronic and electromagnetic decay widths. The universal non-perturbative parameters are expressed in terms of gluonic field-strength correlators, which may be fixed by experimental data or, alternatively, by lattice simulations. Our expressions are expected to hold for most of the charmonium and bottomonium states below threshold. The calculations and methodology are explained in detail so that the evaluation of higher order NRQCD matrix elements in this framework should be straightforward. An example is provided.


I. INTRODUCTION
Heavy Quarkonium is characterized by the small relative velocity v of the heavy quarks in their centre-of-mass frame. This small parameter produces a hierarchy of widely separated scales once multiplied by the mass m of the heavy particle: m (hard), mv (soft), mv 2 (ultrasoft), . . .. In general, we have E ∼ mv 2 ≪ p ∼ mv ≪ m, where E is the binding energy and p the relative three-momentum.
The use of NRQCD [1] allowed a factorization of the physics due to the scale m from the one due to smaller scales. Moreover, it allowed the description of heavy quarkonium inclusive decays into light fermions, photons and leptons, in terms of matrix elements of local 4-quark operators, in a systematic way. These 4-quark operators are of two types: colour-singlet and colour-octet operators. The matrix elements of the colour-singlet operators can be related in a rigorous way with quantum-field theory defined quarkonium wave functions [1]. Intuitively, these wave functions should be related with the wave functions that appear in a Schrödingerlike formulation of the bound-state system, namely two heavy quarks interacting through a potential. On the other hand, the colour-octet ones were thought to have no parallel in such formulation. In either case, even though there had been a lot of relevant work in obtaining the QCD potential in terms of Wilson loops [2], it was not known how to obtain the systematic connection between NRQCD and a Schrödinger-like formulation in the nonperturbative case, or even whether it existed and, if so, under which circumstances. Even in the perturbative case, for which expressions for the potential existed at lower orders in the past [3], a clean and simple derivation of such Schrödinger-like formulation incorporating perturbative ultrasoft gluons was not clear once higher order calculations in α s were required.
The observation that NRQCD still contains dynamical scales that are not relevant to the kinematical situation of the lower-lying states in heavy quarkonium (those with energy scales larger than the ultrasoft scale) [4] (see also [5]), paved the way towards the resolution of the questions above. Indeed, it was realized that further simplifications occur if we integrate them out, and the resulting effective field theory was called pNRQCD [4]. The degrees of freedom of pNRQCD depend on the interplay between the characteristic scales of the given non-relativistic system, namely E, p and the momentum transfer k, and the characteristic scale of non-perturbative physics in QCD, which will be denoted by Λ QCD . Therefore, how a Schrödinger-like formulation develops, and thus how the NRQCD 4-fermion matrix elements will show up within this framework, depends on the specific kinematic situation considered.
When the typical momentum transfer k is much larger than Λ QCD , k ∼ p ≫ E Λ QCD , the pNRQCD Lagrangian [4,6] contains not only the singlet field, which is also present in the Schrödinger-like formulation, but also the octet field, ultrasoft gluons and light quarks. The matching from NRQCD to pNRQCD (integration of the soft scale) can be done in perturbation theory. In nature, this situation is relevant to the Υ(1S) and the t-t production near threshold. If in addition E ≫ Λ QCD , we are entirely in the weak-coupling regime (E ∼ mα 2 s , p ∼ k ∼ mα s ) where non-perturbative effects can be parameterized by local condensates [7]. In this regime pNRQCD has been used to obtain the complete set of logarithmic corrections to the QCD static potential at three loops [8], the complete set of logarithmic corrections to the very heavy quarkonium spectrum at O(mα 5 s ) [9] (see also [10]), the resummation of logs at the same order [11,12] and, very recently, the (almost) complete spectrum of very heavy quarkonium at O(mα 5 s ) [13]. We can still use the same pNRQCD Lagrangian for systems with E Λ QCD . Then, however, some of the calculations in pNRQCD cannot be carried out perturbatively and the non-perturbative effects can no longer be parameterized by local condensates (see [7,6]). When the typical momentum transfer k Λ QCD and the binding energy is small, namely E ≪ Λ QCD , the degrees of freedom of pNRQCD are the singlet field and pseudo-Goldstone bosons (pions), if hybrids and other degrees of freedom associated with heavy-light meson pair threshold production develop a mass gap of O(Λ QCD ), as it is assumed in Refs. [14,15,6] and in what follows. If we ignore Goldstone bosons, which play a negligible role in the present analysis, we recover the celebrated Schrödinger-like picture of quark and antiquark interacting through a potential. Therefore, the pNRQCD Lagrangian reads [14,6]: where h is the pNRQCD Hamiltonian, to be determined by matching to NRQCD. In general, one should be able to obtain the binding energies and the total decay widths from the real and imaginary parts of the complex poles of the propagator. At the accuracy we are aiming at in this paper the total decay width of the singlet heavy-quarkonium state may be defined as: Γ = −2 Im n, l, s, j|h|n, l, s, j , where |n, l, s, j are the eigenstates of the real part of the Hamiltonian h.
In this paper we will be concerned with this situation and will consider in full detail not only the calculation in the general case (A) Λ QCD k (Section III), but also the particular situation (B) Λ QCD ≪ k (Section V): A) Λ QCD is smaller than or of the order of k. In this case, the (non-perturbative) matching to pNRQCD has to be done in a single step. This case has been developed in a systematic way in Refs. [14,15]. As a consequence the complete set of potentials up to order 1/m 2 could be finally calculated [14,15], including a 1/m potential, which had been missed so far and completing (and in some cases correcting) the previous expressions obtained in the literature [2] for the 1/m 2 potential. Most of the charmonium and bottomonium states below threshold are expected to be in this situation. B) Λ QCD is much smaller than the typical momentum transfer k. In this case, the degrees of freedom with energy larger than or similar to k can still be integrated out perturbatively. This leads to an intermediate EFT that contains, besides the singlet, also octet fields and "ultrasoft" gluons (meaning gluons with energies Λ QCD here) as dynamical degrees of freedom [4,6]; it has the same Lagrangian as pNRQCD in the weak coupling regime. We will call this EFT pNRQCD ′ . 1 The octet and "ultrasoft" gluon fields are eventually integrated out by the (non-perturbative) matching to pNRQCD [6].
In either case, it remained to be seen how the matrix elements of the 4-fermion operators are encoded in this formulation. This was especially needed for the octet ones since, as mentioned before, it was thought that they could not be accommodated in a Schrödinger-like formulation. However, in [16], we have shown that, by using pNRQCD, it is indeed possible to relate the matrix elements of the colour-octet operator with the wave function at the origin and additional bound-state independent non-perturbative parameters. This was done for the specific case of P -wave quarkonium decays. Here, we will apply the same method to express all the NRQCD matrix elements relevant to inclusive S-wave quarkonium decays into light hadrons, photons and lepton pairs at O(c(α s (m))mv 3 × (Λ 2 QCD /m 2 , E/m)) (c(α s (m)) being a function of α s (m) computable within perturbation theory). This reduces the number of unknown parameters for the total decay widths of charmonium and bottomonium states below threshold by roughly a factor of 2, which allows us, in turn, to formulate several new model-independent predictions. Particularly important is the fact that our formalism allows the physics due to the solution of the Schrödinger equation, which appears entirely in the wave-function, to be disentangled not only from the short-distance physics at scales of O(m), but also from the gluonic excitations with an energy of O(Λ QCD ). As a consequence, the wavefunction dependence drops out in the ratio of hadronic and electromagnetic decay widths. For this class of observables the reduction in the number of non-perturbative parameters in going from NRQCD to pNRQCD is even more dramatic, since only the (six) non-perturbative universal parameters appearing at this order in pNRQCD are needed.
Finally, we would like to mention the dynamical situation when the binding energy is positive and of the same order of magnitude as the momentum transfer k, namely when E Λ QCD ∼ k. In this case degrees of freedom with energy ∼ Λ QCD cannot be integrated out. States close to and beyond heavy-light meson pair threshold are expected to be in this situation. The results of this paper do not apply, in principle, to this case.
The paper is arranged as follows. Section II reviews some aspects of NRQCD that are relevant to the rest of the paper. Section III provides a detailed description of the computation of the "spectrum" of NRQCD, in particular the ground state, in the 1/m expansion in the general case. It is meant for the reader interested in learning the techniques involved in this type of computations. The description of pNRQCD, its power counting and the relation between the computation of Section III and the Hamiltonian in pNRQCD are given in Section IV. Section V provides a detailed description of the matching between pNRQCD and NRQCD in the particular case E ≪ Λ QCD ≪ k. This section may help the reader who is not willing to go through the general case in Section III, but still wants to get a flavour of the kind of calculations we are performing. Section VI summarizes our results. The reader who is only interested in our final results and wants to skip any computational detail may jump directly to this section. Section VII displays some modelindependent predictions that follow from our results. We finally draw our conclusions in Section VIII. A number of appendices complement the main body of the paper. Appendix A recalls the 4-fermion NRQCD operators at O(1/m 4 ). Appendix B gives the general formula relating an arbitrary NRQCD matrix element with the computation in pNRQCD. Appendix C gives the leading log renormalization group running of the imaginary parts of the 4-fermion NRQCD operators matching coefficients. Appendix D shows how to deal with ill-defined products of distributions within dimensional regularization. Appendix E shows how unitary transformations can relate different forms of the pNRQCD Hamiltonian.

II. NRQCD
NRQCD is obtained from QCD by integrating out the heavy-quark mass scale m [1]. The NRQCD Lagrangian can be written as follows: where L g involves only gluon fields, L light involves light-quark and gluon fields, and L 2n−f are the terms in the Lagrangian with 2n heavy-quark fields. The NRQCD Lagrangian can be organized as a series expansion in α s (m) and in the inverse of the heavy-quark mass 1/m. Powers of α s (m) are encoded into the Wilson coefficients of NRQCD.
In this paper, we aim at a description of heavy quarkonium inclusive decays into light hadrons and electromagnetic decays, whose appearance is due to the imaginary terms of the NRQCD Lagrangian. It is convenient, then, to split the Lagrangian into the Hermitian (real) and the anti-Hermitian (imaginary) part: where and The operators responsible for heavy-quarkonium decays are the NRQCD 4-fermion operators whose matching coefficients carry an imaginary part. For our purposes, it is sufficient to consider either dimension 6 or dimension 8 4-fermion operators: With the superscript EM, we indicate operators responsible for the electromagnetic decays. More explicitly, we have The definitions of the hadronic operators can be found in [1]. For ease of reference, we recall them in Appendix A, where we also give the definitions of the electromagnetic operators. The distinction between hadronic and electromagnetic operators is somewhat artificial. In general the 4-fermion operators listed in Eqs. (A1)-(A18) are all the dimension 6 and 8 operators needed to describe decays into light hadrons and/or hard electromagnetic particles. The information needed in order to describe decays into hard electromagnetic particles is encoded into the electromagnetic contributions to the matching coefficients. The electromagnetic operators defined in [1] arise from singling out the operators accompanying the matching coefficients whose imaginary parts correspond to pure electromagnetic decays and inserting into them the QCD vacuum (|vac vac|). This insertion guarantees that when calculating with these operators in NRQCD, no contamination from soft strong interactions will occur. Hence, the electromagnetic operators encode all the relevant information needed in order to calculate the quarkonium total decay width to electromagnetic particles only. However, one might also be interested in the decays to hard electromagnetic particles and soft light hadrons. In this case, the complementary to the above projector, namely 1 − |vac vac| should be considered. In this paper, however, we will restrict to the processes, and therefore to the operators, originally considered in [1].
The Hermitian piece of the NRQCD Lagrangian can also be written in a 1/m expansion: At order 1/m the different pieces of Eq. (5) read where ψ is the Pauli spinor field that annihilates the fermion and χ is the Pauli spinor field that creates the antifermion, iD 0 = i∂ 0 − gA 0 , iD = i∇ + gA, B i = −ǫ ijk G jk /2; for later use, we also define E i = G 0i and [D·, E] = D · E − E · D. The chromomagnetic matching coefficient c F is known at next-to-leading order and its value can be found in [17]. Concerning the explicit expression of the O(1/m 2 ) Lagrangian, see Ref. [15] for the operators without light quarks and Ref. [18] for the operators including light fermions.
The expansion up to O(H I ) was considered in [14] in order to obtain the 1/m potential. The O(H 2 I ) term was obtained in [15]. The O(H 3 I ) expression is new. A detailed derivation of Eq. (27) will be given in Sec. III D.
The states can also be formally expanded in 1/m: It is convenient to write the above states in terms of some new states |ñ; x 1 , x 2 , defined recursively as (see also Ref. [15]): As a consequence of Eq. (29), it holds that or equivalently (this equation will become crucial in later sections to simplify some calculations) At O(1/m), we obtain At O(1/m 2 ), we obtain where and the second term, due to the normalization of the state, reads (note that N 0 = 1 + N (2) 0 /m 2 + . . . is Hermitian): By using Eq. (15) and the identities obtained in Refs. [14,15], explicit expressions for the above Eqs. (32) and (33) can be worked out. In particular, at order 1/m we obtain (the spin-independent part was first obtained in [14]): where |n (0) stands for a shorthand notation of |n; x 1 , x 2 (0) , the state that encodes the gluonic content of the state |n; x 1 , x 2 (0) and is normalized as (0) n|m (0) = δ nm (for a precise definition, see Eq. (53) and the following discussion). We will use expression (36) in the subsequent sections. In this paper, we are interested in computing Im E n (actually Im E 0 ) with relative accuracy O(1/m 2 ). We will now explain in detail how the different terms of Eq. (27) appear within the quantum-mechanical calculation.
Equations. (24)- (26), as well as the analogous equations in Ref. [15], implicitly assume that the Hamiltonian is Hermitian. This is not true at arbitrary orders and the iteration of imaginary-dependent terms may lead to problems. Nevertheless, at the relative O(1/m 2 ) accuracy we are aiming at in this paper for the imaginary terms and for the n = 0 state, such effects are zero. Therefore, effectively, we only have to compute the expectation value of the imaginary part of the NRQCD Hamiltonian in terms of the O(1/m 2 ) eigenstates of the Hermitian part of the NRQCD Hamiltonian. 3 The reason is that the only imaginary contribution to the states up to O(1/m 2 ) comes from the first line of Eq. (34) and this term is zero for n = 0 because of the subsequent Eq. (69).
The imaginary terms in the NRQCD Lagrangian only appear in the matching coefficients of the 4-fermion operators, i.e. in L 4−f . Therefore, the imaginary part of the NRQCD Hamiltonian has the structure of Eq. (16). Profiting from this structure of the imaginary terms and since the iteration of the leading imaginary terms gives zero, Im E 0 can be computed from Expanding in 1/m the states and Im H, we can identify the different terms of Im E 0 in the 1/m expansion: They read After an explicit calculation, we have Im E since (1) Moreover, we have These results follow from Eq. (31), supplemented by the following argument. The colour structure of Im H (2) is such that, at the gluonic level, the following matrix elements are produced within the total expression: precision is warranted.
(by definition) and In order to deal with this second expression, we note that the lowest excitation, in the limit x 1 → x 2 , has no gluonic content and behaves like |0; x 1 , x 2 (0) = 1 l c / √ N c |vac , so that: where C f = (N 2 c − 1)/(2N c ). The above expressions may appear problematic since they involve the behaviour of the state in the limit x 1 → x 2 . and some regularization could be required in this case. However, we actually only need a weaker condition to ensure that Eq. (44) is zero. What we have is an expression like where O 1 is some unspecified operator. Following Ref. [14], this expression is the spectral decomposition of the Wilson loop (for the definition of a Wilson loop with a number n of operator insertions, see Ref. [15]): where O stays for the insertion of the operator O on a static Wilson loop of spatial extension x 1 − x 2 . In the presence of more operators, the symbol · · · c indicates the connected part (see in particular the erratum of Ref. [15]). One can see that the operator (49) is zero in the limit x 1 → x 2 . In order to obtain this result, it is very important that the delta acts directly on the states 4 . In this situation, and in the limit x 1 → x 2 , one can see that the disconnected piece of the Wilson loop cancels with the connected piece, proving Eq. (44). For the other terms, we have 4 We may have situations where the Wilson-loop operator has the structure In this case the argument does not apply since the delta does not act directly on the Wilson loop.
Indeed, the last two equations hold as well for an arbitrary n and not only for the state n = 0, for which we have explicitly displayed them. It can be easily checked that the imaginary part of Eq. (27) for n = 0 coincides with the above expression (38) supplemented by Eqs. The expressions obtained in the previous section can be rearranged in terms of the pure gluonic content (see Refs. [14,15]). In order to achieve this we have to make the quark field content of the states explicit and use the Wick theorem. There is some freedom in choosing the specific realization of the quark fields under spin transformations. In [14], the following state was chosen In the basis of 4-fermion operators that we are using in this paper (see Appendix A) and in the above basis, the quantum-mechanical operators that naturally appear are 1 l s ⊗ 1 l s and σ i ⊗ σ j , where 1 l s (σ i ) is the identity (sigma-matrix) in spin space acting either on the final-or the initial-spin quark-antiquark state. Analogous definitions can be made for the operators acting on the colour subspace. Another possibility is the state which has been used in Ref. [15]. The quantum-mechanical operators, which naturally appear in this way, are 1 l 1,2 , σ i 1,2 , and they represent the operators acting either on the particle 1 or 2 (in this case we have always a particle interpretation). Analogous definitions can be made for the operators acting on the colour subspace. This representation appears to be more convenient for the calculations of the quantum-mechanical matching. In principle, one could also write the local 4-fermion operators in a basis convenient for these states by using Fierz transformations [20].
In both cases, we assume the state to be properly normalized in the spin sector. Depending on the calculation, one definition turns out to be more useful than the other. In any case, at the end, we are interested to write the quantum-mechanical Hamiltonian relevant to the Schrödinger equation. A way of avoiding ambiguities is to write everything in terms of a definite set of spin operators. We will adopt the operators S i and 1 l acting on a generic 1/2 ⊗ 1/2 spin space and defined as It is possible to transform them in the operators 1 l s ⊗ 1 l s and σ i ⊗ σ j by using the identities: Let us now compute the different matrix elements that appear in Eq. (41). The contribution due to the dimension-8 4-fermion operators reads where Equations. (58)-(61) and (65) provide the explicit expressions of the operators T S and T ij SJ first used in Ref. [16]. The non-perturbative constant E 1 (as well as all the other constants and E (2,norm.) 3 appearing in this section) will be defined in Sec. III F. If we consider the electromagnetic contribution due to H (4) 4−f , we obtain (in this case there are no octet operators) In order to calculate the contribution due to the 1/m correction to the state, we need to know (a 1 l is understood where no spin-operator is displayed): where F j ≡ F (x j ), F being a generic gluonic operator. In particular from the last equation it follows that It is this equation that guarantees that, for the n = 0 quarkonium state, no imaginary contribution is carried by the state (see the discussion at the beginning of Sec. III D). Finally, from the above equations it follows that the contributions due to the 1/m correction to the state read: For the electromagnetic contribution we have the intermediate vacuum, which does not allow an intermediate emission of gluons. This means that The contributions due to the normalization of the state read Exactly the same contribution is obtained from the electromagnetic terms if we change the subscript 1 in the matching coefficients by EM.

F. Gluonic correlators
The non-perturbative constants E n , B n , E 3 , E , which appeared in the previous section, are pure gluonic quantities, since the fermionic fields have been integrated out. Within the quantum-mechanical matching, they are first obtained in terms of gluonic states. For instance, we obtain the expressions · · · · · · · · · .
For the first two equations, there is no need to specify whether the gluonic fields are inserted on the particle or on the antiparticle line since they give the same contribution. We do not give here the complete list of expressions at the quantum-mechanical level, since this section does this in terms of Wilson loop operators. The former may be derived straightforwardly from the latter by spectral decomposition.
Using the techniques of Refs. [14,15], it is possible to express E n , B n , E 3 , E in terms of the more familiar gluonic field correlators. We obtain (traces as well as suitable Schwinger lines connecting the gluon fields are understood if not explicitly displayed) 5 where and similarly for the other structures with four chromoelectric fields that appear in Eqs.
For further use, we also define: IV. PNRQCD

A. Matching to pNRQCD
Expressions (27) and alike are no more than formal expansions in H I , i.e. in 1/m, until some dynamical assumption is made. We will assume a mass gap of order Λ QCD ≫ mv 2 between the lowest-lying excitation and the higher ones. Under this assumption all the excitations (n = 0) decouple from the ground state (n = 0), which is identified as the only degree of freedom of pNRQCD. It corresponds to the singlet state S in the pNRQCD Lagrangian (1). Moreover, the above expansion acquires a dynamical meaning, becoming an expansion in Λ QCD /m and v in the effective field theory.
The above assumption is the same as was made in Refs. [14,15] in the situation without massless fermions. In this work, we are including light fermions. Nevertheless, at least in this paper, we will assume that this does not change the structure of the leading order solution (this was also assumed in Ref. [16]). In other words, we will assume that the size of the typical splittings between the ground state (heavy quarkonium) and the gluonic excitations (hybrids) is much larger than the typical splittings produced by the solutions of the Schrödinger equation for the heavy quarkonium. This is, indeed, supported by lattice simulations where the plots of the static potentials for the heavy quarkonium and hybrids show the same pattern after the inclusion of light fermions [21]. Nevertheless, in principle, a new problem may arise. Once light fermions have been incorporated into the spectrum, new gauge-invariant states appear besides the heavy quarkonium, hybrids and glueballs. On the one hand, we have the states with no heavy quark content. Due to chiral symmetry, there is a mass gap, of O(Λ χ ), between the Goldstone bosons, which are massless in the chiral limit, and the rest of the spectrum. We will consider that the Goldstone bosons are ultrasoft degrees of freedom and that Λ χ ∼ Λ QCD , so that the rest of the spectrum should be integrated out. Besides these, we also have bound states made of one heavy quark and light quarks. In practice, we are considering the Qq-Qq system. The energy of this system is, according to the HQET counting rules [22]: Therefore, sinceΛ ∼ Λ QCD , we will assume that they also have to be integrated out. Problems may appear if we try to study the heavy quarkonium near threshold. In this case there is no mass gap between the heavy quarkonium and the creation of a Qq-Qq pair. Thus, if we want to study the heavy quarkonium near threshold, we should include these degrees of freedom in the spectrum (for a model-dependent approach to this situation see, for instance, [23]). We will not do so in this paper. It may happen, however, that the mixing between the heavy quarkonium and the Qq-Qq is small. Indeed, such a mixing is suppressed in the large N c counting. Summarizing, light fermions contribute within this picture in three ways: 1) hard light fermions, they are encoded into the matching coefficients of the NRQCD Lagrangian and obtained from the computation of perturbative Feynman diagrams at the scale m; 2) soft light fermions, a term that denotes, in a generic way, all the fermions that are incorporated in the potentials; it is expected that their main effects can be simulated by a variation of the value of the parameters in the potentials; 3) ultrasoft light fermions, these are the ones that will become pions and, since they are also ultrasoft degrees of freedom, they should be incorporated in the effective Lagrangian together with the heavy quarkonium. However, we will not consider them in the present paper, even if we do not expect to find conceptual problems in an eventual incorporation.
In conclusion, the matching condition to pNRQCD for the real part reads At O(1/m) the matching has been performed in Ref. [14] and at O(1/m 2 ) in Ref. [15] (for the case without light fermions). We refer to those articles for further details about the structure of the potentials. For the imaginary piece, we have the analogous matching condition: Using the results of the previous sections, we can now write the first two terms in the 1/m expansion of Im h (the P -wave-dependent terms were obtained in Ref. [16]): where The above expressions have been given in 4 dimensions. Therefore, they should be generalized to d dimensions if we want to work in an MS-like scheme in order to use the same scheme as for the NRQCD matching coefficients computation. This becomes relevant when logarithmic ultraviolet divergences appear in the non-perturbative constants. Hence, eventual lattice calculations must be converted to MS in this case. Nevertheless, in several situations, it is not necessary to work in an MS scheme if we only want to obtain the non-perturbative objects from experiment, since the scheme dependence simply goes into a redefinition of the non-perturbative constants. Finally, note also that in addition to the divergences in the non-perturbative constants, which are due to large momentum transfers k, at some point there will also be ultraviolet divergences arising in quantum-mechanical perturbation theory, which are due to large relative momenta p. These must also be regulated in dimensional regularization and MS-subtracted, along the lines worked out in Ref. [24].

B. Power counting in pNRQCD
With the above results, we are in a position to compute the inclusive decays of heavy quarkonium into light particles by using Eq. (2). Before doing so, we have to specify some power-counting rules in order to estimate the importance of the different terms of the pNRQCD Hamiltonian. Previous discussions on this subject, some of which we will repeat here, can be found in Refs. [14,15].
With the results of Sec. IV A and using Eq. (2), the decay width of S-wave quarkonium has schematically the following structure: where c 4−f stands for the NRQCD 4-fermion matching coefficients and R ns0s is the S-wave radial component of the solution of the real piece of the Schrödinger equation: with the normalization (|s spin denotes the normalized spin component): Although E njls coincides with the binding energy of the system at the order we are working at, it will no longer be so when iterations of imaginary parts start playing a role.
From Eq. (90), we can see how the power counting has to be organized. On the one hand, we have an explicit expansion in Λ QCD /m, independent of the details of the bound state. In the most conservative situation (Λ QCD ∼ mv), it would correspond to having the power counting Λ QCD /m ∼ v. We can also find derivatives of the wave function divided by m. They typically scale like ∇/m ∼ v. On the other hand, the normalization condition of the wave function sets the scaling |R njls | 2 ∼ (mv) 3 . This means that a formal O(mv 5 ) accuracy (leaving aside possible α s (m) suppressions due to the NRQCD matching coefficients) is achieved with Eq. (90). At the same order of accuracy, the decay width of P -wave quarkonium has the structure: In the above discussion, we have only considered the leading order power counting of the wave function at the origin ∼ (mv) 3 . This accuracy is sufficient for the P -wave function of Eq. (93), as well as for the wave functions multiplying Λ 2 QCD /m 2 terms or with two ∇ in Eq. (90) but not for the leading order term. In this case, one has to take into account that the wave function at the origin also has subleading contributions in v: |R njls (0)| 2 ∼ (mv) 3 (1 + av + bv 2 + · · ·). Therefore, we have to further specify the solution of Eq. (91), for which we have to set the power counting of the potentials in the Schrödinger equation. Since we do not know the specific dynamics of the different potentials, the only thing we can do is to require consistency of the theory and allow, in principle, the most conservative counting. This would correspond to setting the counting by the largest scale that has been integrated out, i.e. the potentials would scale like (mv) d , d being their dimension. 6 For definiteness, we will also assume α s (m) ∼ v q with q > 1.
Leading order. Consistency of the theory requires the virial theorem to be fulfilled. In other words, the potential at leading order needs to fulfil with the power counting It follows that V (0) ∼ mv 2 (even if, using the most conservative power counting, we would have obtained V (0) ∼ mv). Moreover, in our power counting we have V (1) /m ∼ mv 2 . 7 Therefore, in the most conservative situation, we would have The important point here is that, at this order, the potential is spin-independent (E nl ). Therefore, the leading-order P -wave function reads where |r is the normalized eigenstate of the position and |js stay for J (total angular momentum) and S eigenstates such that where m denotes the third component of the angular momentum and detailed expressions of Y 1 jm (r) can be found in Ref. [26], Appendix B.
Next-to-leading order. The O(1/m 2 ) potential scales utmost as V (2) /m 2 ∼ mv 3 . Therefore, in the most conservative situation, we would have 6 Notice that our power counting rules are different from those of [1,25]. Whereas ours are meant to apply in the situation Λ QCD ≫ mv 2 , the power counting rules in Refs. [1,25] rather follow the counting Λ QCD ∼ mv 2 . Indeed, if we take Λ QCD ∼ mv 2 in our results we obtain a similar power counting for the NRQCD matrix elements. 7 As a consequence, if the potential V (1) is non-perturbative, we have no general argument to consider V (1) /m subleading with respect to V (0) . A lattice simulation or some model-dependent studies are, therefore, highly desirable to discern the issue. Whereas it is difficult to obtain this information from the spectrum structure, the study of the decays may perhaps shed some light on this problem. Finally, we note that, in the perturbative situation, V (1) has an extra α 2 s suppression. Further discussions can be found in Ref. [14].
At this order, spin-dependent contributions start to appear. In particular, the spindependent potential contributing to the S-wave function at the origin reads where [15] Re V (1,1) This potential produces the following correction to the S-wave function: where P n ≡ I − |n0 n0| and r|njls = φ n0 (r)). If the spin-dependent potential (100) is O(mv 3 ), it just provides the leading order spindependent correction to the S-wave function at the origin and one can use the difference between vector and pseudoscalar decays to fix the value of the correction. If the spindependent potential is O(mv 4 ), it provides a correction to the S-wave function squared at the origin, which is of the same order as the O(v 2 ) corrections to the decay width that we have already evaluated. Therefore, in this last situation, Eq. (102) would account for the full difference between the vector and pseudoscalar wave functions at the origin at relative order O(v 2 ), which is the precision we are aiming at in this work. This last counting seems to be supported by the size of the spin-dependent splittings in the bottomonium and charmonium spectra.
For the spin-independent contributions, we will make no assumption at this or higher orders, as their effects will be encoded into the wave functions, which will be left unevaluated. Our results allow for the most conservative counting where V (1) /m ∼ mv 2 and V (2) (spin-independent)/m 2 ∼ mv 3 . We note that, in this power counting, potentials with imaginary part arise in the pNRQCD Hamiltonian at order mα s (m) 2 v 3 (where the powers in α s (m) come from the imaginary part of the 4-fermion matching coefficients in NRQCD). Therefore, corrections due to the iteration of imaginary terms, which could affect the validity of Eq. (2), are far beyond the accuracy of this paper. In fact, the general factorization formula put forward in [1] may not hold beyond a certain order.
In any case, we do not rule out that a different power counting may also lead to consistent equations in the non-perturbative regime for some specific ratios of Λ QCD versus m and versus p and k. This point deserves further investigation and may lead to a different implementation of the matching procedure. We recall that the issue of assessing the power counting in the non-perturbative situation has been addressed before by Beneke [27] and by Fleming et al. [28]. In both cases, the authors have given some freedom to the possible size of the NRQCD matrix elements by introducing a parameter λ that interpolates between the power counting in the perturbative limit and other possible power countings according to the value of λ. In this respect, our formalism may shed more light to clarify this problem, since it incorporates the factorization between the soft and the ultrasoft scales, allowing us to write the NRQCD matrix elements in terms of the wave function at the origin and of some boundstate-independent constants. Another point of concern is whether there are non-perturbative effects that are not accounted for in the 1/m matching.
We conclude this section by giving a useful equality, valid in dimensional regularization: which follows from the fact that we know the behaviour of the potential and the wave function (up to a constant) at short distances and that (see Appendix D) n, j, l, s|V (0) |r = 0 = n, j, l, s|V (1) With this we have discussed the relative importance of the different terms that will appear in our evaluation of the decay widths. The results can be found in Sec. VI.
Although it is not clear whether quarkonia states fulfilling mv ≫ Λ QCD ≫ mv 2 exist in nature (mv ∼ k ∼ p and mv 2 ∼ E will always be understood in the present section), this situation is worth investigating for several reasons. First of all, the calculation in the general case of the Sec. III is non-standard and, hence, any independent check of it, even if it is in a particular case, is welcome. Secondly, the calculation in this case can be divided into two steps. The first step can be carried out by a perturbative calculation in α s , which involves far more familiar techniques. The second step, even if it is non-perturbative in α s , admits a diagrammatic representation, which makes the calculation somewhat more intuitive. Third, the more detailed information on the potential allows us to make important tests on how the terms in the potential can be consistently reshuffled by means of unitary transformations [14], as is illustrated in the example provided in Appendix E.
A. pNRQCD ′ As mentioned in the Introduction, we shall call pNRQCD ′ the EFT for energies below mv. Since mv ≫ Λ QCD , the integration of the energy scale mv, namely the matching between NRQCD and pNRQCD ′ can be carried out perturbatively in α s . This is done following Refs. [4,6]. A tree-level matching is sufficient, but higher orders in the multipole expansion will be needed. We only display below the terms eventually required in the calculation: where the traces are over colour space only. S and O are chosen here to transform as a 1/2 ⊗ 1/2 representation in spin space (hence σ 1 − σ 2 = σ 1 ⊗ 1 l 2 − 1 l 1 ⊗ σ 2 ); h s and h o read as follows (again we only display terms eventually required in the calculation): The Feynman rules associated to this Lagrangian are displayed in Fig. 1.

B. Matching pNRQCD to pNRQCD ′
The matching of pNRQCD ′ to pNRQCD can no longer be done perturbatively in α s , but it can indeed be done perturbatively in the following ratios of scales: Λ QCD /mv (multipole expansion), Λ QCD /m and mv 2 /Λ QCD . The diagrams contributing to the calculation are displayed in Figs. 2 to 9.
We have focused on contributions to S-wave states involving imaginary parts. Since the imaginary parts, which are inherited from NRQCD, sit on local (δ (3) (r), ∇δ (3) (r)∇, etc.) terms in the pNRQCD ′ Lagrangian, they tend to cancel when multiplied by the r's arising from the multipole expansion. Hence, for an imaginary part to contribute, it is necessary to have a sufficient number of derivatives (usually arising from the mv 2 /Λ QCD expansion) as to kill all the r's. Since derivatives are always accompanied by powers of 1/m, it implies that at a given order of 1/m, only a finite number of terms in the multipole expansion contribute. In our case a fourth order in the multipole expansion is sufficient. The natural way to organize the calculation in our case would be to assign a size mv p to Λ QCD , 1 < p < 2 and v q to α s , 1 < q < 2, and to carry out the calculation at the desired order in v. However, our main goal here is not the phenomenological relevance of the situation mv ≫ Λ QCD ≫ mv 2 , but providing an independent calculation to support the results of Sec. IV A. Hence, irrespectively of what p and q may be, we will only be interested in fishing up the imaginary pieces that contribute to S-wave states up to order 1/m 4 .
The two diagrams in Fig. 2 correspond to the leading contribution in the Λ QCD /mv and Λ QCD /m expansion, respectively. Figure 3 displays the evaluation of each of them in the mv 2 /Λ QCD expansion. The diagrams in Fig. 4 correspond to the next-to-leading order contributions in the Λ QCD /mv expansion, and Figs. 5-9 display their evaluation in the mv 2 /Λ QCD expansion. It is then clear that the basic skeleton of the calculation consists of the x = (Λ QCD /mv) 2 and y = (Λ QCD /m) 2 expansions, which suggests writing the pNRQCD Hamiltonian as: The interpolating fields of pNRQCD ′ and pNRQCD will be related by: The matching calculation reads: The right-hand side of the matching calculation has the following structure (up to a global i factor, which is dropped): Hence, once we have made sure that, up to contact terms, the left-hand side of Eq. (110) has exactly this structure, we can easily identify the contributions to the pNRQCD Hamiltonian from the second term of the expression (111).

C. Calculation
Let us then proceed to the calculation of the left-hand side of Eq. (110) (in order to match Eq. (111) a global i factor will also be dropped).
Diagram (a) of Fig. 2 gives: The fact that mv 2 /Λ QCD is small is implemented by expanding the exponential. This guarantees that we will eventually get usual, energy-independent, potentials.
The first contributions arise at O (mv 2 /Λ QCD ) from the O(1/m 4 ) P -wave (Fig. 3a) and S-wave (Fig. 3e) terms in the octet potential of (107): where T ij SJ are defined in Eqs.
Illdefined expressions arise in the calculation, from products of distributions (both products of two delta functions and products of delta functions with non-local potentials, which explode as r → 0). It is most convenient to use dimensional regularization in this case, which sets all these terms to zero. This is shown in Appendix D, where the relation with other regularizations is also discussed. Having this in mind, it is clear that, at the order we are interested in, Im (V o − V s ) 2 = 0 and Im (V o − V s )r = 0. Hence, we only have to consider: If we decide to take one power (E − h s ) to the right and one to the left we have: which does not produce any imaginary part. However, an equally acceptable expression is: which does produce an imaginary part. This apparent paradox only reflects the fact that expression (115) by itself (as well as some of the expressions we will find below) does not determine uniquely its contribution to the potentials. This expression always leads to contact terms, wave-function normalization and potential, as is apparent in (116) and (117), but depending on how we decide to organize the calculation, the terms associated to each of these pieces change. For instance, when matched to (111), (116) gives: whereas (117) gives: This should not be a surprise. It has already been discussed in Ref. [14] that this ambiguity exactly corresponds to the freedom of making unitary transformations in a quantummechanical Hamiltonian, and does not affect physical observables. This is discussed in detail in Appendix E for the decay widths of the S-wave states we are concerned with. In order to fix the contribution to the potential of any term once forever, we will use the following prescription. If we have an expression with singlet propagators 1/(E − h s ) only in the external legs, and an even number of powers of (E − h s ), we will take the one closest to the left propagator to the left and the one closest to the right propagator to the right, and repeat until no power is left except in contact terms. Accordingly, in the intermediate steps, when terms with a single external leg 1/(E − h s ) and several powers of (E − h s ) are produced, one should take these powers towards the 1/(E − h s ) leg until no power is left except in contact terms. If the number of powers of (E − h s ) is odd, we use the same prescription until a single power is left. We then write (E − h s ) = (E − h s )/2 + (E − h s )/2 and take one half to the right and one half to the left. Expressions with an internal singlet propagator also appear, which require a more careful treatment as will be discussed after (128) below. Note that this prescription to organize the calculation needs not coincide with the prescription for fixing the wave-function normalization in Sec. IV A. Hence, we only expect to agree with the results of that section up to a unitary transformation. Anyway, with this prescription, (115) gives rise to the potential obtained in (116) and hence to no imaginary part.
.. contribute, giving rise to: It is the first term in the fifth line that renders the contribution depicted in Fig. 3d: At O(m 4 v 8 /Λ 4 QCD ) and higher, only imaginary parts beyond 1/m 4 are produced. Consider next the diagram Fig. 2b. Since the chromomagnetic moment already provides two powers of 1/m, only the linear term in the expansion of the exponential contributes (Figs. 3b and 3c). This gives: Consider next Fig. 4a. Because of the four r's in the expression, only the following term in the expansion (E − h o ) 3 = (E − h s ) 3 + ... contributes. We obtain: For the symmetric diagram, we have: In fact both contributions are the same, adding up to (see formula i) above Eq. (15) of Ref. [14]): Consider next Fig. 4b. The only contributions come from .. in one octet propagator and 1 in the other. We obtain (Fig. 5b): Then consider Fig. 4c. From here we get several contributions. Because of the four r's we need a total of three powers of (E −h o ). When all the powers come from the octet propagator in the middle, we get contributions from ( The ones from the second term read (Fig. 6): When a power of (E − h o ) does not come from the octet propagator in the middle, all the powers can be substituted by (E − h s ). If we put these contributions together with the first term before (126), we obtain (Fig. 7): Consider next Fig. 4d. Clearly this diagram contains the iteration of lower order potentials, which must be isolated. This is achieved by adding and subtracting the projection operator into the gluonic ground state 1 = (1−|0 0|)+|0 0|. The piece (1−|0 0|) contains new contributions to the potential only, whereas the piece |0 0| contains both the iteration of lower order potentials and new contributions to the potential. Consider first the piece (1 − |0 0|). It is identical to Fig. 4c by taking V o → V s in the expression before (126) and changing the chromoelectric field correlators accordingly. We then have (Fig. 8): Consider next the contribution from |0 0|. The vacuum insertion leads to an internal singlet propagator. To be specific, we have: The exponentials of (E −h o ) will be expanded. In order to be consistent with the calculation of the lower order potentials and subtract only their iteration, we must treat the powers of (E − h o ) at each side of the internal singlet propagator exactly as we did in the calculation We can easily identify the contributions that match the following terms in (111): (131) We also see that, apart from the terms above, there are additional terms in (130) that may (and do) eventually lead to new contributions to the potential (none of them with imaginary parts for this example). For them we use the same prescription as stated at the beginning of the section. The contributions to the imaginary parts come from the following terms in the expansion only: (i) an (E − h o ) 4 from an octet propagator and a 1 from the other one ( Fig.  9, first diagram), and (ii) an (E − h o ) 3 from an octet propagator and an (E − h o ) from the other one (Fig. 9, all of them). They read:

D. Results
Combining all the above calculations we obtain the same result as in Sec. IV A, except for the terms proportional to Im( 2S+1 S S ). With the mere replacement: where we have defined: and the same expressions apply. As mentioned before, the difference is due to the different prescription to fix the wavefunction normalization in Sec. III B. In Appendix E we show that there exists an unitary transformation such that our results can be taken to the form of Sec. IV A, and hence they are equivalent for all purposes.
In fact, it is somewhat surprising that the two calculations lead to identical results (up to a unitary transformation). On general grounds, one could only expect that the result in this section be a particular case of the general results of Sec. III. In fact the real parts of the potentials in h are indeed particular cases of the potentials in [15]. However, since we did not need their specific form at any stage we have not lost generality in our final expressions. More surprising is the fact that the matching coefficients of the terms in the multipole expansion in pNRQCD ′ (105) were only calculated at tree level here, whereas the expressions in Sec. III correspond to an all-loop result. This indicates that there must be a symmetry protecting these terms against higher-loop corrections, which may (or may not) be an extension of reparametrization invariance [29] or Poincaré invariance itself [30]. 8 In summary, we have presented in this section an alternative derivation of (141)-(146), which does not rely so heavily on the 1/m expansion. The matching from NRQCD to pNRQCD ′ , which can be done perturbatively in α s , can indeed be implemented in the 1/m expansion, as originally proposed [4], but it can also be done entirely in the framework of the threshold expansion [31,13], where the kinetic term is kept in the denominator for potential loop contributions and the on-shell condition is used (the results obtained in either way are related by local field redefinitions). The matching between pNRQCD ′ and pNRQCD is done in the Λ QCD /mv, Λ QCD /m and mv 2 /Λ QCD expansions. The approaches taken in these two steps are quite different from the strict 1/m expansion of Sec. III, and the coincidence of the results strongly supports their correctness.

VI. RESULTS
In this section we list our expressions for S-wave decays up to O(c(α s (m))mv 3 × (Λ 2 QCD /m 2 , E/m)) and for P -wave decays up to O(c(α s (m))mv 5 ). The P -wave decay widths were first obtained in [16] and are given here for completeness. The S-wave decay widths are new. In order to help the reader and for further convenience, we will start by recalling, at the same level of accuracy, the expressions of the decay widths as they are known from NRQCD. In the following we define the radial part of the vector S-wave function as ) and the radial part of the pseudoscalar S-wave function as is the derivative of the leading order P -wave function. The symbols V and P stand for the vector and pseudoscalar S-wave heavy quarkonium and the symbol χ for the generic P -wave quarkonium (the states χ(n10) and χ(nJ1) are usually called h((n − 1)P ) and χ J ((n − 1)P ), respectively).

A. Decay Widths in NRQCD
Including up to the NRQCD 4-fermion operators of dimension 8, the inclusive decays of heavy quarkonia are given by: At the same order the electromagnetic decays are given by:

B. Decay Widths in pNRQCD
Up to O(c(α s (m))mv 3 × (Λ 2 QCD /m 2 , E/m)) for S-wave and O(c(α s (m))mv 5 ) for P -wave, the inclusive decays of heavy quarkonia are given in pNRQCD by: At the same order the electromagnetic decays are given by:

C. NRQCD Matrix elements
By comparing the decay widths in NRQCD and pNRQCD we obtain the following dictionary between the matrix elements of NRQCD and the non-perturbative constants of pNRQCD, valid up to (once normalized to m) O(v 3 × (Λ 2 QCD /m 2 , E/m)) for the S-wave matrix elements and up to O(v 5 ) for P -wave matrix elements: (2,EM) 3 Any other S-wave dimension-6 matrix element is 0 at NNLO and any other S-wave dimension-8 matrix element is 0 at LO.
Equation (152) is worth emphasizing. It is of the singlet type but, because of the term proportional to E 1 , its leading contribution is not only proportional to what one would expect from a pure singlet potential model. In Ref. [32] the authors have also elaborated on Eq. (152). Within the context of NRQCD [1], they use the leading equations of motion 9 , the power-counting rules of [1,25] and some arguments to neglect some mass-like terms, which could be generated under regularization. They get where the term proportional to E 1 is missing. Nevertheless, this does not necessarily reflect any inconsistency in any of the derivations since, according to the (perturbative-like) powercounting rules of [1,25], the term due to E 1 would be subleading. In any case, it would be very interesting to see how a term proportional to E 1 would appear in the derivation of Ref. [32]. Here, we only would like to point out the possibility that an E 1 /m term may show up as a correction to the neglected mass-like term in Ref. [32]. Finally, let us note that in the n0 , both terms on the right-hand side of Eq. (152) are of the same order and contribute to the decay width at order c(α s (m))mv 5 . Phenomenologically this is particularly relevant to the case of pseudoscalar decays into light hadrons and to the electromagnetic decays. In the case of vector decays into light hadrons the contribution coming from the operator V Q (nS)|P 1 ( 3 S 1 )|V Q (nS) may not be so important since the matching coefficient Im g 1 ( 3 S 1 ) ∼ α s (m) 3 is suppressed by a factor α s (m) with respect to the other ones (with the exception of Im f 1 ( 3 S 1 ) and Im f 8 ( 3 P 1 ), which are also of order α s (m) 3 ).

D. Evolution equations
In [1] evolution equations for the 4-fermion operators have been obtained. If we focus on the states that we are studying in this paper, the following evolution equations for the NRQCD matrix elements are obtained: Since we have, at O(α s ) and leading log accuracy, Eqs. (147)-(150) are compatible with the evolution equations (159)-(162) at leading log accuracy. Note that at this order there is no ν dependence in the states, and hence the derivatives with respect to ν can be taken out of the expectation values. In Ref. [16] it was proved that Eq. (163) gives the correct running for the octet operator of Eq. (156). In Appendix C, the reader can find the evolution equations and their leading-order solutions for the imaginary parts of all the 4-fermion matching coefficients needed in this work.

VII. MODEL-INDEPENDENT PREDICTIONS
The inclusive decays of the heavy quarkonium (either hadronic or electromagnetic) are usually considered up to, and including, NRQCD matrix elements of 4-fermion operators of dimension 8. This means to consider the O(1/m 2 , 1/m 4 ) local 4-fermion operators of the NRQCD Lagrangian. With this accuracy, the decay into light hadrons of a vector Swave state is described in NRQCD by the matrix elements of two singlet operators (O 1 ( 3 S 1 ) and P 1 ), and three octet operators (O 8 ( 3 S 1 ), O 8 ( 1 S 0 ) and O 8 (P )). The corresponding pseudoscalar S-wave state decay needs, at the same level of accuracy, the additional knowledge of the matrix element of the singlet operator O 1 ( 1 S 0 ). The electromagnetic decays of the same S-states need the additional knowledge of the matrix elements of the singlet electromagnetic operators O EM ( 3 S 1 ) and O EM ( 1 S 0 ) respectively. The decay of a P -wave quarkonium state into light hadrons and the corresponding electromagnetic decay are described in NRQCD with the above accuracy by the matrix element of a singlet (O 1 (P )) and an octet (O 8 ( 1 S 0 )) operator. If we consider that in the bottomonium system in principle 14 S-and P -wave states lie below threshold (Υ(nS) and η b (nS) with n = 1, 2, 3; h b (nP ) and χ bJ (nP ) with n = 1, 2 and J = 0, 1, 2) and that in the charmonium system this is the case for 8 states (ψ(nS) and η c (nS) with n = 1, 2; h c (1P ) and χ cJ (1P ) with J = 0, 1, 2), all the bottomonium and charmonium S-and P -wave decays into light hadrons and into photons or e + e − are described in NRQCD up to O(1/m 4 ) by 46 unknown NRQCD matrix elements (40 for the S-wave decays and 6 for the P -wave decays). These matrix elements have to be fixed either by lattice simulations [33] or by fitting the data [34]. Only in the specific case of matrix elements of singlet operators does NRQCD allow an interpretation in terms of quarkonium wave functions and one can resort to potential models.
At the same level of accuracy S-and P -wave bottomonium and charmonium decays are described in pNRQCD, under the dynamical assumption Λ QCD ≫ mv 2 , by only 19 nonperturbative parameters. These are the 13 wave functions (one for each of the 10 S-wave quarkonium states below threshold, for which we need to distinguish different spin states, and a total number of 3 for the P -wave quarkonium states) and 6 universal non-perturbative parameters, which do not depend on the flavour and on the state (E 1 , 3 , E ).
In the above discussion we have counted NRQCD matrix elements by their dimensionality only. A more refined discussion would require that a maybe less conservative power counting be assigned to the NRQCD matrix elements as well as that the α s (m) suppression due to the short-distance NRQCD matching coefficients be taken into account. As we have already mentioned throughout the paper, the power counting of the NRQCD matrix elements is an open issue. To consider all the possibilities and phenomenological consequences goes beyond the scope of the present paper, whose aim is to set the theoretical framework. However, we would like to mention a few things. In the standard NRQCD power counting [25], the octet matrix elements are O(v 4 ) suppressed for S-wave decays if compared with the leading order. This is not so within our framework where, assuming the counting Λ QCD ∼ mv, they would only be O(v 2 )-suppressed. This is potentially relevant to Γ(V → LH) since Im f 1 ( 3 S 1 ) is O(α s (m))-suppressed with respect Im f 8 (S). In other words, the octet matrix element effects could potentially be much more important than usually thought for these decays. It would be interesting to analyse this possibility further.
The dramatic reduction in the number of parameters makes it possible, in the framework of pNRQCD, to formulate several new predictions with respect to NRQCD. In particular it is possible to relate information gained from decay widths of quarkonium with a specific flavour and principal quantum number to decay widths of quarkonium with different flavour and/or principal quantum number. Following this strategy in [16] the non-perturbative parameter E 3 has been fixed on the charmonium P -wave decay data and used to predict ratios of Pwave decay widths for the bottomonium system (in this case and at leading order there is no ambiguity on the relative size between the singlet and the octet contribution). Here we will concentrate on some exact model-independent relations valid for S-wave decays.
Let us consider the ratios of hadronic and electromagnetic decay widths for states with the same principal quantum number: Ten of these ratios exist, ten being the number of bottomonium and charmonium states below threshold. As we discussed above, in NRQCD, and if one includes all the NRQCD operators up to O(1/m 4 ), these 10 ratios are described by 40 non-perturbative parameters. It is a specific prediction of pNRQCD that, for the states for which the assumption Λ QCD ≫ mv 2 holds, the wave-function dependence drops out from the right-hand side of Eqs. n0 . In principle, if all the 10 bottomonium and charmonium S-wave states below threshold belonged to the dynamical regime Λ QCD ≫ mv 2 , then, in the framework of pNRQCD, the ratios of hadronic and electromagnetic decay widths would be described by the 6 non-perturbative universal parameters listed above only.
Particularly simple is, in pNRQCD, the expression of the ratios between R V n and R P n with different principal quantum number. We obtain up to order v 2 (with the counting , M(nS) being the meson mass): It is to be stressed that the octet-type contributions cancel (otherwise they would be 1/α s (m) enhanced in the vector case). This prediction should be compared with the one expected in NRQCD. Within the standard (perturbative-like) power counting, the same prediction is obtained in NRQCD. However, if one counts α s (m) ∼ v 2 as was done in [32], the contribution due to the octet matrix elements is of the same order as the corrections obtained above and it should be taken into account in the vector case. Therefore, in principle, one is able to check the theory and/or the power counting. As an example, taking m b = 5 GeV we get for the Υ(2S) and Υ(3S) states of the bottomonium system, R Υ 2 /R Υ 3 ≃ 1.3, which is close (within a 10% accuracy) to the experimental central value of about 1.4 that one can get from [35]. Let us also notice that, since Im m)), up to corrections of order v 3 we find that R P n , i.e. the ratio between hadronic and electromagnetic decay widths for pseudoscalar quarkonium, is the same for all radial excitations. However, it is not the purpose of this work to carry out a comprehensive and detailed phenomenological analysis, which is left to a subsequent publication.

VIII. CONCLUSIONS
We have obtained the imaginary part of the pNRQCD Hamiltonian up to O(1/m 4 ) in the non-perturbative regime (k Λ QCD ≫ mv 2 ). The expressions are given in Eqs. (87)-(89). As for any quantum-mechanical Hamiltonian, also the pNRQCD Hamiltonian is defined up to a unitary transformation. An alternative expression, related to the previous one by a unitary transformation, can be found in Sec. V D.
We have applied our results to calculate the inclusive decay widths to light hadrons, photons and leptons up to O(c(α s (m))mv 3 ×(Λ 2 QCD /m 2 , E/m)) for S-wave heavy quarkonium and up to O(c(α s (m))mv 5 ) for P -wave heavy quarkonium. These are given in Eqs. (141)-(146) and are the main result of the paper. An alternative way to present it is given in Sec. VI C, where all the NRQCD matrix elements entering in quarkonium decays up to this order are expressed in terms of the quarkonium wave functions at the origin and 6 nonperturbative gluonic correlators, which are flavour-and state-independent, and for this reason may be called universal. The wave-function dependence factorizes in all these expressions. It is particularly remarkable that this is also true for the octet matrix elements.
We have derived our expressions in two different ways: in Sec. III under the general assumption Λ QCD < ∼ k and in Sec. V under the particular assumption k ≫ Λ QCD . In the first case, we have matched directly NRQCD to pNRQCD in an entirely non-perturbative one-step procedure, based on the Hamiltonian formulation of NRQCD. In the second case, we have matched NRQCD to pNRQCD in a two-step procedure, the first perturbative, the second non-perturbative, but still with a clear diagrammatic interpretation based on the multipole expansion. The fact that these two completely different ways of deriving the pNRQCD Hamiltonian give the same answer up to a unitary transformation can be considered a stringent test on the correctness of the result. In Sec. VI D we have also checked that the evolution equations of our universal parameters are consistent at leading log accuracy with the known evolution equations of the NRQCD matrix elements.
In Sec. VII we have considered the phenomenological implications of our results. There exist 14 charmonium and bottomonium states below threshold. We expect our results to be applicable to most of these states. The exceptions are, on the one hand, the Υ(1S), which is commonly understood as a weak coupling state (i.e. k ≫ E Λ QCD ), and, on the other hand, states that are too close to the D-D threshold for charmonium or to the B-B threshold for bottomonium, like, maybe, the ψ(2S). Going from NRQCD to pNRQCD reduces the number of non-perturbative parameters needed to calculate the inclusive decay widths associated with these states by about a factor of 2. The situation is even better if we consider ratios of hadronic and electromagnetic decay widths. Since the wave-function dependence factorizes, it drops out in the ratios. It follows that only 6 universal parameters, which depend only on the light degrees of freedom of QCD, are needed. The already known data will be sufficient to fix all these parameters, to allow checks and to make new predictions. Moreover, suitable combinations of ratios give rise to novel parameter-free, model-independent predictions. We have considered some of them in Sec. VII.
The non-perturbative universal parameters that we have introduced do not necessarily need to be fitted to the experimental data. We have provided expressions for them in terms of correlators of gluonic fields. This allows for an eventual evaluation on the lattice. These parameters may also be obtained from QCD vacuum models [36]. We note that, once they become fixed, our results make the evaluation of NRQCD octet matrix elements possible from properties of the wave functions at the origin. Hence, any potential model that leads to definite wave functions [37] will provide definite results for these matrix elements. Nevertheless, it should be pointed out that, if we wish to obtain the NRQCD matrix elements given in Eqs. (147)-(150) with the aforementioned accuracy, any potential model to be used here must be consistent with the structure of the potential derived from NRQCD in terms of Wilson loops in Refs. [14,15]. In fact, the wave functions defined in this paper can also be computed in a model-independent way without resorting to data fitting. This is so because our wave functions correspond to the solution of a Schrödinger equation where the potentials are given in terms of expectation values of Wilson loops with suitable field insertions. Therefore, once lattice simulations are provided for the potentials [38], the wave function can be obtained unambiguously without any model dependence.
Since our method reduces the number of unknown parameters with respect to NRQCD, we expect it to become increasingly relevant as the number of needed NRQCD matrix elements increases. This seems to be necessary in the calculation of charmonium decay widths, where the non-relativistic expansion converges slowly. Indeed, higher order operators have been considered recently in Refs. [39,40]. In Appendix B, we give the general matching formula for the NRQCD matrix elements to the pNRQCD results without going through the whole matching procedure outlined in the main body of the paper.
We have also addressed, mainly in Sec. IV B, the issue of the power counting in NRQCD in the non-perturbative case. We believe that our formalism provides a suitable theoretical framework to study it. The power counting of NRQCD is not known a priori in the nonperturbative regime and it could, in principle, be different, depending on each dynamical system. This is particularly transparent in pNRQCD. There, the potentials are functions of r and Λ QCD . Therefore, as the typical value of r changes from system to system, one should accordingly assign a different size to each given potential. Moreover, having expressed the NRQCD matrix elements in terms of wave functions and universal correlators, we have disentangled the soft scale k, now entering in the wave function square, from the Λ QCD /m and E/m corrections. In fact, this is why we can construct ratios of convenient decay rates where the k-dependence drops, providing a more constrained set of relations. For these ratios the fixing of the power counting reduces to the evaluation of the correlators, while taking into account possible enhancement effects due to the NRQCD matching coefficients.
Finally, although in the present paper we have focused on inclusive decays to light hadrons, there should be no conceptual problem, a priori, in considering the NRQCD matrix elements that appear in heavy quarkonium production. We also expect there a significant reduction in the number of non-perturbative parameters. In particular, our formalism may shed some light on the power counting problems that appear in the heavy quarkonium polarization data [28].

APPENDIX A: 4-FERMION OPERATORS
Here we list the relevant 4-fermion operators of dimension 6 and 8, as taken from Ref. [1], where we use the conventional notation T (ij) ≡ (T ij +T ji )/2−T kk δ ij /3. The electromagnetic operators are defined as follows: where |vac is the vacuum state of QCD.

APPENDIX B: DIRECT MATCHING TO PNRQCD OF NRQCD MATRIX ELEMENTS
In principle, it is possible to match directly to pNRQCD matrix elements of NRQCD that involve operators different from the Hamiltonian H. In this way NRQCD matrix elements can be expressed in terms of non-local correlators without going through the full matching procedure outlined in the main body of the paper. This is useful if no iteration of these NRQCD operators are necessary in the matching calculation. In order to do this it is necessary to have an explicit expression for the state |0; x 1 , x 2 . Up to O(1/m) it can be found in Eq. (36). This way of doing will be particularly useful in order to work out higher order operators that will appear in going beyond O(mv 5 ) in the expansion of the heavy quarkonium decay width. Higher order operators appear to be necessary for charmonium decays, where the non-relativistic expansion converges slowly, assuming v 2 c ∼ 0.3. The master equation is (|H represents a generic heavy-quarkonium state at rest, P = 0, with quantum numbers n, j, l and s as defined in Ref. [1]): where r = x 1 − x 2 , r ′ = x ′ 1 − x ′ 2 , R = (x 1 + x 2 )/2 and R ′ = (x ′ 1 + x ′ 2 )/2 (note that R ′ |P = 0 = 1 and P = 0|P = 0 = d 3 x). As an example, let us consider here the NRQCD matrix element of the dimension-9 operator F EM ( 3 P 0 ) = 1 6 ψ † σ · gEχ|vac vac|χ † σ · Dψ + H.c., which is relevant to describe the electromagnetic decay χ c0 → γγ at order mv 7 accuracy. Owing to spin symmetry, the same matrix element enters into the χ c2 → γγ decay. These contributions have recently been considered in [39]. In the Hamiltonian formalism of Sec. III the matrix element (B2) is written as where |n011 is the Schrödinger wave function of the state χ Q (n01). Now we expand the state 0; x 1 x 2 | according to Sec. III C. The first non-vanishing contribution comes from the 1/m correction given in Eq. (36). Inserting that expression into Eq. (B4) and keeping in mind that only the term with the derivative projects onto the |n011 state, we obtain In the second equality we have made use of Eq. (53) and of the Wick theorem. Finally from the fact that δ (3) (r)|0 (0) = δ (3) (r)1 l c / √ N c |vac and from Eq. (73) we get or equivalently, using Eq. (151), Similar considerations may in principle also be applied to the matrix elements needed at relative order v 4 for S-wave decays. For a complete set of these operators and for considerations concerning their relevance in phenomenological studies, see Ref. [40].

APPENDIX C: RUNNING EQUATIONS OF THE MATCHING COEFFICIENTS
The running equations obtained in Appendix B.3 of Ref. [1] for the NRQCD 4-fermion operators give us information on the running of their matching coefficients. The running equations read as follows ν Im g 8 ( 1 S 0 )(ν) = Im g 8 ( 1 S 0 )(m) − Consider next: since, upon the translation p ′ → p ′ + p ′′ , the integral over p ′ has no scales. Alternatively, if we use a cut-off regularization, for instance by smoothing the delta in momentum space, like: we obtain: which can be removed by local counterterms. Hence DR implements nothing but a suitable subtraction prescription. Analogously, it is easy to see that δ (3) (r)/r s for s = 0, 1, 2, ... reduces to local terms.

APPENDIX E: UNITARY TRANSFORMATIONS
It is well known that quantum-mechanical Hamiltonians, which are related by unitary transformations, lead to the same physics. This fact is particularly relevant to quantummechanical Hamiltonians, which are derived from a field theory, which is our case. It is also well known that the quantum-mechanical potentials that are obtained from QED depend on the gauge one uses in the calculation (this is also so for QCD in perturbation theory), but physical observables computed with either potential turn out to be the same. It is perhaps not so well known that the potentials obtained in one gauge can be related to the ones obtained with a different gauge by means of a unitary transformation. In fact the arbitrariness in the form of the potentials is not only due to gauge dependence. It depends in general on the way one carries out the matching calculation. Any correct result is related to any other one by means of a unitary transformation.
We shall use this fact here to prove that the result obtained in Sec. V D is equivalent to the one obtained in Sec. IV A.
Consider the following unitary transformation: Consider next a delta function in the Hamiltonian: which on S-wave states reduces to:  (Fig. a)) and £ É Ñ ¾ (Fig. b)). All corrections not contained in ´¾ Øµ ¿ and ´¾µ ¿ arise from them after expanding the internal propagators. This generates all terms exhibited in Figure 3.  Figure 3: Diagrams generated by those in Fig. 2. a) corresponds to a È wave octet correction; b) and c) give rise to a chromomagnetic two-field correlator accompanying both a spin-flip/octet and a non-flip/singlet imaginary coefficient, respectively; d) produces the term proportional to ¿ times the binding energy; e) shows the structure introduced by the ÁÑ ½´¾ Ë·½ Ë Ë µ-proportional contact interaction. a) ¤ ❶ ❶ ❶ ❶ Figure 4: Diagrams corresponding to the next-to-leading contributions in the £É ÑÚ ¾ expansion. After expansion of the internal propagators, as explained in the text, they produce the series of graphs presented in figures 5 to 9, which originate the terms proportional to £ ❶ ❶ À ❶ ❶ ¤ ❶ ❶ À ❶ ❶ Ö ¾ Ñ ¤ AE ÁÑ ½´¾ Ë·½ Ë Ë µ Ñ ¾ AE´¿ µ´Ö µ À:= ´ × µ Figure 9: Diagrams contributing to the potential generated by the vacuum insertion in Fig. 4d). a) causes the structure ¼ to appear. The four of them are responsible for the ¿ ½ term. The operators and ¤ act through suitable commutations, which are not reflected in the figures, on the vertices. À must be taken left and right according to the prescription given in the text.