Power expansion for heavy quarkonium production at next-to-leading order in $\rm e^+e^-$ annihilation

We study heavy quarkonium production associated with gluons in $\rm e^+e^-$ annihilation as an illustration of the perturbative QCD (pQCD) factorization approach, which incorporates the first nonleading power in the energy of the produced heavy quark pair. We show how the renormalization of the four-quark operators that define the heavy quark pair fragmentation functions using dimensional regularization induces"evanescent"operators that are absent in four dimensions. We derive closed forms for short-distance coefficients for quark pair production to next-to-leading order ($\alpha_s^2$) in the relevant color singlet and octet channels. Using non-relativistic QCD (NRQCD) to calculate the heavy quark pair fragmentation functions up to $v^4$ in the velocity expansion, we derive analytical results for the differential energy fraction distribution of the heavy quarkonium. Calculations for ${}^3S_1^{[1]}$ and ${}^1S_0^{[8]}$ channels agree with analogous NRQCD analytical results available in the literature, while several color-octet calculations of energy fraction distributions are new. We show that the remaining corrections due to the heavy quark mass fall off rapidly in the energy of the produced state. To explore the importance of evolution at energies much larger than the mass of the heavy quark, we solve the renormalization group equation perturbatively to two-loop order for the ${}^3S_1^{[1]}$ case.

Abstract: We study heavy quarkonium production associated with gluons in e + e − annihilation as an illustration of the perturbative QCD (pQCD) factorization approach, which incorporates the first nonleading power in the energy of the produced heavy quark pair. We show how the renormalization of the four-quark operators that define the heavy quark pair fragmentation functions using dimensional regularization induces "evanescent" operators that are absent in four dimensions. We derive closed forms for short-distance coefficients for quark pair production to next-to-leading order (α 2 s ) in the relevant color singlet and octet channels. Using non-relativistic QCD (NRQCD) to calculate the heavy quark pair fragmentation functions up to v 4 in the velocity expansion, we derive analytical results for the differential energy fraction distribution of the heavy quarkonium. Calculations for 3 S [1] 1 and 1 S [8] 0 channels agree with analogous NRQCD analytical results available in the literature, while several color-octet calculations of energy fraction distributions are new. We show that the remaining corrections due to the heavy quark mass fall off rapidly in the energy of the produced state. To explore the importance of evolution at energies much larger than the mass of the heavy quark, we solve the renormalization group equation perturbatively to two-loop order for the 3 S

Introduction
Heavy quarkonium production is a subject of continuing interest [1][2][3]. The production of a heavy quarkonium state always involves an intrinsic hard scale, the heavy quark mass, m Q . In the presence of an even larger hard scale, E H , such as transverse momentum or particle energy, large logarithms ln(E 2 H /(2m Q ) 2 ) can appear, which are fully perturbative but require resummation. The perturbative QCD (pQCD) factorization procedure developed in refs. [4,5] is an organization of high-energy quarkonium production into subprocesses of different characteristic regions in momentum and coordinate space to make this resummation possible. Closely-related results have been derived from an effective theory viewpoint in refs. [6,7]. Recently, the formalism also has been applied to light meson production in deep-inelastic scattering [8].
For quarkonium production at large transverse momentum, this pQCD factorization approach provides a unified framework for leading-power (LP) and next-to-leading power (NLP) behavior in momentum transfer for production cross sections. Suppressing convolutions associated with the initial state, such factorized cross sections for the production of heavy quarkonium H can be represented as (1.1) The convolutions in this expression are in hadronic momentum fraction z, and variables, u and v, which represent the fraction of the hadron's total momentum carried by the heavy quark in the amplitude and in its complex conjugate, respectively 1 . The first term represents the LP production of a single parton i = q,q, g at short distances. The second term describes the NLP production of a heavy quark pair in a specific spin and color state κ at short distances. The full expression for the factorized cross section in eq. (1.1) with explicit convolutions is [4] E P dσ A+B→H+X d 3 P (P ) = i=q,q,g with the same plus component, wheren µ is the light-cone vectorn µ = (1, 0, 0 ⊥ ). The p + component can be projected out by the oppositely directed light-cone vector n µ = (0, 1, 0 ⊥ ), as p + = p · n. As indicated in eq. (1.2), the short-distance coefficients dσ depend only on the lightlike momenta p Q and pQ, related to p µ by p = z(p Q + pQ) = zp c . That is, the heavy quark is treated as massless in the short-distance coefficients. A proof of this factorization at NLP was provided in ref. [4]. It should be noted that, rather than the full set of NLP non-perturbative functions, the formalism includes only the heavy quark pair fragmentation functions. This is justified under the reasonable assumption that they give the dominant contributions at NLP when a heavy quarkonium is produced. As is characteristic of factorized power corrections [9][10][11][12][13][14][15][16], the fragmentation functions in this expression involve not only the momentum fraction carried by the pair, but also the relative momentum of the constituents of the pair, the heavy quark and antiquark. The short-distance functions, therefore, produce non-diagonal products of partonic states, in which the heavy quark pair has the same total momentum, p c in eq. (1.2), shared differently by the heavy quark and antiquark. The fragmentation function mediates between these products of states, non-diagonal in u and v, and the full set of states that include the heavy quarkonium, of momentum P . Operator definitions for the fragmentation functions can be found in ref. [4], given as Fourier transforms of matrix elements of four-quark operators, here suppressing gauge links, ij,lk (p) C [I] ab,cd X 0|ψ c,l (y − 1 )ψ d,k (0)|H(P )X H(P )X|ψ a,i (y − )ψ b,j (y − + y − 2 )|0 , (1.4) where C (I) and P (s) are members of a set of projections in color and Dirac spin space, respectively. Explicit forms of projections will be given in section 2. We will apply this formalism to heavy quarkonium production in association with gluons for e + e − annihilation. For these processes, we provide the first closed expressions for the NLP short-distance coefficients at next-to-leading order (NLO). We will encounter a number of features that will recur in all the processes that can be factorized in this manner, and will thus play a role in any program to develop a phenomenology of heavy quark production that includes NLP mechanisms. In addition, once e + e − annihilation short-distance coefficients are calculated, the asymptotic behaviors of exact calculations in NRQCD at O(α 2 α 2 s ) allow us to check these short-distance coefficients explicitly, using the formalism for NRQCD fragmentation functions developed in refs. [17,18]. Previously, the two formalisms were compared numerically for hadronic scattering in ref. [19], with good results, but there is nothing like an analytical comparison.
The calculation of short-distance coefficients depends, of course, on the scheme used to define the matrix elements in eq. (1.4), which involve four quark fields relatively on the light cone. As in the case of single-parton fragmentation functions, such operator configurations require renormalization as "cut vertices" [20], as in ref. [4], and thus depend on the renormalization scheme. In this paper, we will use a modified minimal subtraction (MS) approach in dimensional regularization. The presence of the four-quark operator then makes it necessary to include the effects of operators that are absent in a purely four-dimensional calculation. This phenomenon is of great importance in the treatment of amplitudes mediated by effective four-quark operators [21,22]. In our case, we will encounter "evanescent" fragmentation functions, associated with projections P (s) that disappear in four dimensions. We will see that the comparison to exact NRQCD calculations, even at NLO, is ambiguous without a clear understanding of how these new functions enter into the calculation of factorized cross sections.
In this work, we study the energy fraction, x H = 2E H /Q, distributions of heavy quarkonium H produced in gluon-associated processes, which only contributes at NLP up to O(α 2 α 2 s ). Note that the Belle experiment is able to separate heavy quarkonium production associated with and without cc, σ(e + e − → J/ψ + cc) and σ(e + e − → J/ψ + X cc ) [23], and that the contributions without cc are mainly gluon-associated production at Belle's energy within the NRQCD formalism [24][25][26][27][28][29]. The LO color octet cross section begins at order O(α 2 α s ), where it would appear with an end-point enhancement [30] at x H ∼ 1. This, however, has not been observed in the data [24][25][26][27][28][29]. We focus our study in the region x H away from 1 and thus only calculate the real contributions of the NLO shortdistance coefficients. For NRQCD calculations of the gluon-associated processes, there exist NLO numerical calculations for both color-singlet channels and color-octet channels [26][27][28]. Closed expressions for the x H differential distributions are also available for 3 S [1] 1 [29,31,32] and 1 S [8] 0 [33] to order O(α 2 α 2 s ). We will use these results to make explicit comparison to the results of pQCD factorization for these cases and see that the fixed order expression of pQCD factorization can approach the high-energy limit of NRQCD even at moderate energies.
Perturbative QCD factorization can also go beyond NRQCD fixed order calculations by resumming logarithms like ln(E 2 H /(2m Q ) 2 ) through evolution equations [4]. We will solve the evolution equations to two loops to test the significance of evolution for these functions.
We begin in section 2 with a discussion of the factorization formalism of ref. [4], identifying additional evanescent fragmentation functions associated with dimensional regularization. We calculate the relevant short-distance coefficients for heavy quark pair production in association with gluons at NLO in section 3. In section 4, using the fragmentation functions of refs. [17,18], we compare our results to the asymptotic behavior of NRQCD calculations in the channels for which the O(α 2 α 2 s ) results are available and find agreement. Lastly, in section 5, we study the approach to asymptotic behavior and evolution for the 3 S 1 channel. We also show NLO pQCD predictions for different NRQCD channels at various center-of-mass (CM) energies, and include comparison to data from the Belle collaboration [23].

