Semi-inclusive radiative decays of Upsilon 1S

We discuss in detail the photon spectrum of radiative Upsilon 1S decays taking into account a number of results that have recently appeared in the literature. In particular, we show how to consistently combine expressions which are valid in the upper end-point region, where NRQCD factorization breaks down, with those of the central region, where NRQCD factorization holds. An excellent description of data is achieved, but theoretical errors are large.


I. INTRODUCTION
Semi-inclusive radiative decays of heavy quarkonium systems (see [1] for a review) to light hadrons have been a subject of investigation since the early days of QCD [2,3]. It was thought for some time that a reliable extraction of α s was possible from the photon spectrum normalized, for instance, to the decay into muon pairs. However, the upper end-point region of the spectrum (namely z → 1, z being the fraction of the maximum energy the photon may have) appeared to be poorly described by the theory even when Sudakov resummations were carried out [4]. This led to some authors to claim that a non-vanishing gluon mass was necessary in order to describe the data [5], even when relativistic corrections were taken into account [6]. Later on, with the advent of Non-Relativistic QCD (NRQCD) [7], these decays could be analyzed in a framework where short distance effects, at the scale of the heavy quark mass m or larger, could be separated in a systematic manner [8]. These short distance effects are calculated perturbatively in α s (m) and encoded in matching coefficients whereas long distance effects are parameterized by matrix elements of local NRQCD operators. Even within this framework, a finite gluon mass seemed to be necessary to describe data [9]. However, about the same time it was pointed out that in the upper end-point region the NRQCD factorization approach breaks down and shape functions, namely matrix elements of non-local operators, rather than NRQCD matrix elements, must be introduced [10] . Early attempts to modeling color octet shape functions produced results in complete disagreement with data [11], and hence later authors did not include them in their phenomenological analysis [12,13]. Notwithstanding this region has received considerable attention lately, as it was recognized that the so called Soft-Collinear Effective Theory (SCET) [14,15] may help in organizing the calculation and in performing resummations of large (Sudakov) logs [12,13,16,17]. In fact, the early resummation of Sudakov logarithms [4] has been recently corrected [17] within this framework, and statements about the absence of Sudakov suppression in the color singlet channel [18] have been clarified [13]. For the Υ(1S) state, the bound state dynamics is amenable of a weak coupling analysis, at least as far as the soft scale (mv, v ∼ α s (mv) << 1, the typical velocity of the heavy quark in the quarkonium rest frame) is concerned [19,20,21,22,23,24]. These calculations can most conveniently be done in the framework of potential NRQCD (pNRQCD), a further effective theory where the contributions due to the soft and ultrasoft (∼ mv 2 ) scales are factorized [25,26,27] (see [28] for a review). Recently a calculation of the color octet shape functions, which combines SCET and pNRQCD, has become available [29].
It is the aim of this work to put together all known theoretical ingredients for these decays in order to see if a good description of data is achieved in the whole range of z, without the introduction of a finite gluon mass [30]. The theoretical calculation of the so called direct contributions is under good parametric control in the central region and in most of the upper end-point region. Indeed, if one uses the original NRQCD velocity counting rules [7], namely α s (m) ∼ v 2 and α s (mv) ∼ v, together with existing calculations at weak coupling, a complete NLO expression can be put forward in the central region. For the upper end-point region a complete LO expression, which includes color octet contributions and takes into account both Sudakov and Coulomb resummations, is also available. The merging of the central and the upper end-point regions will be discussed in detail. The fragmentation contributions, i. e. those for which the photon originates from the decay products of the heavy quarks, are parametrically of the same order as the direct photon contributions in the central region [31] and overweight the direct photon contributions in the lower end-point region (z → 0). They play a minor, but non-negligible, role in our analysis and are subject to large theoretical uncertainties.
We distribute the paper as follows. In the next section we separate the contributions to the decay width into direct and fragmentation. Sections III and IV are devoted to either contributions respectively. In section V we carry out the phenomenological analysis and section VI is devoted to the conclusions.

II. THE PHOTON SPECTRUM
The contributions to the decay width can be split into direct ( dir ) and fragmentation ( f rag ) We will call direct contributions to those in which the observed photon is emitted from the heavy quarks and fragmentation contributions to those in which it is emitted from the decay products (light quarks). This splitting is correct at the order we are working but should be refined at higher orders. z ∈ [0, 1] is defined as z = 2E γ /M (M is the mass of the heavy quarkonium state), namely the fraction of the maximum energy the photon may have in the heavy quarkonium rest frame.

III. DIRECT CONTRIBUTIONS
The starting point is the QCD formula [10] dΓ where J µ (x) is the electromagnetic current for heavy quarks in QCD and we have restricted ourselves to 3 S 1 states. q is the photon momentum, which in the rest frame of the heavy quarkonium is q = (q + , q − , q ⊥ ) = (zM, 0, 0). We have used light cone coordinates q ± = q 0 ± q 3 . The approximations required to calculate (2) are different in the lower end-point region (z → 0), in the central region (z ∼ 0.5) and in the upper end-point region (z → 1). We will denote Γ dir by Γ c and Γ e in the central and upper end-point regions respectively (the expressions for the lower end-point region will not be necessary).

A. The central region
For z away from the lower and upper end-points (0 and 1 respectively), no further scale is introduced beyond those inherent of the non-relativistic system. The integration of the scale m in the time ordered product of currents in (2) leads to local NRQCD operators with matching coefficients which depend on m and z. At leading order one obtains where and e Q is the charge of the heavy quark. The α s correction to this rate was calculated numerically in ref. [32]. The expression corresponding to (4) in pNRQCD is obtained at lowest order in any of the possible regimes by just making the substitution where ψ n0 (0) is the wave function at the origin. The final result coincides with the one of the early QCD calculations [2,3]. We will take the Coulomb form ψ 10 (0) = γ 3 /π for the LO analysis of Υ(1S) (γ is defined in (A13)). The NLO contribution in the original NRQCD counting [7] is v 2 suppressed with respect to (3). It reads In the original NRQCD counting or in the weak coupling regime of pNRQCD the new matrix element above can be written in terms of the original one [33] The matching coefficient can be extracted from an early calculation [6] (see also [35]). It reads where (ξ = 1 − z) The contributions of color octet operators start at order v 4 . Furthermore, away of the upper end-point region, the lowest order color octet contribution identically vanishes [8]. Hence there is no 1/α s enhancement in the central region and we can safely neglect these contributions here.
If we use the counting α s (µ h ) ∼ v 2 , α s (µ s ) ∼ v (µ h ∼ m and µ s ∼ mv are the hard and the soft scales respectively) for the Υ(1S), the complete result up to NLO (including v 2 suppressed contributions) can be written as The first term consist of the expression (3) with the Coulomb wave function at the origin (5) including corrections up to O (α s (µ s )) 2 [36,37], the second term is given in (6), and the third term consists of the radiative O (α s (µ h )) corrections to (3) which have been calculated numerically in [32]. Let us mention at this point that the O (α s (µ s )) 2 corrections to the wave function at the origin turn out to be as large as the leading order term. This will be important for our final results. Note that the standard NRQCD counting we use does not coincide with the usual counting of pNRQCD in weak coupling calculations, where α s (µ h ) ∼ α s (µ s ) ∼ α s (mv 2 ). The latter is necessary in order to get factorization scale independent results beyond NNLO for the spectrum and beyond NLO for creation and annihilation currents. However, for the Υ(1S) system (and the remaining heavy quarkonium states) the ultrasoft scale mv 2 is rather low, which suggests that perturbation theory should better be avoided at this scale [20]. This leads us to standard NRQCD counting. The factorization scale dependences that this counting induces can in principle be avoided using renormalization group techniques [38,39,40,41,42]. In practice, however, only partial NNLL results exists for the creation and annihilation currents [43,44] (see [45] for the complete NLL results), which would fix the scale dependence of the wave function at the origin at O(α 2 s (mv)). We will not use them and will just set the factorization scale to m.

B. The lower end-point region
For z → 0, the emitted low energy photon can only produce transitions within the non-relativistic bound state without destroying it. Hence the direct low energy photon emission takes place in two steps: (i) the photon is emitted (dominantly by dipole electric and magnetic transitions) and (ii) the remaining (off-shell) bound state is annihilated into light hadrons. It has a suppression ∼ z 3 with respect to Γ 0 (see [46,47] for a recent analysis of this region in QED). Hence, at some point the direct photon emission is overtaken by the fragmentation contributions [8,31].
In practice this happens about z ∼ 0.4, namely much before than the z 3 behavior of the low energy direct photon emission can be observed, and hence we shall neglect the latter in the following.