NLP factorization and fragmentation functions
As shown for LP in ref. [34] and NLP in ref. [4], the perturbative short-distance coefficients of the factorized single-particle inclusive cross section, eq. (1.2) are sensitive only to the large lightlike momentum components of the parton(s) whose fragmentation produces the observed hadron. For NLP involving a heavy quark pair, the corresponding contributions to the cross section are convolutions in the longitudinal momenta of the heavy quark and antiquark of the pair between the short-distance coefficients and the fragmentation functions. The short-distance coefficients and fragmentation functions describe the production of the pair, and its subsequent evolution into a heavy quarkonium, respectively. This NLP contribution is only one of many, but it is reasonable to assume that it dominates the class of NLP corrections simply because of the presence of the pair in heavy quarkonium production.

Pinch surfaces and factorization
In this subsection, we review the underlying structure that allows us to separate and normalize short-distance coefficients and fragmentation functions for NLP heavy quark pair production. This requires the construction of the heavy quark pair fragmentation functions at the partonic level, which begins with diagrammatic analysis.
The observation that underlies factorization procedures is that long-distance contributions to cross sections for hard-scattering processes arise from well-defined regions in phase space and loop momenta. These regions are associated with so-called "pinch surfaces" [35][36][37] in the massless limit. Each such region is associated with a specific power behavior in the high energy limit, LP, NLP and so forth. The LP and NLP contributions in particular can be factorized into perturbatively-computable hard parts (short-distance coefficients), which determine the power behavior, and universal parton distributions for incoming hadrons and fragmentation functions for observed hadrons, which absorb mass and other long-distance dependence. At large momentum transfer and to NLP, even the dependence on a heavy quark mass factorizes from the short-distance scattering [38]. In this approximation, LP and NLP contributions can be significantly simpler to determine than by derivation from full fixed-order calculations.
We recall first the procedure at LP. A generic pinch surface contributing to LP is illustrated in figure 1a for gluon fragmentation. The momentum of the gluon passing the horizontal dashed line in the figure represents a propagator near its mass shell. In the region in question, the subdiagram below the cut describes the short-distance production of a nearly on-shell gluon from a hard scattering. Subsequent fragmentation of the nearly on-shell gluon into a heavy quarkonium over long times is represented by the subdiagram above the line. Every pinch surface that contributes to the cross section at LP takes this form (for a gluon or quark), with a short-distance piece, which is perturbatively calculable, and a non-perturbative but universal function. This is what makes factorization possible.
Given this factorization, it is possible to improve systematically the calculation of the shortdistance coefficients, and as a result to determine systematically universal non-perturbative fragmentation functions by comparison to experiment or other non-perturbative input. This program relies on the assumption that the calculation of the short-distance coefficients is by construction independent of long-time dynamics, and can be carried out in an infraredregulated version of QCD. In particular, for LP fragmentation, we normally carry out the computation of the short-distance coefficient in dimensionally-regularized QCD, taking the observed final-state hadron to be a parton itself. from higher order calculations, even at high p T , suggest that we need to supplement the simplest perturbative expansion in the powers of α s for the CSM, and by implication for the NRQCD factorization of Eq. (1), to take into account radiation from heavy quark pairs produced at short, and intermediate, time scales.
We propose to expand the production cross section of heavy quarkonia at high p T in powers of 1/p T first, when p T ≫ m H , and only then to expand perturbatively factorizable hard parts in powers of α s . In the remainder of this section, we argue that the cross section for producing a heavy quarkonium at large transverse momentum at collider energies can be expanded as a power series of 1/p 2 T , and that the leading power term and the first subleading power terms can be perturbatively factorized into infrared safe short-distance hard parts in convolution with nonperturbative but universal long-distance fragmentation functions [25,26]. Schematically, all contributions to all pinch surfaces at LP for the production of a hadron at high energy can be represented as where we recall from eq. (1.3) that p µ is the light-like projection of hadron momentum P µ . The function H e + e − →aX represents all diagrammatic contributions to the production of parton a, below the dashed line in figure 1a, and T a→H(P ) those above. Here for simplicity of notation, we take 0 < x H < 1 as the fractional momentum of the observed hadron H in e + e − annihilation, integrated over phase space dΠ(P ). The starting form for our discussion for the contributions from generic NLP pinch surfaces involving the production of a heavy quark pair at short distances [4], illustrated by figure 1b, is × T ab,cd;ij,kl (z, u, v; P ) dΠ(P ) . (2.2) In the second expression we exhibit sums over two pairs of color indices (ab, cd) and two pairs of spin indices (ij, kl) linking the long distance part T and the short-distance part H in both the amplitude (ab, ij) and complex conjugate (cd, kl). The form of eq. (2.2) is already suggestive of the factorized NLP cross section in eq. (1.2). Bridging the gap between the two requires the identification of the dominant momentum dependence and the separation of color and spin traces between the long-and short-distance coefficients. These steps are described in ref. [4], and we will recall them here for their detailed implementation at NLO. Consider first the treatment of momentum flow. As in eq. (2.1), the limit on the z integral in eq. (2.2) refers to the production of a hadron integrated over phase space at fixed fractional momentum x H , while the u and v integrals go from zero to one, representing the fraction of the heavy pair's momentum carried by the heavy quark in the amplitude and complex conjugate, respectively.
The hard functions H depend only on the light-like projections of the heavy quark and antiquark momenta on either side of the (vertical) final state cut. Specifically, in terms of the parameters u . . .v introduced above, these momenta are given in terms of p µ from eq. (1.3) the light-cone projection of the final state momentum, P µ , by That is, in addition to the pair spin-color state κ, we must specify the momentum fractions u and v carried by the quarks in the amplitude and complex conjugate in the partonic final state. The factorized cross section will be a triple convolutions in z, u, and v. We now turn to the separation of the sums over color and spin indices. Here, we follow the method of ref. [4], modified for spin projections to accommodate dimensional regularization. Quite generally, at fixed values of u and v, we can expand the perturbative longdistance and short-distance functions into color singlet and octet matrices (components ab) times Dirac matrices (components ij). For example, for the short-distance function, H, we find [4] H e + e − →[QQ]X (p Q , pQ, p Q , p Q ) ab,ij;cd,kl =δ ab δ cd where H a,I are coefficients of the generators of spin state I and color state a. The two Dirac matrices correspond to the amplitude and complex conjugate side of the short-distance coefficients and they are diagonal [4] for unpolarized initial and final states. The t A are SU(3) generators in fundamental representation. In four dimensions, the Γ I can be taken from the usual basis of Dirac matrices in four dimensions: , γ µ , σ µν , γ µ γ 5 , γ 5 } . (2.5) We will use the structure of eq. (2.4) as a guide in the construction of factorizing projections in color and spin.