C. The upper end-point region
In this region the standard NRQCD factorization is not applicable [10]. This is due to the fact that small scales induced by the kinematics enter the problem and have an interplay with the bound state dynamics. In order to study this region, one has to take into account collinear degrees of freedom in addition to those of NRQCD. This can be done using SCET as it has been described in [13,16]. In this region, the color octet contributions are only suppressed by v 2 or by 1 − z. Since their matching coefficients are enhanced by 1/α s (µ h ), they become as important as the color singlet contributions if we count α s (µ h ) ∼ v 2 ∼ 1 − z. We will write where CS and CO stand for color singlet and color octet contributions respectively.

Color singlet contributions
For the color singlet contribution we shall use the expression with the Sudakov resummed coefficient in ref. [17] where the definitions for the different functions appearing in (12) are collected in the Appendix A.

Color octet contributions
For the color octet contributions we use 2 The matching coefficients provided in this reference become imaginary for extremely small values of z − 1, a region where our results do not hold anyway. We have just cut-off this region in the convolutions.
The (tree level) matching coefficients (up to a global factor) and the various shape functions are encoded in S S+P (x), The definitions of all the functions appearing in (14) and (15) are collected in the Appendix A. The shape functions S S , S P 1 and S P 2 may become S MS S , S MS P 1 and S MS P 2 or S sub S , S sub P 1 and S sub P 2 depending on the subtraction scheme employed. The procedure used to renormalize the shape functions is explained in the Appendix B. In figure 1 we plot the end-point contribution (11) with the shape functions renormalized in an M S scheme (blue dashed line) and in the sub scheme, which makes additional subtractions (red solid line), together with the experimental data [48] (we have convoluted the theoretical curves with the experimental efficiency, the overall normalization of each curve is taken as a free parameter. For the details of the scale setting see section V). We see that both schemes are equally good for the description of the shape of the experimental data in the end-point region.
Note (from Appendix A) that we use the octet shape functions calculated in ref. [29], which are crucial in order to have a good description of data in the upper end-point region. We would like to comment on the validity of those formulas. This is limited by the perturbative treatment of the US gluons. The typical momentum of these gluons in light cone coordinates turns out to be: Note that the typical k ⊥ is not fixed by the bound state dynamics only but by a combination of the latter and the end-point kinematics. Hence, the calculation is reliable provided that k ⊥ 1GeV., which means z < 0.92. Note also that the typical three momentum of the heavy quarks in the shape function is given by This means that the multipole expansion always holds in the end-point region, regardless that M (1 − z) is bigger or smaller than mα 2 s . This important point was not sufficiently emphasized in [29] and it is the ultimate reason why such a good description of data was obtained there for z ∈ [0. 7,1]. Recall that at z < 0.7 or so we are at the border of the end-point region and the standard NRQCD factorization formulas should hold. This means that the contributions of the octet shape functions should merge suitable contributions of the NRQCD factorization formulas, as we discuss below. One way to proceed is the following. If we assume that the expressions for the end-point contain the ones of the central region up to a certain order in (1 − z), we could just subtract from the expressions in the central region the behavior when z → 1 at the desired order and add the expressions in the end-point region. Indeed, when z → 1 this procedure would improve on the central region expressions up to a given order in (1 − z), and when z belongs to the central region, they would reduce to the central region expressions up to higher orders in α s . This method was used in ref. [13] and in ref. [49]. In ref. [13] only color singlet contributions were considered and the end-point expressions trivially contained the central region expressions in the limit z → 1. In ref. [49] color octet contributions were included, which contain terms proportional to (1 − z). Hence, the following formula was used Even though a remarkable description of data was achieved with this formula (upon using a suitable subtraction scheme described below), this method suffers from the following shortcoming. The hypothesis that the expressions  [48] for the end-point contain the ones for the central region up to a given order in (1 − z) is in general not fulfilled. As we will see below, typically, they only contain part of the expressions for the central region. This is due to the fact that some α s (µ h ) in the central region may soften as α s (M (1 − z)), others as α s (M √ 1 − z) and others may stay at α s (µ h ) when approaching the end-point region. In a LO approximation at the end-point region, only the terms with the α s at low scales would be kept and the rest neglected, producing the above mentioned mismatch. We shall not pursue this procedure any further.
Let us look for an alternative. Recall first that the expressions we have obtained for the upper end-point region are non-trivial functions of M (1 − z), M √ 1 − z, mα s (mv) and mα 2 s (mv), which involve α s at all these scales. They take into account both Sudakov and Coulomb resummations. When z approaches the central region, we can expand them in α s (M √ 1 − z), α s (M (1 − z)) and the ratio mα s (mv)/M √ 1 − z. They should reduce to the form of the expressions for the central region, since we are just undoing the Sudakov and (part of) the Coulomb resummations. Indeed, we obtain where . The details of this derivation are given in the Appendix C. The color singlet contribution reproduces the full LO expression for the central region in the limit z → 1. The color octet shape functions S P 1 and S P 2 give contributions to the relativistic corrections (6), and S P 2 to terms proportional to (1 − z) in the limit z → 1 of (3) as well. We have checked that, in the z → 1 limit, both the (1 − z) ln(1 − z) of (3) and the ln(1 − z) of the relativistic correction (6) are correctly reproduced if µ c ∼ M √ 1 − z, as it should. All the color octet shape functions contribute to the O(α s (µ h )) correction in the first line of (20). There are additional O(α s (µ h )) contributions coming from the expansion of the (Sudakov) resummed matching coefficients of the color singlet contribution and of the S P 2 color octet shape function. The α s log(1 − z) in (19) reproduces the logarithm in dΓ c LO,αs /dz.
We propose the following formula This formula reduces to the NRQCD expression in the central region. When we approach the upper end-point region the second terms in each of the parentheses are expected to cancel corresponding terms in the z → 1 limit of the expression for the central region up to higher order terms (in the end-point region counting). Thus, we are left with the resummed expressions for the end-point (up to higher order terms).
There are of course other possibilities for the merging. For instance, one may choose a z 1 below which one trusts the calculation for the central region and a z 2 above which one trusts the end-point region calculation, and use some sort of interpolation between z 1 and z 2 (see for instance [50]). This would have the advantage of keeping the right approximation below z 1 and beyond z 2 unpolluted, at the expense of introducing further theoretical ambiguities due to the choice of z 1 and z 2 , and, more important, due to the choice of the interpolation between z 1 and z 2 . We believe that our formula (22) is superior because it does not introduce the above mentioned theoretical ambiguities. The price to be paid is that the expressions from the central region have an influence in the end-point region and vice-versa. This influence can always be chosen to be parametrically subleading but large numerical factors may make it noticeable in some cases, as we shall see below.