Color projections
The long-and short-distance parts in eq. (2.2) are connected by two heavy quark pairs, each of which can be in a singlet or octet state. The projection in eq. (2.2) is standard, and the process may be represented explicitly by the relations a b ,c d + C [8] ab,cdC Here we take the color projections C for the heavy quark pairs that enter T and isolate ones that give singlet or octet pair production. Gauge invariance ensures that these projections are diagonal between the amplitude and complex conjugate. Acting on T (above the dashed line), these are [4] C [1] ab,cd = Correspondingly, projectionsC acting on the hard subdiagrams arẽ We now turn to the factorization of the spin degrees of freedom.

Spin projections and dimensional regularization
The spin analog of the color decomposition given in eq. (2.6) in four dimensions is a Fierz decomposition for each pair, with a normalization suitable to the basis in eq. (2.5), We shall see below that in calculations of T and H using dimensional regularization, we need to expand the basis of Dirac structures to reflect the arbitrary number of dimensions involved. The regularized cross section will include terms that are nonvanishing resulting from the product of poles associated with collinear radiation, a long-time process, times short-distance factors that vanish when the number of dimensions is taken to four. The introduction of new, evanescent projections [21,22] will enable us to identify such terms, and organize them appropriately. In particular, this analysis will be necessary to reconcile NLP factorization at next-to-leading order with existing NRQCD calculations in selected channels, and to provide unambiguous definitions for short-distance coefficients for channels in which explicit calculations have not been carried out. The full extended Dirac algebra in D dimensions is given by linear combinations of the elements of the sets where all of the (distinct) hatted indicesm 1 =m 2 · · · =m j are outside the usual four dimensions. We assume that all our continued dimensions are spacelike and that the trace of any product of Dirac matrices with all different indices vanishes, including a single matrix. Even though we usually think of dimensional continuation as infinitesimal, as long as = 2 − D/2 remains variable, there is no limit to the size of the algebra generated this way. Nevertheless, the number of elements of the algebra that can be realized in a given calculation is finite, since there are no more Dirac matrices than the number of vertices plus the number of propagators along the fermion lines in any given diagram. We will use the full set of elements of the sets above, {Γ I } including {ΓÎ j } to realize an orthonormal basis with which to expand both T and H. As a result, the traces that link the long-and short-distance parts will be nonzero only when the Dirac structure of the two factors matches exactly for the two functions. As we expand the short-distance part to higher orders, the expansion has more and more terms from {ΓÎ j }. For the project at hand, however, we only encounter three Dirac matrices with the evanescent projections. This makes the calculations described below more manageable than they might otherwise have been.
To facilitate the discussion of the new evanescent projections for our NLO calculations, we introduce the notation for any value of µ and ν, Of special interest for us in this calculation are the matrices with one and two hatted indices, wherem =n again. From this point on, we will call states associated with these evanescent projections, evaA and evaB, respectively. As observed in ref. [4], in the frame where the heavy quarkonium moves in the zdirection, as in eq. (1.3), the leading contributions of the long-distance functions T are proportional to matrices in with large + components dominate, i.e. terms projected out with Γ I ∼ n /. Therefore, in the limit m Q /E H → 0, we require in eq. (2.9) Γ I ∼ n / above the cut and Γ I ∼ p / below the cut. This is the case in arbitrary numbers of dimensions. The observed momentum p µ is kept in four dimensions, so that the full set of NLP matrices below the cut is given by suitably normalized elements of the set p / , p / γ 5 , p /γ mn , p /γmn , (2.13) in both for the amplitude and complex conjugate side. In the absence of polarizations in the initial state, the expansion of the hard scattering in terms of Dirac matrices is diagonal in these matrices, as anticipated in eq. (2.4) above. This diagonality is then carried over to the general spin decomposition given by eq. (2.9), which ensures that the fragmentation functions are also diagonal in the projections. New evanescent fragmentation functions based on these Dirac matrices must also be added to achieve exact NLP factorization in dimensional regularization. The first two matrices in eq. (2.13), which appear with both singlet and octet color, correspond to longitudinal and scalar polarization configurations for the pair. They are referred to as vector (v) and axial (a) projections. We can think of the additional matrices with hatted indices as extensions of p /γ 5 into D dimensions. For us, γ 5 remains in four dimensions. Of course, for D = 4, only the first two of these structures are necessary, and they appear in the Fierz identity, eq. (2.9). Transverse projections are also possible, of course, but are associated with even powers of Dirac matrices, and vanish at NLP for the gluon-associated processes we consider [4].
The resulting, normalized spin projection operators for the v, a, evaA, and evaB fragmentation functions are given by where repeated indices are summed. The corresponding projection operators for hard, short-distance coefficients are given byP In principle, there are separate evanescent functions for each combination of hatted indices, but for four-dimensional final spin and momentum states, they are all related by rotations. Using rotational invariance in dimensions beyond four, the functions are the same for every value of the hatted indices. As a result we sum over all hatted indices, dividing by the number of extra dimensions. Of course, this counting is only possible in integer numbers of dimensions, but it has an analytic continuation to any complex value of = 2 − D/2.
The overall factors for fragmentation projections in eq. (2.14) evaA and evaB come with overall counting normalizations (for the collinear fragmentation functions) , respectively. This counting normalization ensures the orthonormality condition with s, s = v, a, evaA, evaB. As long as such orthonormality condition is satisfied, the exact normalizations of P andP are not important. However, since both 2 D−4 and are ∼ 1 , depending on whether we choose to put this overall normalization with the hard projections or fragmentation projections we change the order of in the hard coefficients or fragmentation functions. We choose to include these counting factors in the projectors for the fragmentation functions, and use MS scheme for the evanescent as well as four-dimensional fragmentation functions. This compensates for the suppression induced by summing over the number of dimensions beyond four in eq. (2.14). For fixed hatted indices and independent normalizations, the corresponding fragmentation functions are proportional to 1/ , just as for four dimensional projections. Of course, our final results are scheme independent, but our choice is physically consistent with the above remark that evanescent contributions come from order short-distance parts convolved with poles in long-distance fragmentation functions.
With these normalizations, we can generalize the Fierz identity in four dimensions to where the primed indices are contracted with the hard scattering, and the unprimed with the fragmentation functions, and where the remainder, R, contributes only beyond NLP. R includes the transverse projections, which may also be NLP for other processes. Again, we can identify a finite list of projections because we are working at finite order in the hard scattering at NLP. Now we can separate both color and spin traces in the integrand of eq. (2.2), using eq. (2.6) and eq. (2.18), respectively, which, with standard approximations on loop momenta, give us the NLP part of the factorization found in eq. (1.2) [4].

Operator definitions for partonic fragmentation functions
Each pair of projections, eqs. (2.7) for color and (2.14) for spin, defines a fragmentation function. The result is given above in eq. (1.4). Compared to ref. [4], we now include in the set of projections the full list of D-dimensional Dirac structures in eq. (2.14).
As in the single-parton case, we will compute the short-distance coefficients for the heavy quark pair states from IR regulated, D-dimensional partonic cross sections, using perturbative and regulated partonic fragmentation functions. In these calculations, the heavy quark mass is set to zero. In contrast to the single-parton case, however, although the total momentum of the pair must be the same in the amplitude and complex conjugate of figure 1b, the relative momenta of the quark and antiquark within the two pairs need not be the same. This simply means that there can be different histories on how a heavy quark pair reaches the final heavy quarkonium state. The partonic implementation of the pair-to-hadron fragmentation function, eq. (1.4), is therefore slightly different than for single-parton fragmentation.
To compute partonic short-distance coefficients, we define fragmentation functions that can be directly applied to products of amplitudes in the infrared-regulated theory. These partonic fragmentation functions depend on fractional momenta of the heavy quarks in the final state. They are given by, again suppressing the gauge links, where p + is the lightlike projection of the heavy quark pair momentum, and in both final states, p Q + pQ = p. In these matrix elements, the "heavy" quarks are treated as massless, since they are designed to match the long-time behavior of the products of amplitudes in the infrared-regulated theory. The produced pair(s) QQ(κ , u ) are labelled explicitly by the relevant heavy quark momentum fraction. In each of these Fock states we take the final heavy quark and antiquark momenta to be lightlike and parallel to the projection of the heavy quark pair's momentum, p (2.20) Normalizing the quark-pair states to the projection operators in eq. (2.14) we have at LO, This simple result makes the relation between the LO short-distance coefficients and the Born cross section eq. (1.2) direct, as we will see below.
The partonic fragmentation functions in eq. (2.19) were employed in ref. [4] to identify evolution kernels. Here, we use them to compute partonic fragmentation functions at order α s , and it is useful to have an explicit definition. In principle, it is possible to combine these off-diagonal fragmentation functions into a more direct analog of the hadronic functions in eq. (1.4). This would be carried out by combining linear combinations of states with fixed u and v into some "pure" partonic final state characterized only by the total momentum of the pair. This is not necessary for our calculations here, however, so we will not develop the formalism in this direction.

Partonic fragmentation functions at order
Real emission diagrams for computing partonic fragmentation functions at order α s for light-cone gauge. We label the flow of + components of momentum to make z, u, v, u , v variable dependence explicit. For the full momentum dependence and details of the computation, see ref. [4].
At order α s , the partonic fragmentation functions for the four dimensional basis were calculated in ref. [4] for the four-dimensional projections s = v, a of eq. (2.14). Here, we give the set of partonic fragmentation functions needed for the factorized NLP cross section at O(α s ). The only difference is in the Dirac projections, and we do not give the details of the computation here. Recalling that our final states are always four-dimensional, we only need fragmentation functions with κ = v, a in color octet or singlet for NLO calculations of our process. The MS scheme results can be written as where the kernels P κ→κ for the non-evanescent states are given in the appendix of [4].
We extend the results of [4] with the relevant evanescent fragmentation functions of our process, namely κ = evaA [8], evaB [8] and κ = v, a in color singlet or octet state. Since these are all off-diagonal splitting functions, at O(α s ), they are given by the real emission diagrams in figure 2 for light-cone gauge. We find where the functions S ± and ∆ [1,8] ± are defined by As we will demonstrate below, these evanescent fragmentation functions are crucial in our factorization to maintain scheme independence and agree with NRQCD calculation.