Merging at LO
If we wish to use only the LO expressions for the central region, we should take (20) at LO, namely and substitute them in (22). Unexpectedly, the results obtained with this formula in the central region deviate considerably from those obtained with formula (3) (see Fig. 2). This can be traced back to the fact that the α s √ 1 − z corrections in (20) are enhanced by large numerical factors, which indicates that the merging should better be done including α s (µ h ) corrections in the central region, as we discuss in the next section. Alternatively, we may change our subtraction scheme in order to (partially) get rid of these contributions, as discussed in the Appendix B. With the new subtraction scheme (sub) the situation improves, although it does not become fully satisfactory (see Fig. 2). This is due to the fact that some α s √ 1 − z terms remain, which do not seem to be associated to the freedom of choosing a particular subtraction scheme. In spite of this the description of data turns out to be extremely good. In figure 3 we plot, using the sub scheme, the merging at LO (solid red line) and also, for comparison, equation (18) (blue dashed line). We have convoluted the theoretical curves with the experimental efficiency and the overall normalization is taken as a free parameter.

Merging at NLO
If we wish to use the NLO expressions for the central region (10), we should take all the terms displayed in (19)- (20) and substitute them in (22). Unlike in the LO case, for values of z in the central region the curve obtained from are an accident and that the α 3 s (µ s ) ones (see [51,52] for partial results) will again be small, one should use these α 2 s (µ s ) corrections. We consider below the two cases.
If we stay at LO (or including α s (µ s ) corrections) for the wave function at the origin, the curve we obtain for z → 1 differs considerably from the expressions for the end-point region (11) (see Fig. 4). This can be traced back to the α s √ 1 − z term in (20) again. This term is parametrically suppressed in the end-point region, but, since it is multiplied by a large numerical factor, its contribution turns out to be overwhelming. This term might (largely) cancel out against higher order contributions in the end-point region, in particular against certain parts of the NLO expressions for the color singlet contributions, which are unknown at the moment.
If we use the wave function at the origin with the α 2 s (µ s ) corrections included, the curves we obtain for z → 1 become much closer to the expressions for the end-point region (11) (see Fig 5). Hence, a good description of data is obtained with no need of additional subtractions, as shown in figure 6 (as usual experimental efficiency has been taken into account and the overall normalization is a free parameter).