Cross sections and coefficient functions
With the order α s pair fragmentation functions in hand, we are ready to derive the shortdistance coefficients at LO and NLO. We set the stage by exhibiting the short-distance coefficients in terms of partonic cross sections. We then compute the partonic cross sections for the states of interest up to NLO. Finally, we use the partonic fragmentation functions found in section 2 to derive expressions for the NLO short-distance coefficients.

Coefficient functions from cross sections
To compute the NLP short-distance coefficients dσ (n) e + e − →[QQ(κ)] , we replace the heavy quarkonium H of eq. (1.1) by [QQ(κ)] and compute the left-hand and right-hand sides to n-th order in α s . We outline the procedure, in schematic notation, for the process under consideration.
For an arbitrary differential cross section, we write in the condensed notation of eq. (1.1), where p is the final heavy quark pair's momentum. In this expression, the arguments of the functions and the corresponding convolution symbols are taken as in eq. (1.2). We will suppress the explicit arguments of fractional momenta below. At lowest order for the gluon-associated production process shown in figure 3, at O(α s ) and with the final state [QQ(κ)]g, only the second term contributes, (3.2) Then, using eq. (2.21) for the zeroth order fragmentation functions, we arrive at 3) The short-distance coefficient is therefore equal to the partonic cross section at lowest order. Of course, for LO leptonic annihilation, κ is always an octet configuration. Evaluating eq. (3.1) at O(α 2 s ) in the same way, but with two gluons in the final state, gives (3.4) Note again that we only compute the real diagrams at O(α 2 s ) as we focus on the region x H away from 1. Then using eq. (2.21) for D (0) , we solve directly for the short-distance coefficient at order α 2 s , dσ (2) e + e − →[QQ(κ)]gg = dσ where now κ can be singlet or octet. 2 The subtraction term on the right-hand side removes the long-distance behavior in the O(α 2 s ) cross section, which is due to the collinear emission of a gluon by the pair, a process that requires times that are large compared to 1/E H .

Partonic projections and phase space
We next establish the notation and projections necessary to calculate the partonic cross sections on the right-hand side of eqs. (3.3) and (3.5). As in the discussion of partonic fragmentation functions, eq. (2.19), we will use the projection operators of eq. (2.15) to define our partonic cross sections, where the superscript m = 1 or 2 denotes the order of α s , which is also equal to the number of gluons radiated, corresponding to the m + 1-particle phase space factor dΠ m+1 .
Here, the amplitudesM are computed perturbatively, by stripping the final-state spinors for the pair and replacing them by a combination of the spin and color projections shown above. In explicit calculations to produce spin state s, and color state I, we thus make the replacement for s = v, a, evaA, evaB, and similarly for singlet and octet color states. The third line of eq. (3.6) defines the "squared" amplitude in this context. As for the perturbative computation of partonic fragmentation functions, only the total momentum of the pair is the same in the amplitude and complex conjugate. Finally, Lorentz-invariant phase space in eq. (3.6) is given in D = 4 − 2 dimensions, not including the symmetry factors, by where k i = (E k i , k i ) and p i = (E i , p i ) label the momenta of the incoming particles and the outgoing particles, respectively. For the LO diagrams of figure 3 and the NLO real diagrams in figure 4, eq. (3.6) then becomes and respectively. The extra 1/2 in eq. (3.10) is the symmetry factor associated with the two gluons in the final state. We parameterize phase space in terms of fractional energies, labelling the pair's total energy as E 1 . In particular, for three-particle phase space we use For an unpolarized initial state, the D-dimensional phase space factors can be integrated over angles to give where we define (3.14) The limits y = 0 and y = 1 correspond to one of the two final state gluons in figure  4 carrying away maximal energy, i.e. x 2 = 1 and x 3 = 1, respectively. This gives the configuration in which the other gluon becomes collinear to the observed heavy quark pair, giving rise to collinear singularities. We now turn to the calculation of matrix elements and cross sections.

Dirac traces and cross sections
We continue by separating the squared S-matrix into the leptonic and hadronic tensors as where the superscripts, (m), indicate the perturbative order in the strong coupling. Here, the leptonic part, taken at lowest order in QED coupling with spin averaging, given in terms of electron and positron momenta k i is As noted above, we integrate over angles at fixed fractional momentum, x 1 3 . For our x 1 distribution, using current conservation, q µ H µν = 0, with q µ = (k 1 + k 2 ) µ , we can simplify the hadronic tensor as (3.17) We find for quarks with fractional change e Q , where, in the second expression, we have factored out the scale, µ, and coupling dependence from −g µν H µν,(m) to define F (m) κ . We now combine the leptonic and hadronic tensors from eq. (3.16) and (3.18), respectively, in the squared matrix element in eq. (3.15). These results are combined with twoand three-particle phase space, eqs. (3.12) and (3.13), in eqs. (3.9) and (3.10) to derive expressions for the two-and three-particle cross sections, Here, we have introduced a convenient normalization, which reduces to the inclusive Born cross section for e + e − → µ + µ − in D = 4 dimensions. As follows from our use of the projection, eq. (2.18) in the computation of partonic fragmentation functions above, we take γ 5 to be strictly 4-dimensional to carry out the Dirac trace in computation of F (m) κ . This is consistent with Breitenlohner-Maison-'t Hooft-Veltman (BMHV) γ 5 scheme [39,40]. For LO, the heavy quark pair state κ can only be octet and we find We recall that, by eq. (3.3), the LO cross sections are identical to the LO short-distance coefficients that we use in the computation of the NLO short-distance coefficients (see eq. (3.5)). Going on to NLO, we compute F κ for κ = v, a in color singlet or octet state. Although the F (2) κ may be computed for the evanescent states, they will be subleading. To be specific, to compute non-evanescent dσ To clarify the structure of the y integral, we rewrite F (2) κ for κ = v, a in color singlet or octet state as where we have expanded to first order in as the higher order terms in do not contribute to the NLO cross section in D = 4. For the term χ (1) , we note that there are contributions coming from the additional D − 4 components of the momentum p 2 , denoted asp 2 . When there is γ 5 in the D-dimensional Dirac trace, such D − 4 components naturally appear in the BMHV scheme. As demonstrated in Appendix A, they can simply be replaced bŷ To be explicit, F κ (x 1 , y, u, v, ) hasp 2 2 dependent terms with some coefficient b that can be replaced by eq. (3.31) as The term −b(x 1 , u, v)(1 − x 1 ) contributes to χ (1) κ in eq. (3.30) above. Using our decomposition of F (2) κ and also using the symmetry y ↔ 1 − y, we can explicitly separate eq. (3.20) into finite and pole pieces as To define the remaining functions on the right in eq. (3.33), we introduce the notation, in terms of which they are given by

Explicit short-distance coefficients
We can now present explicit results of short-distance coefficients differential in x 1 . From eq. (3.3), LO short-distance coefficients differential in x 1 are identical to the corresponding partonic cross section, i.e.
At NLO, we use eq. (3.5) to derive the x 1 differential results, where we now make the MS scheme dependence explicit. To compute eq. (3.45) we make use of the results in eqs. (3.26) -(3.29) for the Born cross sections at fixed u and v, and MS partonic fragmentation functions given in eqs.
(2.23) -(2.30) and [4]. We give explicit subtraction terms for κ = v, a in color singlet or octet state 46) where N (x 1 , u, v) and X κ (x 1 , u, v) are defined above in eq. (3.34). The functions Z κ (x 1 , u, v) define the extra finite pieces of the subtracton terms, and are given by Then, combining the cross section in eq. (3.33) with its subtraction given by eq. (3.46), according to eq. (3.45), our explicit NLO short-distance coefficients for vector or axial states in MS scheme take the form (3.51) We note that omission of evanescent fragmentation functions from the sum over ζ in eq. (3.46) would change the values of Z κ (x 1 , u, v), but would not change X κ (x 1 , u, v) since dσ (1) ζ dx for evanescent state ζ is linear in . Therefore, at this order, evanescent subtractions are not necessary in subtracting the poles. However, they subtract the finite terms that are sensitive to long-distance dynamics, resulting from a long-distance pole times a term proportional to in the hard-scattering function. In fact, at higher orders they would be needed even to subtract the IR poles consistently. At higher loops, there are terms with the same LO hard part proportional to , multiplied by multiple poles of higher order evanescent fragmentation functions.

Comparison to NRQCD
The determination of short-distance coefficients in the previous section is consistent with any model of factorized long-distance behavior. As noted in the introduction, fragmentation functions for the heavy quark pair have been computed in refs. [17,18], assuming the applicability of NRQCD factorization to the corresponding matrix elements. In this section, we shall make use of this assumption, which is certainly valid at NLO. By combining the fragmentation functions of refs. [17,18] with the pQCD short-distance coefficients derived in the previous section, eqs. (3.44) and (3.51), we find cross sections for the relevant NRQCD channels that can be compared directly to cross sections in NRQCD at the O(α 2 α 2 s ). For two of the channels, 3 S [1] 1 [29,31,32] and 1 S [8] 0 [33], we are lucky enough to have explicit NRQCD calculations at the O(α 2 α 2 s ) with which to compare. To this order, these channels begin at NLP, and we will see that the comparison is relatively direct. We expect and will confirm that the cross sections derived as above, by combining pQCD shortdistance coefficients and NRQCD fragmentation functions in the 3 S 1 and 1 S [8] 0 channels, have high-energy behaviors that agree precisely with those of the NRQCD calculations. This equality, while in principle straightforward, requires consistent renormalization procedures along with systematic treatment of the evanescent partonic fragmentation functions introduced in section 2. We begin our discussion by reviewing the relationship between NLP factorization and fixed-order NRQCD calculations.

Relating NLP factorization to fixed-order NRQCD
Non-relativistic QCD can be applied with or without the presence of a perturbative scale beyond the heavy quark mass. NRQCD treats the heavy quark mass, m Q , as a hard scale, separating amplitudes and cross sections for these processes into relatively short-distance coefficients, associated with hard scales at and above O(m Q ), and universal long-distance matrix elements (LDMEs) associated with soft scales at and below O(m Q v). Schematically, the cross section can be factorized at NRQCD factorization scale µ Λ ∼ m Q as wheref QQ(ν) is a perturbative coefficient describing the production of a heavy quark pair in NRQCD state ν, and O is the LDME describing the formation of the observed heavy quarkonium H from the heavy quark pair state ν. These LDMEs are scaled in powers of the heavy quark pair's relative velocity v 1. Applications of NRQCD have resolved tensions between theoretical predictions and experimental measurements, using a limited number of LDMEs [2,3]. Puzzles remain, however, especially involving polarization and the tension between LDME values required to fit different production processes [2,41,42]. Since these tensions often involve some larger perturbative scale, E H , where heavy quark pair fragmentation is important [19], it is natural to reconsider these cross sections in the extended pQCD formalism discussed above.
We can carry out matching for our pQCD fragmentation functions in eq. (1.4) in terms of the same universal LDMEs of NRQCD at the input scale µ 0 ∼ 2m Q , [4,5,17,18] as where κ and ν label a relativistic pQCD state and a non-relativistic NRQCD state, respectively, also appearing in eq. (1.2) and in eq. (4.1). Here, the matching coefficientŝ d QQ(κ)→QQ(ν) (z, u, v; m Q , µ 0 , µ Λ ) are perturbatively computable. This matching is analogous to the NRQCD treatment of single-parton fragmentation, as developed in refs. [43,44], for i = q,q, g. Using these models for fragmentation functions at µ = µ 0 , we can express the single-particle inclusive cross section in the notation of eq. (1.1) at any fixed order 4 in terms of a limited number of non-perturbative NRQCD LDMEs as Once again, we emphasize that the true motivation for NLP pQCD factorization comes from resummation of large logarithms of ln where we suppress dependence on variables other than the factorization scale µ = µ 0 . Expanding in relative velocities up to v 4 , the full set of LDMEs associated with heavy quark pair states is (4.5) This set provides predictions in terms of only a few NRQCD parameters. As emphasized in refs. [17,18], computing fragmentation functions using NRQCD in principle enables us to go from fitting several non-perturbative LP and NLP fragmentation functions to determining a small number of LDMEs. Our discussion here checks the consistency of this procedure.
Comparing the basic NRQCD relation, eq. (4.1) with the high-E H factorization with NRQCD input for the fragmentation functions, eq. (4.4), one can expect an order-by-order relation between NRQCD and NLP factorization short-distance coefficient functions, (4.6) Clearly, the full NRQCD calculation contains more information than the NLP factorized cross section at the fixed order in α s . All such information, however, appears beyond NLP. This was confirmed numerically for high-p T production at hadronic colliders in selected channels by ref. [19], using LO pQCD short-distance coefficients. Such comparisons serve both as a test of the NLP formalism, and as a tool for studying the approach to high-energy behavior. In the remainder of this section, we use the results of section 3 to analytically confirm eq. (4.6) for the channels 3 S