IV. FRAGMENTATION CONTRIBUTIONS
The fragmentation contributions can be written as where C a represents the partonic kernels and D aγ represents the fragmentation functions. The partonic kernels can again be expanded in powers of v [8] The leading order term in v is the color singlet rate to produce three gluons  The color octet contributions start at order v 4 but have a 1 αs enhancement with respect to (26) Then the color singlet fragmentation contribution is of order α 3 s D g→γ and the color octet fragmentation are of order v 4 α 2 s D g→γ ( 1 S 0 and 3 P J contributions) or v 4 α 2 s D q→γ ( 3 S 1 contribution). We can use, as before, the counting v 2 ∼ α s to compare the relative importance of the different contributions together with the existing models for the fragmentation functions [53]. The latter tell us that D q→γ is much larger than D g→γ . This causes the O(v 4 α 2 s D q→γ ) 3 S 1 octet contribution to dominate in front of the singlet O(α 3 s D g→γ ) and the octet O(v 4 α 2 s D g→γ ) contributions. In fact, α s D q→γ is still larger than D g→γ , so we will include in our plots the α s corrections to the color octet contributions (27) proportional to D q→γ , which have been calculated in [8]. In addition, the coefficients for the octet 3 P J contributions have large numerical factors, causing these terms to be more important than the color singlet contributions. Let us finally notice that the α s corrections to the singlet rate will produce terms of O(α 4 s D q→γ ), which from the considerations above are expected to be as important as the octet 3 S 1 contribution. These α s corrections to the singlet rate are unknown, which results in a large theoretical uncertainty in the fragmentation contributions.
For the quark fragmentation function we will use the LEP measurement [54] where and for the gluon fragmentation function the model [55]. These are the same choices as in [13]. However, for the O 8 ( 1 S 0 ) and O 8 ( 3 P 0 ) matrix elements we will use our estimates in [29] Υ(1S) The above numbers are obtained in an M S scheme from dimensionally regularized US loops only. The value we assign to the S-wave matrix element is compatible with the recent (quenched) lattice determination (hybrid algorithm) [56]. Notice that we do not assume that a suitable combination of these matrix elements is small, as it was done in [13]. The O 8 ( 3 S 1 ) matrix element can be extracted from a lattice determination of the reference above [56]. Using the wave function at the origin with the α 2 s (µ s ) corrections included, we obtain, which differs from the estimate using NRQCD v scaling by more than two orders of magnitude: (we have taken v 2 ∼ 0.08), which was used in ref. [13]. The description of data turns out to be better with the estimate (33). However, this is not very significant, since, as mentioned before, unknown NLO contributions are expected to be sizable.
In the z → 0 region soft radiation becomes dominant and the fragmentation contributions completely dominate the spectrum in contrast with the direct contributions [31]. Note that, since the fragmentation contributions have an associated bremsstrahlung spectrum, they can not be safely integrated down to z = 0; that is 1 0 dz dΓ f rag dz is not an infrared safe observable. In any case we are not interested in regularizing such divergence because the resolution of the detector works as a physical cut-off.

V. SCALE SETTING AND ERROR ANALYSIS
Formula (22) requires dΓ e /dz for all values of z. The color octet shape functions, however, were calculated in the end-point region under the assumption that M √ 1 − z ∼ γ, and the scale of the α s was set accordingly. When z approaches the central region M √ 1 − z ≫ γ, and hence some α s will depend on the scale M √ 1 − z and others on γ (we leave aside the global α s (µ u ), which will be discussed below). In order to decide the scale we set for each α s let us have a closer look at the formula (20). We see that all terms have a common factor γ 3 . This indicates that one should extract γ 3 factors in the shape functions, the α s of which should stay at the scale µ s . This is achieved by extracting γ 3/2 in I S and I P . If we set the remaining α s to the scale µ p = m(M (1 − z)/2 − E 1 ), we will reproduce (20) when approaching to the central region, except for the relativistic correction, the α s of which will be at the scale µ p instead of at the right scale µ s . We correct for this by making the following substitution Notice that the replacements above are irrelevant as far as the end-point region is concerned, but important for the shape functions to actually (numerically) approach the expressions (20)  Our final plot in Fig. 7 is obtained by using the merging formula (22) at NLO with the α 2 s (µ s ) corrections to the wave function at the origin included for the direct contributions plus the fragmentation contributions in section IV including the first α s corrections in C q and using the estimate (33) for the Υ(1S)|O 8 ( 3 S 1 )|Υ(1S) matrix element. The error band is obtained by replacing µ c by √ 2 ±1 µ c . Errors associated to the large α 2 s (µ s ) corrections to the wave function at the origin, to possible large NLO color singlet contributions in the end-point region and to the fragmentation contributions are difficult to estimate and not displayed (see the corresponding sections in the text for discussions). The remaining error sources are negligible. As usual experimental efficiency has been taken into account and the overall normalization is a free parameter.