Input fragmentation functions and factorized cross sections
For a given NRQCD channel, we compute the fixed-order pQCD prediction of the x H = 2E H /Q differential version of eq. (4.4) to O(α 2 α 2 s ). To do so, we need matching coefficients for fragmentation functions,d, in eq. (4.3) to O(α s ). As our gluon-associated process will only involve NLP scaling to the order we consider, we only need to consider NRQCD factorization of the pair fragmentation functions.
It will be convenient below to factor out mass and coupling dependence compared to thed in eq. (4.2), defining perturbative coefficientsd (i) by Note that the input fragmentation function does not have u and v dependence that is present in the partonic fragmentation functions in eq. (2.22). This is because the final state in the input fragmentation function is a heavy quarkonium. Of course, unlike the partonic fragmentation function, the input fragmentation function is also a non-perturbative object with perturbative matching coefficient extracted. We now apply the notation of eq. (4.7) to the general cross section eq. (4.4) for heavy quarkonium production at fixed x H . We can then isolate the contribution from a fixed intermediate NRQCD state ν to a heavy quarkonium production cross section, at fixed x H at order O(α 2 α 2 s ), as (4.8) The second relation defines the functionsk , which result from the convolutions ofσ (1) withd (1) andσ (2) withd (0) , respectively. Their normalization is set by separating the overall factor, σ 0 , with σ 0 given in eq. (3.21) above.
To facilitate the comparison of these results with direct NRQCD calculations, we rescale NRQCD coefficients,f in eq. (4.1), with the same overall normalization as in eq. (4.8), We conclude that for each channel ν, This is the O(α 2 s ) version of eq. (4.6). Here, the measure of higher-power corrections to the large-E H behavior is given by We will test eq. (4.10) below for channels in which an explicit NRQCD calculation is available.
The first step in verifying eq. (4.10) is to recall explicit results ford (0) andd (1) . These can be found in the appendices of [17,18] 5 . There, however, thed (1) were presented using the Kreimer γ 5 scheme of refs. [45][46][47] in MS subtraction scheme. To be consistent with our factorization procedure for the pQCD short-distance coefficients above, we must recompute these coefficients in the BMHV γ 5 scheme in MS subtraction scheme. To be self-contained, we present all the functionsd used in this paper, some of which are different from [17,18] due to the difference in γ 5 scheme. We present the ones needed for 3 S [1] 1 and 1 S [8] 0 here, but list all the other relevant ones in Appendix B. Suppressing the arguments, they are given asd Note that ∆ [1,8] ± are identical to those found in eqs. (2.32) and (2.33) with u = v = 1/2. Note also that we have dropped δ(1 − z) dependent terms in these NLO matching coefficients. This is because the LO cross sections dσ (1) /dx 1 are all proportional to δ(1−x 1 ). Then in eq. (4.8), any δ(1 − z) terms ind (1) contribute only at x H = 1, and thus are not included in this study.
We also do not need to compute input fragmentation functions for evanescent intermediate states. We recall that the LO short-distance coefficient functions given in eqs. (3.28) and (3.29) for evanescent states are proportional to . These terms contribute to the partonic cross section at order 0 because they multiply the infrared pole of the evanescent partonic fragmentation functions. All input fragmentation functions calculated from NRQCD, however, evanescent or four-dimensional, are finite after renormalization, and the poles of the partonic calculation are replaced by finite logarithms. The corresponding terms thus remain of order , and vanish in four dimensions.
Convolving the short-distance coefficients and thed presented here, we can write explicit expressions for thek (i) ν (x H ; µ 0 , µ Λ ) in eq. (4.8) for the relevant NRQCD channels. We find for 3 S [1] 1 and 1 S [8] 0 , and we take µ 0 = E H to remove large logarithms from the short-distance coefficients. Notice that there are additional logarithms coming from the threshold limit x H → 1. Meaningful comparison to the data near such endpoints can only be made after resumming these large logarithms, which can be done in principle by combining threshold resummation techniques to our work [48][49][50][51][52]. Aside from such threshold logarithms, now that large logarithms only remain in the fragmentation functions, one can try to evolve them to E H and resum the large logarithms of ln(E H /2m Q ). In this way, QCD factorization demonstrates how the natural scale choice appears.
Of course, both NRQCD and QCD calculations also have a renormalization scale, µ r . In NRQCD calculations of gluon-associated e + e − processes, the renormalization scale is chosen to be µ r = 2m Q or µ r = Q/2. For example, refs. [26][27][28] calculate the total cross section, integrated over energy E H , and are thus left with only m Q and Q as the relevant scales. For our energy fraction distribution, we will identify the renormalization scale with the factorization scale, µ r = µ 0 = E H , for both NLP predictions and NRQCD calculations. 1 and more recently from [33] for 1 S [8] 0 , we then find for x H < 1,

Comparisons for 3 S
(4.34) In summary, for the cases where it can be checked, the pQCD factorization formalism successfully reproduces the correct, and reasonably non-trivial, high-energy behavior of the full calculation. The results of Appendix C can also be used in eq. (4.8) to give new closed expressions for the high energy behavior of the channels for which explicit calculations do not exist.

Numerical results
In the following, we carry out a few numerical investigations of the results of the preceding sections. We begin by studying the x H distribution in the 3 S 1 channel, to illustrate the approach of the fixed order pQCD cross section to the full NRQCD result as E H increases relative to m H . In section 5.2, we study energy fraction x H distributions of different NRQCD channels for H = J/ψ, and compare to the x H distribution of Belle data for J/ψ production at Q = 10.6 GeV. In section 5.3, we explore the significance of the logarithmic corrections associated with the evolution of the heavy quark pair fragmentation functions. Ratio    ) are dropped, and thus the ratio produces a divergence. This is obviously a region where such power corrections are dominant. Although the sum (k (1) +k (2) )/f does not change under different factorization scale choice µ, we take µ = E H scale choice to remove large logarithms fromk (2) . We observe from figure 5 that the sum (k (1) +k (2) )/f can reproduce the NRQCD result quite accurately over a large range of x H already by Q = 30 GeV. According to eq. (4.33), the ratio of the sum (k (1) +k (2) )/f approaches unity as the energy increases at any fixed x H . These plots show how accurate the pQCD results are at these energies. To quantify further how closely the pQCD result approaches the NRQCD result, we define E down and E up as the minimum and the maximum value at which the ratio lies within 0.95 < (k (1) +k (2) )/f < 1.05, respectively.

Approach to full NRQCD result
In figure 6, we plot E down as a function of Q 2 . We find it more natural to plot it against Q 2 , rather than Q, since the maximum energy a heavy quarkonium can carry, As shown in figure 6, the pQCD result approaches NRQCD In figure 7, we plot (E up − E down )/(E max − E min ) to show the percentage of the energy range that describes the NRQCD result within 5 percent as a function of Q 2 . When Q 2 ≈ 25 GeV, more than 90 percent of the available range describes NRQCD result within 5 percent.
It is interesting to note that in the case of pp collisions at Tevatron energies, ref. [5] carried out a similar factorized NLP analysis based on an order α s short-distance coefficient, followed by order α s fragmentation. This is the analog of ourk (1) term above alone. In ref. [5], this term was sufficient to reproduce the O(α 2 s ) numerical NRQCD result reasonably well over a large range of p T . The results here show a similar qualitative agreement from k (1) , when the factorization scale is chosen as µ = E H .

Distribution in x J/ψ
Next, we compare the x H distributions of the various NRQCD channels for H = J/ψ, as computed from pQCD factorization through eq. (4.8). As noted in section 4, there are full calculations only for 3 S [1] 1 and 1 S [8] 0 channels in the literature, so that the remaining pQCD-based curves we will exhibit are in this sense new. Our intention is to compare the shapes, rather than the magnitudes of these contributions to the inclusive cross sections.
Nevertheless, evaluations of the cross sections due to the NRQCD channels require specific values for the long-distance matrix elements. We make the following choices. For the singlet, we adopt the value from ref. [53] For the octet channels, we take the values suggested as maxima for the treatment of this process in ref. [27], except that we also keep a nonzero 3 S [8] 1 matrix element for the purposes of comparison,   Figure 8 gives the x J/ψ -distributions found from eq. (4.8) for CM energies of 10.6, 30 and 100 GeV for these choices of matrix elements. For 10.6 GeV, we can compare to Belle data [23]. 7 From figure 8, we find that the color-singlet distribution has the same general shape as the data, decreasing gently toward zero as x J/ψ approaches unity. In contrast, all the non-negligible color-octet curves provide end-point enhancements. The curves retain these features at the higher energies. Meaningful comparisons to the data near such endpoints, however, can only be made after organizing large logarithms there [48][49][50][51][52].

Leading logarithms at two-loop order for
The results of section 5.1 suggest that pQCD factorization cross sections approach those of NRQCD rather quickly. The relative ease of computation for pQCD factorization hard parts, where m Q is taken to be zero, can facilitate the systematic computation of fixed order 6 Note that there is an additional factor 1/2 difference in our choice of normalization in the definition of singlet LDME. 7 We have transformed their momentum distribution to xH distribution  x J/ψ (GeV) NRQCD results at high energy. The pQCD factorization approach, however, provides not only the asymptotic behavior for NRQCD calculations, but also evolution equations [4] that organize logarithms of ln(E 2 H /(2m Q ) 2 ) to all orders. To be specific, the heavy quark pair fragmentation functions satisfy the evolution equations 1 to study the importance of resumming the logarithms, ln(E 2 H /(2m Q ) 2 ). We restrict ourselves to the color-singlet channel case as a test. Note that at the order we consider for gluon-associated processes, the heavy quark pair can only be created at short distances, and we need only the evolution of the heavy quark pairs among themselves by eq. (5.4). The general problem includes mixing between the heavy quark pairs and the single partons [4]. Equation (5.4) has the following two-loop order solution for the fixed heavy pair NRQCD state 3 S which can be easily checked perturbatively [56,57]. As appropriate to our factorization, we set our input scale µ 0 = 2m Q and the hard scale µ = E H . These choices remove large logarithms from the short-distance coefficients and keep them in the fragmentation functions. Their ratio specifies the logarithms we would like to resum, ln(E 2 H /(2m Q ) 2 ). The fragmentation functions at the input scale µ 0 are found as in eq. (4.7), where we now specify the NRQCD state ν there to be 3 S As we only concern ourselves with the leading logarithms, unsuppressed by further powers of α s , we use the tree level matchingd (0) . Since the only nonzerod is κ = v [1], the final pQCD state in the kernel must always be v [1].
To get the x H distribution, we must do a further convolution ⊗ z;u,v of the two-loop solution in eq. (5.6) with short-distance coefficients computed in section 3. Again, because we do not want further suppression in α s , we only need to convolve with the leading-order short distance coefficients found in eq. (3.26) and (3.27), and thus κ in eq.(5.6) must be either v [8] or a [8]. From now on, we suppress the arguments of the short-distance coefficients, splitting functions, and fragmentation functions for simplicity. Calculating the single logarithms from the perturbative solution gives an expected result, in agreement with eq. (4. 23 GeV, E H and 2m Q do not create a strong hierarchy and the resummation of such logarithms is not so important. However, as can be seen from figure 9, the leading logarithms at twoloop order modify the results up to ∼ 30% already at Q = 30 GeV. Such a large correction implies that solving the evolution equation to resum the logarithms ln(E 2 H /(2m Q ) 2 ) is crucial for a reliable prediction at higher energies than at Belle.

Conclusions
In this paper, we presented the first NLO calculation of short-distance coefficients in the context of NLP perturbative QCD factorization for heavy quarkonia, extending the quark pair fragmentation formalism developed in ref. [4]. We calculated short-distance coefficients as closed expressions for the physically-relevant x H distribution of heavy quarkonia in e + e − annihilation. We showed that it is useful to include evanescent operators, which are absent in the four-dimensional theory, when using dimensional regularization. At NLO, the contribution from the evanescent operators organizes finite pieces that are associated with long-distance dynamics. This organization is important to clarify comparisons with fixed-order NRQCD calculations.
Combining our calculation of NLO short-distance coefficients with NRQCD factorization of the heavy quark pair fragmentation functions [17,18], we derived the high energy behavior of the x H distributions in fixed-order NRQCD for various channels. We presented numerical results of these channels at various CM energies, which exhibit end-point enhancements for the octet channels and a large contribution from the singlet channel away from the tail. Using the singlet channel, we numerically illustrated that fixed-order NRQCD results rapidly approach their high energy behavior and demonstrated the im-portance of resumming logarithms of ln(E 2 H /(2m Q ) 2 ) by studying the numerical impact of including the leading logarithms at two-loop order.
We also showed explicit analytical agreement between our derived high energy behavior and the fixed-order NRQCD calculations in the literature for the 3 S [33] channels, illustrating the need for including the evanescent operators identified. The formalism developed in our work provides the groundwork for explicit NLO and higher order calculations of short-distance coefficients for many other factorized cross sections at next-to-leading power. We anticipated applications to angular distributions in e + e − annihilation, and to many single-particle cross sections for heavy quarkonia in hadronhadron scattering.

B Matching coefficients of input heavy quark pair fragmentation functions in BMHV scheme
In this appendix, we summarize all the matching coefficients of input heavy quark pair fragmentation functions used in this paper with BMHV γ 5 and MS subtraction scheme. The calculation of the same coefficients with Kreimer γ 5 and MS subtraction scheme can be found in the appendices of [17,18]. The additional terms in BMHV γ 5 scheme are indicated by the square brackets below. 2 )],z =1 + ∆  k (1)