VI. CONCLUSIONS
We have analyzed the photon spectrum in radiative Υ(1S) decays within an Effective Field Theory framework. For the direct contributions, the merging of the results for the central and upper end-point regions has been discussed in detail. We have shown how to consistently combine the complete LO results for the upper end-point region with the complete NLO ones for the central region. We have seen that the large α 2 s (µ s ) corrections to the wave function at the origin are important in order to get a good description of data. Otherwise, parametrically subleading large contributions in the end-point region would be necessary. We would like to emphasize that our final results for the direct contributions are essentially parameter free: only the mass of the bottom quark m, the strong coupling constant α s , and the proper choice of subtraction scales (which appear in logarithms) are used as an input. For the fragmentation contributions, we have pointed out that if the commonly used model for the gluon fragmentation into a photon is appropriated, α s corrections to the LO color singlet matrix element giving rise to a light quark which fragments into a photon may be as important as the LO results. Hence, fragmentation contributions suffer from large theoretical uncertainties. Nevertheless, if we put together the available theoretical results for these contributions with the ones for the direct contributions, an excellent description of data is achieved for the whole part of the spectrum where experimental errors are reasonable small. Clearly, our results indicate that the introduction of a finite gluon mass [30] is unnecessary. One should keep in mind, however, that in order to have the theoretical errors under control higher order calculations are necessary both in the direct (end-point) and fragmentation contributions.
Before closing, let us mention that the inclusion of color octet contributions in the end-point region together with the merging with the central region expression described in this work may be useful for production processes like inclusive J/ψ production in e + e − machines [50,58,59].  (33) and (32) for Υ(1S)|O8( 3 S1)|Υ(1S) respectively. The grey shaded region is obtained by varying µc by √ 2 ±1 µc. The green shaded region on the right shows the zone where the calculation of the shape functions is not reliable (see section III C 2). The pink dashed line is the result in [13], where only color singlet contributions were included in the direct contributions.

Acknowledgments
We are grateful to Dave Besson and to Michael Kramer for providing us with the data of ref. [48] and the numerical results of ref. [32] respectively. Thanks are also given to Antonio Pineda for important discussions, and to Antonio Vairo for asking the right questions. We acknowledge financial support from a CICYT-INFN 2004 collaboration contract, Acciones Integradas HI2003-0362 (Spain-Italy), the MCyT and Feder (Spain) grant FPA2001-3598, the CIRIT (Catalonia) grant 2001SGR-00065 and the network EURIDICE (EU) HPRN-CT2002-00311. X.G.T. thanks the Institute for Nuclear Theory at the University of Washington for its hospitality and the Department of Energy for partial support during the completion of this work. X.G.T. also acknowledges financial support from the Departament d'Universitats, Recerca i Societat de la Informació of the Generalitat de Catalunya and the Fons Social Europeu.

APPENDIX A: DEFINITIONS
In this appendix we collect the definitions for the formulas that appear in the paper.
is the number of colors, n f = 4 is the number of light flavors.

APPENDIX B: REGULARIZATION AND RENORMALIZATION
The shape functions (A10)-(A12) are ultraviolet (UV) divergent and require regularization and renormalization. In ref. [29] it was pointed out that using dimensional regularization (DR) for the US loop only was enough to regulate them. In fact, the expressions (A10)-(A11) implicitly assume that DR is used, otherwise linearly divergent terms proportional to ψ 2 10 (0) would appear (which make (A10)-(A11) formally positive definite quantities). In addition an M S scheme was used to subtract the poles. Since it turns out that the final outcome strongly depends on the details of this subtraction, let us spell out the procedure carried out in ref. [29]. In order to isolate the 1/ε poles, I S and I P were expanded up to O(1/z ′ 2 ) . The result was subtracted and added to the integrand of (A10)-(A11) (for (A12) this is not necessary since the only divergent piece is independent of I P ). The subtracted part makes the shape functions finite. The added part contains linear and logarithmic UV divergencies. The 1/ε (D = 4 − 2ǫ) poles displayed in formulas (16) of ref. [29], and eventually subtracted, were obtained by making dx → dx(x/µ) −ε in (A10)-(A12). This was motivated by the fact that x ∼ k 2 ⊥ (k ⊥ being the transverse momentum of the US gluon) but differs from a standard M S scheme. Linear divergences are set to zero as usual in DR.
We have used here a regularization and renormalization scheme which is closer to the standard one in pNRQCD calculations. We have regulated both the US loop and the potential loops (entering in the bound state dynamics) in DR. We have identified US divergencies by taking the limit D → 4 in the US loops while leaving the potential loops in D dimensions [60]. Potential divergencies are identified by taken D → 4 in the potential loops once the US divergencies have been subtracted. It turns out that all divergencies in S P 2 are US and all divergencies in S S are potential. S P 1 contains both US and potential divergencies. The potential divergences related with the bound state dynamics can be isolated using the methods of ref.
[61]. The formulas corresponding to (16) of ref. [29] in this scheme read For simplicity, we have set D = 4 everywhere except in the momentum integrals. µ p is defined in section V. µ c and µ pc are the subtraction points of the US and potential divergencies respectively. If we subtract the 1/ε poles and set µ c = M √ 1 − z and µ pc = √ mµ c we obtain exactly the same result as in ref. [29] for what the potential divergences is concerned 3 . For the US divergences there is a factor ln µc 2k+ of difference with respect to the previous scheme.
In ref. [49] an additional subtraction related to linear divergencies was made. This subtraction was necessary in order to merge smoothly with the results in the central region. We will also need this subtraction here when merging at LO, as discussed in subsection III D. We use , which differ from the M S scheme by the subtraction of the second term in the square brackets.

APPENDIX C: THE SHAPE FUNCTIONS IN THE CENTRAL REGION
When z approaches the central region from the upper end-point, the shape functions should reduce to matrix elements of NRQCD operators multiplied by the corresponding matching coefficients. We will see here that this is indeed the case.
Let us first consider the S-wave octet shape function as defined in [29] I S ( k + 2 + x) := d 3 xψ 10 (x) 1 − 2N c |r|). When z approaches the central region, k + ∼ M (1 − z) ≫ −E 1 and the larger three momentum scale is M √ 1 − z ≫ γ, the typical three momentum in the bound state. Therefore we can treat the Coulomb potential in (C1) as a perturbation when it is dominated by this scale. It is convenient to proceed in two steps. First we write h o = h s + (V o − V s ), where h s = p 2 /m + V s , V s = −α s C f /|r|, and expand V o − V s . This allows to set h s − E 1 to zero in the left-most propagator and makes explicit the cancellation between the first term in the series and the first term in (C1). It also makes explicit that the leading term will be proportional to α s (M √ 1 − z). Second, we expand V s in h s = p 2 /m + V s . In addition, since M √ 1 − z ≫ γ, the wave function can be expanded about the origin. Only the first term in both expansion is relevant in order to get (20).
Consider next the P -wave shape functions as defined in [29] I P ( In order to proceed analogously to the S-wave case, we have first to move the x i away from the wave function For the left-most propagators we can now proceed as before, namely expanding V o − V s . Note that the leading contribution in this expansion of the second term above exactly cancels against the first term. Of the remaining contributions of the second term only the next-to-leading one (O(α s )) is relevant to obtain (20). Consider next the leading order contribution in this expansion of the last term. It reads Now we proceed as before with the left-most propagators, namely expanding V o − V s . The leading order contribution of the first term above produces the relativistic correction O(v 2 ) of (20). The next-to-leading contribution of this term and the leading order one of the second term are O(α s ) and also relevant to (20). The next-to-leading order contribution of the last term in (C3) in the V o − V s expansion of the left-most propagator is also O(α s ) and relevant to (20).