Gluon fragmentation into quarkonium at next-to-leading order

We present the first calculation at next-to-leading order (NLO) in $\alpha_s$ of a fragmentation function into quarkonium whose form at leading order is a nontrivial function of $z$, namely the fragmentation function for a gluon into a spin-singlet S-wave state at leading order in the relative velocity. To calculate the real NLO corrections, we introduce a new subtraction scheme that allows the phase-space integrals to be evaluated in 4 dimensions. We extract all ultraviolet and infrared divergences in the real NLO corrections analytically by calculating the phase-space integrals of the subtraction terms in $4-2\epsilon$ dimensions. We also extract the divergences in the virtual NLO corrections analytically, and detail the cancellation of all divergences after renormalization. The NLO corrections have a dramatic effect on the shape of the fragmentation function, and they significantly increase the fragmentation probability.


I. INTRODUCTION
The production of a hadron in a high energy collision is in general an extremely complicated problem dominated by nonperturbative aspects of QCD. There are several ways to simplify the problem in order to make a theoretical analysis more tractable. One way is to consider inclusive production of the hadron, summing over all possible additional hadrons in the final state. Another simplification is to consider the production of the hadron with transverse momentum p T that is much larger than the momentum scale Λ QCD of nonperturbative effects in QCD, so that there are aspects of the problem that involve the small coupling constant α s (p T ). Another simplification is to consider a hadron whose constituents include a heavy quark whose mass m is much larger than Λ QCD , so that there are aspects of the problem that involve the small coupling constant α s (m). If the hadron is a heavy quarkonium, whose constituents are a heavy quark and antiquark, there are further simplifications from the typical relative velocity v of the constituents being small compared to 1. The most theoretically tractable problem is the one in which all these simplifying features are combined: the inclusive production of quarkonium at large p T .
A rigorous factorization theorem for the inclusive production of a single hadron at large p T was derived by Collins and Soper in 1981 [1]. It states that in the inclusive cross section for producing a hadron H at p T ≫ Λ QCD , the leading power in the expansion in powers of Λ QCD /p T can be expressed as a sum of perturbative QCD (pQCD) cross sections for producing a parton convolved with fragmentation functions: The sum extends over the types of partons (gluons, quarks, and antiquarks). The pQCD cross sections dσ are essentially inclusive cross sections for producing the parton i, which can be expanded in powers of α s (p T ), convolved with parton distributions if the colliding particles are hadrons. The nonperturbative factors D i→H (z) are functions that give the probability distribution for the longitudinal momentum fraction z of the hadron H relative to the parton i. The symbol "⊗" in Eq. (1) represents an integral over z. Evolution equations for the fragmentation functions can be used to sum large logarithms of p T /Λ QCD to all orders in α s . The factorization formula in Eq. (1) applies equally well to heavy quarkonium with m ≫ Λ QCD . A proof of this factorization theorem that deals with issues specific to quarkonium production was first sketched by Nayak, Qiu, and Sterman in 2005 [2]. It gives the leading power (LP) in the expansion in powers of m/p T and it applies only at p T ≫ m. We will refer to this factorization theorem as the LP factorization formula. In the case of cross sections summed over quarkonium spins, the corrections are suppressed by a power of m 2 /p 2 T multiplied by logarithms of p T /m. The LP factorization formula has limited predictive power, because the nonperturbative factors D i→H (z) are functions of z that must be determined from experiment.
In 1994, Bodwin, Braaten, and Lepage proposed the NRQCD factorization formula, which uses an effective field theory called nonrelativistic QCD to separate momentum scales of order m and larger from momentum scales of order mv and smaller. The theoretical status of the NRQCD factorization conjecture is discussed in Ref. [3]. The NRQCD factorization formula states that the inclusive cross section for producing a quarkonium state H can be expressed as the sum of pQCD cross sections for producing a QQ pair with vanishing relative velocity multiplied by NRQCD matrix elements: The sum extends over the color and angular-momentum channels of the QQ pair. The pQCD cross sections dσ are essentially inclusive cross sections for producing the QQ pair, which can be expanded in powers of α s (m), convolved with parton distributions if the colliding particles are hadrons. The nonperturbative factors O H n are multiplicative constants that can be expressed as vacuum expectation values of four-fermion operators in nonrelativistic QCD [4]. They scale as definite powers of the typical relative velocity v of the Q orQ in H. The NRQCD matrix element O H n is essentially the probability for a QQ pair created in the state n to evolve into a final state that includes the quarkonium H. Through heroic calculations over the past decade, the inclusive pQCD cross sections for producing a QQ pair in all the most phenomenologically relevant channels have been calculated to next-to-leading order (NLO) in α s , even for the most difficult case of hadron collisions [5][6][7]. In many cases, the NLO corrections are very large, which suggests that higher order corrections may also be important. However the NLO calculations are sufficiently difficult that further improvement of the accuracy to next-to-next-to-leading order seems to be out of the question.
The predictive power of the LP factorization formula in Eq. (1) can be increased by applying the NRQCD factorization conjecture to the fragmentation functions. This reduces the nonperturbative factors from functions of z to multiplicative constants. The fragmentation function for the parton i to produce the quarkonium H is expressed as a sum of functions of z that can be calculated using pQCD multiplied by NRQCD matrix elements: The sum extends over the color and angular-momentum channels of a nonrelativistic QQ pair. The pQCD functions d i→(QQ)n (z) can be expanded in powers of α s (m). The nonperturbative factors are NRQCD matrix element O H n . If the LP/NRQCD factorization formula obtained by inserting Eq. (3) into Eq. (1) is expanded in powers of α s (m), it should reproduce the leading power in the expansion in powers of m/p T of the NRQCD factorization formula in Eq. (2). The usefulness of the LP/NRQCD factorization formula has proved to be limited at present collider energies. Explicit calculations using the NRQCD factorization formula have revealed that, in some channels, the LP cross section is not the largest contribution until p T is almost an order of magnitude larger than m [8].
An important recent development is the derivation of a factorization theorem that extends the LP factorization formula in Eq. (1) to the next-to-leading power (NLP) of m 2 /p 2 T . This factorization theorem was proven by Kang, Qiu, and Sterman [9,10]. A similar factorization formula has been derived by Fleming, Leibovich, Mehen, and Rothstein using soft collinear effective theory [11,12]. In the NLP factorization formula, the terms suppressed by m 2 /p 2 T are expressed as a sum of pQCD cross sections for producing a collinear QQ pair convolved with double-parton fragmentation functions, which are nonperturbative probability distributions in the longitudinal momentum fraction of the quarkonium H relative to the QQ pair. The predictive power of the NLP fragmentation formula can be dramatically increased by applying the NRQCD factorization formula to the double-parton fragmentation functions as well as the single-parton fragmentation functions [13,14]. In the case of cross sections summed over quarkonium spins, the corrections are suppressed by a power of m 4 /p 4 T multiplied by logarithms of p T /m. If the resulting NLP/NRQCD factorization formula is expanded in powers of α s (m), it should agree with the first few terms in the expansion in powers of m/p T of the NRQCD factorization formula.
The NLP/NRQCD factorization formula opens the door to dramatic improvements in the accuracy of theoretical predictions for quarkonium production at very large p T . The factorization formula can be expressed as a triple expansion in powers of α s , v, and m/p T . NLP factorization incorporates subleading powers of m/p T . The NRQCD expansion includes subleading powers of v. Accurate predictions also require including subleading powers of α s . The aspects of the problem that are perturbative are the pQCD cross sections for producing single partons, the pQCD cross sections for producing collinear QQ pairs, the coefficient functions of NRQCD matrix elements in the NRQCD expansions of both the single-parton fragmentation functions and the double-parton fragmentation functions, and the evolution kernels for both sets of fragmentation functions. It would be desirable to have all these ingredients calculated to NLO in α s .
The NRQCD-expanded LP fragmentation formula was actually first applied to quarkonium production at large p T back in 1993, when the first fragmentation functions for S-wave quarkonium states were calculated to leading order (LO) in α s for channels that are leading order in v. The fragmentation functions for a gluon into spin-singlet and spin-triplet Swave states at leading order in α s and v were calculated by Braaten and Yuan [15,16]. The fragmentation function for a heavy quark into spin-singlet and spin-triplet S-wave states at leading order in α s and v were calculated by Braaten, Cheung, and Yuan [17]. The fragmentation functions have since been calculated at LO in α s for all the color and angular-momentum channels that are predicted by NRQCD factorization to be most phenomenologically relevant. For the color-octet 3 S 1 channel, in which the LO fragmentation function is proportional to δ(1 − z), the fragmentation function has been calculated to nextto-leading order (NLO) in α s [13,18]. In this paper, we present the first NLO calculation of a fragmentation function whose form at LO is a nontrivial function of z: the fragmentation function for a gluon into a spin-singlet S-wave state η Q at leading order in v.
The outline of our paper is as follows. In Section II, we calculate the LO fragmentation function using Feynman rules introduced by Collins and Soper. We also define some quantities that are useful in the NLO calculations. In Section III, we explain how the calculation of the real NLO corrections can be simplified by introducing subtraction terms that cancel all the ultraviolet and infrared divergences, allowing the phase-space integrals to be calculated in 4 dimensions. We calculate the phase-space integrals of the subtraction terms analytically using dimensional regularization. In Section IV, we describe the calculation of the virtual NLO corrections, and we present analytic results for the poles that arise from the dimensionally regularized loop integrals. In Section V, we show how all the poles from phase space integrals and from loop integrals are cancelled by renormalization of the parameters of QCD and renormalization of the operator whose matrix element defines the fragmentation function. Some numerical illustrations of our results are presented in Section VI. We discuss the prospects for the NLO calculation of all the other phenomenologically relevant fragmentation functions in Section VII. Some details of the calculation of integrals at NLO are presented in appendices. In Appendix A, we derive parametrizations of massless twoparton phase space integrals that are used to integrate the subtractions terms for the real NLO corrections. In Appendix B, we present the pole terms in the loop integrals with an eikonal propagator that appear in the virtual NLO corrections.

II. LEADING-ORDER FRAGMENTATION FUNCTION
In this section, we calculate the perturbative fragmentation function for g → QQ, with the QQ pair in a color-singlet 1 S 0 state, at leading order in α s . We also introduce some related expressions that are useful in the calculation at next-to-leading order in α s .

A. Feynman rules
Gluon fragmentation functions can be calculated using Feynman rules derived by Collins and Soper in 1981 [1]. The fragmentation function is expressed as the sum of all possible cut diagrams of a particular form. The diagrams include an eikonal line that extends from a gluon-field-strength operator on the left side of the cut to a gluon-field-strength operator on the right side. Single virtual gluon lines are attached to the operators on the left side and the right side. The two virtual gluon lines from the operators are connected to each other by gluon and quark lines produced by QCD interactions, with possibly additional gluon lines attached to the eikonal line. The cut passes through the eikonal line, the line for the particle into which the gluon is fragmenting, and possibly additional gluon and quark lines. An example of a cut diagram with the cut passing through the lines of a heavy quark and antiquark and an additional gluon is shown in Figure 1.
The Feynman rules for the cut diagrams are relatively simple [1]. The 4-momentum K of the gluon that is fragmenting enters the diagram through the operator vertex on the left side of the eikonal line and it exits through the operator on the right side. Some of that momentum flows through the single virtual gluon attached to the operator and the remainder flows through the eikonal line. The operator at the left end of the eikonal line is labelled by a Lorentz index µ and a color index a. The operator at the right end of the eikonal line is labelled by a Lorentz index ν and a color index b. The Feynman rules can be summarized as follows: • If the single virtual gluon line attached to the operator at the left end of the eikonal line has momentum q, Lorentz index λ, and color index c, the Feynman rule for the operator is −i(K.ng µλ − q µ n λ )δ ac , where n is a light-like 4-vector.
• The attachment of an additional gluon line with Lorentz index β and color index c to an eikonal line with color indices d and e to the left and right of the attachment has the Feynman rule g s f cde n β .
• The propagator for an eikonal line carrying momentum q is i/(q.n + iǫ). The Feynman rule for a cut eikonal line carrying momentum q is 2πδ(q.n).
• The remaining Feynman rules are those of QCD.
The particle into which the gluon is fragmenting has a specified 4-momentum. In the case of fragmentation of a gluon into quarkonium, it is convenient to express that 4-momentum as 2p. The longitudinal momentum fraction z of the quarkonium is z = (2p).n/K.n.
The fragmentation function is the sum of all cut diagrams contracted with −g µν and δ ab and multiplied by the Collins-Soper prefactor [1] The factors in the denominator include the number of color and spin states of a gluon in D = 4 − 2ǫ dimensions. The factor of z D−3 arises from an integral over a transverse momentum.

B. NRQCD Factorization
The NRQCD factorization formalism [4] can be used to expand the fragmentation function for producing a quarkonium state into a sum of matrix elements of NRQCD operators multiplied by perturbatively calculable coefficients. The NRQCD matrix elements scale as definite powers of the relative velocity v of the heavy quark in the quarkonium. For a 1 S 0 quarkonium state η Q , the matrix element that is leading order in v was denoted [4]. Within the vacuum-saturation approximation, it can be interpreted as proportional to the square of the wavefunction at the origin for the quarkonium: The NRQCD factorization conjecture asserts that its coefficient, which is a function of z, can be calculated as a power series in α s (m), where m is the heavy quark mass. If we keep only the O 1 ( 1 S 0 ) η Q term, the fragmentation function for g → η Q can be expressed as The function d LO (z) in the leading-order term was calculated by Braaten and Yuan in 1993 [15]. Our goal is to calculate the function d NLO (z) in the next-to-leading order term. This function also depends on renormalization and factorization scales that have been suppressed in Eq. (7). The formation of the quarkonium η Q from the fragmentation of a gluon involves nonperturbative effects that are represented by the sum of infinitely many Feynman diagrams. Thus the fragmentation function D g→η Q (z) can not be calculated directly using perturbative QCD. However the coefficient of O 1 ( 1 S 0 ) η Q in D g→η Q (z) can be determined from the perturbative calculation of the fragmentation function for producing an appropriate QQ state at a fixed order in α s . The simplest choice is a QQ pair in a color-singlet spin-singlet state with zero relative momentum. Its angular momentum quantum numbers are therefore 1 S 0 . The perturbative fragmentation function D g→QQ (z) for producing the QQ pair has the same form as in Eq. (7) but with a different prefactor: It can be calculated from cut diagrams for the gluon fragmentation function in which the cut lines include Q andQ. Given the normalization of the NRQCD operator O 1 ( 1 S 0 ) defined in Ref. [4], the NRQCD matrix element for the QQ pair is If dimensional regularization is used to regularize ultraviolet and infrared divergences, this matrix element has no NLO corrections. By dividing the perturbatively calculated fragmentation function D g→QQ (z) by 2N c , we can obtain the coefficient of The perturbative fragmentation function D g→QQ (z) for producing a color-singlet spinsinglet QQ pair with zero relative momentum can be conveniently obtained by replacing the spinors from cuts through the heavy quark and antiquark lines in a cut diagram by projection matrices. The replacement rule for the product of the spinors to the left of the cut is This projection matrix is the product of a color matrix with explicit indices i and j that projects the QQ pair into a color-singlet state and a Dirac matrix that projects it into a 1 S 0 state. The projection matrix Γ ij on the right side of Eq. (10) satisfies Tr(Γ ij γ 0 Γ † ij γ 0 ) = 4m, which is the standard relativistic normalization for a particle of mass 2m in its rest frame. With the projection matrix in Eq. (10), the Dirac structure on each side of the cut reduces to the trace of Dirac matrices that include a single factor of γ 5 .
We use dimensional regularization in D = 4 − 2ǫ dimensions to regularize ultraviolet and infrared divergences. Since the conventional definition of γ 5 is specific to 4 dimensions, there is the possibility of an incompatibility between the definition of γ 5 and dimensional regularization. One property of γ 5 that we will use is that the trace of a product of γ 5 and fewer than four gamma matrices is 0. In the LO and NLO diagrams for the fragmentation function, this property can be used to reduce the Dirac trace in the amplitude on the left side of the cut to Tr([γ µ , γ λ , γ ρ ]p /γ 5 ), where [γ µ , γ λ , γ ρ ] is the antisymmetrized product of three gamma matrices whose 6 terms have coefficients +1 or −1. The Dirac trace on the right side of the cut can similarly be reduced to Tr([γ ν , γ σ , γ τ ]p /γ 5 ). After integrating over the momentum of the radiated gluon, the only independent tensors that can be contracted with the product of these Dirac traces to give a scalar are g µν g λσ g ρτ and g µν g λσ n ρ n τ . In 4 dimensions, the Dirac trace from the left side of the cut is The Dirac trace on the right side of the cut gives a similar expression with Levi-Civita tensor ǫ νστ β . In 4 dimensions, the product of ǫ µλρα and ǫ νστ β can be expressed as an antisymmetrized sum of products of four metric tensors with 24 terms. With some of the more common prescriptions for γ 5 , these metric tensors can be interpreted as those for D dimensions. In this case, the contractions of the two independent tensors with the product of the two traces reduces to The study of alternative prescriptions for γ 5 can ultimately be reduced to its effects on these two expressions.

C. Born fragmentation function
The fragmentation function for g → QQ can be calculated perturbatively from the cut diagrams in which the cut lines include Q andQ. At leading order in α s , the cut diagrams are the diagram in Figure 1 and three other diagrams. One of the other diagrams is obtained by interchanging the vertices where the gluon from the operator and the final-state gluon attach to the quark line on the left side of the cut. The other two are obtained by making a similar interchange on the right side of the cut. The cut lines are those for the Q andQ, the final-state gluon, and the eikonal line. The final-state Q andQ are on-shell with equal momenta p and total longitudinal momentum fraction z. The final-state gluon is on-shell with a momentum q whose phase space must be integrated over. The cut through the eikonal line gives a factor of 2πδ(K.n − (2p + q).n).
The amplitude corresponding to the sum of the two diagrams on the left side of the cut can be written down using the Feynman rules: where i, j, and c are the color indices of the final-state Q,Q, and gluon. After replacing the spinors by the projector in Eq. (10) and using the fact that the trace of the product of γ 5 and fewer than four gamma matrices is 0, the amplitude can be reduced to The cut diagram is obtained by multiplying this amplitude, whose free indices are µ and a, by its complex conjugate with indices ν and b, integrating over the phase space of the gluon, and summing over its color and spin states. The fragmentation function D g→QQ (z) is then obtained by contracting the cut diagram with δ ab (−g µν ) and multiplying by the Collins-Soper prefactor in Eq. (5).
We denote the product of the differential phase space for the final-state gluon with momentum q and the factor 2πδ(K.n − (2p + q).n) from the cut through the eikonal line by dφ Born . We use dimensional regularization with D = 4 − 2ǫ dimensions to regularize both ultraviolet and infrared divergences. The integral over q.n can be evaluated using the delta function from the cut through the eikonal line. After integrating over the angles of the transverse components of q, dφ Born reduces to a single differential: where s is the invariant mass of the QQg system: In Eq. (15), there is an implied Heavyside theta function that imposes the constraint s > 4m 2 /z. The fragmentation function for g → QQ at leading order in α s can be expressed as where N CS is the Collins-Soper prefactor in Eq. (5) and the function A Born in the integrand is We will refer to this function as the Born squared amplitude. The LO fragmentation function in D dimensions is Since the integral in Eq. (19) has no divergences, we can set ǫ = 0. The LO fragmentation function in 4 dimensions reduces to After evaluating the integral over s, the final result for the LO fragmentation function is The NRQCD matrix element, which is given by Eq. (9), can be inserted by multiplying by . Comparing with Eq. (8), we can read off the function d LO (z) in the fragmentation function for g → QQ: This same function d LO (z) appears in the fragmentation function for g → QQ in Eq. (7). The leading-order fragmentation function calculated by Braaten and Yuan in 1993 [15] can be reproduced by inserting the expression for

D. Born squared amplitudes with uncontracted Lorentz indices
To facilitate the calculation of the NLO corrections to the fragmentation function, it is convenient to generalize the integration measure for the LO fragmentation function in Eq. (19) by allowing q to be an arbitrary light-like vector. The Collins-Soper prefactor in Eq. (5) can be generalized to a function of p and q: The Born phase-space measure in Eq. (15) generalizes to where s = (2p + q) 2 . The product of N Born , dφ Born , and the function A Born (p, q) in Eq. (18) defines a LO differential fragmentation function with general light-like vector q: where s = (2p + q) 2 and z is the longitudinal momentum fraction z = (2p).n (2p + q).n .
If this measure is multiplied by a function of s and integrated over s from 4m 2 /z to ∞, it defines a function of z.
In the calculation of the real NLO corrections to the fragmentation function, it is convenient to have expressions for the Born squared amplitude with a pair of uncontracted Lorentz indices. They will be used to construct subtraction terms that cancel the ultraviolet and infrared divergences in the NLO corrections point-by-point in the phase space. Such amplitudes with uncontracted indices cannot be expressed as a linear combination of the contracted tensors in Eqs. (12b), but our prescription to extend γ 5 in D dimensions can still be used. As will become clear later, contributions from subleading terms in ǫ always originate from the Laurent expansion of subtraction terms involving the Born squared amplitude with no Lorentz indices. Hence the study of alternative prescriptions for γ 5 can indeed ultimately be reduced to its effects on the two expressions in Eqs (12b). There are two useful choices for the uncontracted indices µ and ν. One choice is the Lorentz indices associated with the ends of the eikonal line. The other choice is the Lorentz indices associated with the polarization vectors of the cut gluon line. We will refer to those expressions as the Born tensors.
The Born tensor with Lorentz indices associated with the eikonal line is where l µ and T µν are They satisfy l.n = 0 and T µν n ν = 0. Upon contracting A µν eikonal with −g µν , we recover the Born squared amplitude in Eq. (18): The Born tensor with Lorentz indices associated with the final-state gluon is where the tensors are Their coefficients are The argument z of C i (z, p.q) is the longitudinal momentum fraction in Eq. (26), which depends on p and the arbitrary light-like vector q. The tensors in Eqs. (31) satisfy T µν i q ν = 0 because of gauge invariance. Upon contracting A µν gluon with −g µν , we recover the Born squared amplitude in Eq. (18):

III. REAL NLO CORRECTIONS
The real NLO corrections to the perturbative fragmentation function for g → QQ, with the QQ pair in a color-singlet 1 S 0 state, come from cut diagrams with two real partons in the final state. The two partons can be two gluons or a light quark-antiquark pair (qq). Cut diagrams with two real gluons can be obtained from the four LO cut diagrams with a single real gluon, such as the diagram in Figure 1, by adding a gluon line that crosses the cut and runs from any of the 6 colored lines on the left side of the cut to any of the 6 colored lines on the right side of the cut. The additional gluon line can also be attached to the operator vertex, with the fragmenting gluon attached to the eikonal line. The cut diagrams with a light qq pair can be obtained from the four LO cut diagrams by replacing the real gluon line that crosses the cut by a virtual gluon that produces a qq pair that crosses the cut.

A. Subtraction procedure
Each of the cut diagrams involves an integral over the phase space of the two real partons in the final state. The integrals diverge in several phase-space regions, yielding poles of both infrared (IR) and ultraviolet (UV) nature. Five (overlapping) boundaries in the phasespace can be associated with the singular behaviour of the integrand. The boundaries can be defined in terms of Lorentz invariants. We denote the momenta of both the Q andQ by p and the momenta of the final-state partons (which can be gluons or a light quark and antiquark) by q 1 and q 2 . The invariant mass of the four particles in the final state is s = (2p + q 1 + q 2 ) 2 . The boundaries of the singular regions are represented in Figure 2: 1. the integration up to the boundary (2p + q 1 ) 2 /s = 0 yields a UV pole, 2. the integration up to the boundary (2p + q 2 ) 2 /s = 0 yields a UV pole, 3. the integration up to the boundary (q 1 + q 2 ) 2 /m 2 = 0 yields an IR pole, 4. the integration up to the boundary q 2 .n/(q 1 .n + q 2 .n) = 0 yields an IR pole, 5. the integration up to the boundary q 1 .n/(q 1 .n + q 2 .n) = 0 yields an IR pole.
The phase-space regions X i,j connecting two of the above boundaries are associated with double poles, which can either be of pure infrared nature (in the case of X 3,4 and X 3,5 ) or of mixed nature (in the case of X 1,4 and X 2,5 ).
The real NLO contribution to the fragmentation function can be expressed as where N CS is the Collins-Soper prefactor in Eq. (5), dφ real is the product of the differential phase space for final-state partons with momenta q 1 and q 2 and the factor 2πδ(K.n − (2p + q 1 + q 2 ).n) from the cut through the eikonal line, and A real is the squared amplitude. Our strategy to extract the poles in the expression in Eq. (34) is to design a subtraction term D i for each of the five singular regions in Figure 2 whose integral over that region has poles that match those of the integral of A real . The contribution to the fragmentation function in Eq. (34) can be expressed as S u b t r a c t i o n The subtraction terms D i are designed so that the integral in the first term is finite and can be evaluated in D = 4 dimensions. The integrals in the second term are evaluated in D = 4 −2ǫ dimensions, so the UV and IR divergences appear as poles in ǫ. The construction of the subtraction terms D i is described in Sections III B and III D, where we follow closely the subtraction procedure introduced by Catani and Seymour [19]. The analytic integration of the subtraction terms to obtain the poles in ǫ is described in Sections III C and III E, where we again follow closely the procedure introduced in Ref. [19].

B. Subtractions for UV and mixed poles
The UV poles in the real NLO contribution to the fragmentation function are matched by the integrals of the subtraction terms D 1 associated with the limit (2p + q 1 ) 2 /s → 0 and D 2 associated with the limit (2p + q 2 ) 2 /s → 0. The total invariant mass s, the invariant mass s i for the system consisting of the QQ pair and the parton of momentum q i , and the longitudinal momentum fraction y i for that system are Our subtraction term D i associated with the limit s i /s → 0 includes a factor of A µν eikonal (p, q i ), where A µν eikonal is the Born tensor defined in Eq. (27) whose Lorentz indices µ and ν are associated with the eikonal line. The factor A µν eikonal can be interpreted as arising from the fragmentation of a gluon with longitudinal momentum y i K.n into a QQ pair with longitudinal momentum zK.n via the radiation of a gluon of momentum q i .
The subtraction terms D 1 and D 2 are given by where the kernel V UV µν (y i , l i ) is defined by The 4-vectors l 1 and l 2 appearing as the second argument of V UV µν (y i , l i ) in Eq. (37) are defined by These 4-vectors are orthogonal to n: l i .n = 0.

C. Integrals with UV and mixed poles
Explicit expressions for the poles in the integral of the subtraction term D i (i = 1, 2) can be obtained by carrying out the integration over the (3 − 2ǫ)-dimensional slice associated with a fixed value of s i = (2p + q i ) 2 . A convenient decomposition of the phase-space measure is derived in Appendix A: The prefactor N Born (p, q i ) is defined in Eq. (23). The factor dφ Born (p, q i ), which is differential in s i , is defined in Eq. (24). The measure dφ (i) for integration over the slice with fixed s i is where dΩ ⊥ is the transverse angular measure whose integral is 2π 1−ǫ /Γ(1−ǫ). The differential variables s, s i , and y i are defined as functions of p, q 1 , and q 2 in Eqs. (36). The range of y i is from z to 1. The range of s i is from 4m 2 /(z/y i ) to ∞, and the range of s is from s i /y i to ∞.
To carry out the integration over the transverse angles in Ω ⊥ , we observe that the 4vectors l µ 1 and l µ 2 defined in Eq. (39) are orthogonal to n µ , so Lorentz invariance implies where A and B are functions of s i , y i and u. Because of gauge invariance, the amplitude A µν eikonal is orthogonal to n µ and n ν , so that only the term A(−g µν ) survives after contracting the tensor on the right side of Eq. (42) with A µν eikonal . We can determine the coefficient A by contracting both sides of Eq. (42) by g µν : After integrating over the angles in Ω ⊥ , one can make the replacement whereP (real) gg (y) is the real-gluon contribution to the Altarelli-Parisi splitting function for g → g without any regularization of the pole at y = 1: The contraction of −g µν in Eq. (44) with the Born tensor A µν eikonal (p, q i ) gives the Born squared amplitude A Born (p, q i ) obtained from Eq. (18) by replacing q by q i . The UV pole can be made explicit by integrating analytically over the variable s: where NdφA Born (p, q i ) is the LO differential fragmentation function obtained from Eq. (25) by replacing q by q i . The variables s and z in Eq. (25) are replaced by s i and z/y i . The IR pole associated with the y i = 1 endpoint can be extracted by applying the plus prescription: where P (real) gg (y) is the real-gluon contribution to the Altarelli-Parisi splitting function for g → g: With the use of Eq. (47), the expression in Eq. (46) can be expressed as the sum of a term with a double pole, a term with a single pole, and a finite remainder: where the functions I n (z) are g→QQ (z) is defined in Eq. (19). The functions D log (z) and D log 2 (z) are defined by D log (z) = NdφA Born (p, q) log(s/4m 2 ), (51a) where the measure NdφA Born (p, q), which is differential in s = (2p + q) 2 , is given in Eq. (25). These functions appear in Eqs. (50b) and (50c) with argument z/y. In Eqs. (50a) and (50b), there are terms of order ǫ 0 from the expansion of D (LO) g→QQ (z) to order ǫ 2 and the expansion of D log (z) to order ǫ. These expansions are not actually needed, because the canceling poles in ǫ will also be expressed in terms of the functions D (LO) g→QQ (z) and D log (z).

D. Subtractions for IR poles
The IR poles in the real NLO contribution to the fragmentation function are matched by the subtraction terms D 3 , D 4 , and D 5 associated with the limits q 1 .q 2 → 0, q 2 .n → 0, and q 1 .n → 0, respectively. The expressions for these subtraction terms can be made more compact by introducing a light-like 4-vectorq that has the same longitudinal momentum as It satisfiesq 2 = 0 andq.n = (q 1 + q 2 ).n. It is also convenient to introduce variabless, u, and λ defined bys Our subtraction terms D 3 , D 4 , and D 5 include a factor of A µν gluon (p,q), where A µν gluon is the Born tensor defined in Eq. (30) whose Lorentz indices µ and ν are associated with the finalstate gluon. The factor A µν gluon can be interpreted as arising from the fragmentation of a gluon with longitudinal momentum K.n into a QQ pair with longitudinal momentum zK.n via the radiation of a gluon of momentumq.
The subtraction terms D 4 and D 5 are defined by where the kernels V IR(i) (p, q 1 , q 2 ) are The subtraction term D 3 matches the IR poles originating from collinear partons in the final state. It can be expressed in the form The kernel V gg µν + V qq µν has been split into two terms associated with collinear gluons and collinear quarks. It is convenient to introduce a 4-vectorq(u) whose components arȇ The kernels associated with collinear gluons and collinear quarks are where T F = 1 2 is the trace of the square of a generator for the fundamental representation.

E. Integrals with IR poles
Explicit expressions for the poles in the subtraction terms D 3 , D 4 , and D 5 can be obtained by carrying out the phase-space integration over the (3 − 2ǫ)-dimensional slice associated with a fixed value ofs = (2p +q) 2 . A convenient decomposition of the phase-space measure is derived in Appendix A: N CS dφ real (p, q 1 , q 2 ) = N Born (p,q)dφ Born (p,q) dφ(p, q 1 , q 2 ). (59) The prefactor N Born (p,q), which is defined in Eq. (23), coincides with N CS . The factor dφ Born (p,q), which is differential ins, is defined in Eq. (24). The measure dφ for integration over the slice with fixeds is where dΩ ⊥ is the transverse angular measure. The differential variabless, u, and λ are defined as functions of p, q 1 , and q 2 in Eqs. (53). The range ofs is from 4m 2 /z to ∞. The range of λ is from 0 to ∞, and the range of u is from 0 to 1. In the expressions for D 4 and D 5 in Eq. (54), the contraction of −g µν with A µν gluon (p,q) gives the Born squared amplitude A Born (p,q) obtained from Eq. (18) by replacing q byq. In the expression for D 3 in Eq. (56), the Born tensor A µν gluon (p,q) is contracted with the tensors V gg µν and V qq µν defined in Eq. (58). A factor of A Born (p,q) appears only after integrating over the transverse angles in Ω ⊥ . To carry out that integration, we observe that the 4-vectorq(u) defined in Eq. (57) is orthogonal toq, so Lorentz invariance implies where the coefficients C and D are functions of q 1 .q 2 and u. Because of gauge invariance, the Born tensor A µν gluon (p,q) is orthogonal toq µ andq ν , so only the term C(−g µν ) survives after contracting the tensor on the left side of Eq. (61) with A µν gluon (p,q). We can determine the coefficient C by contracting both sides of Eq. (61) by g µν : The IR poles can be made explicit by integrating over the variables u and λ. The integral over λ can be evaluated analytically, and it gives a pole in ǫ. The integral over u gives a second pole in ǫ. After isolating the term that gives the pole, the integrand can be expanded in powers of ǫ and then integrated over u. The resulting expressions for the integrals of D 4 and D 5 are the same: where NdφA Born (p,q) is the LO differential fragmentation function obtained from Eq. (25) by replacing q byq. The function V(p,q) includes all the poles in ǫ: After integrating overs, the integral of D 4 + D 5 reduces to i=4,5 The LO fragmentation function D where the measure NdφA Born , which is differential in s = (2p + q) 2 , is given in Eq. (25). In Eq. (65), there are terms of order ǫ 0 from the expansion of D (LO) g→QQ (z) to order ǫ 2 and the expansion of D log (z) to order ǫ. These expansions are not actually needed, because the canceling poles in ǫ will also be expressed in terms of the functions D (LO) g→QQ (z) and D log (z). In the expression for the integral of D 3 , the IR poles appear in a multiplicative constant factor that can be separated into contributions from gluons and quarks: The factors V gg and V qq are integrals over λ and u that can be evaluated analytically. Their Laurent expansions to order ǫ 0 are After integrating overs, the integral of D 3 reduces to The poles in ǫ in V gg + V qq will be cancelled by another term proportional to D (LO) g→QQ (z), leaving only the terms of order ǫ 0 .

IV. VIRTUAL NLO CORRECTIONS
The virtual NLO corrections to the perturbative fragmentation function for g → QQ, with the QQ pair in a color-singlet 1 S 0 state, come from cut diagrams with one loop on either the right side or the left side of the cut. Loop diagrams on one side of the cut can be obtained from the LO diagrams by adding a gluon line connecting any pair of the 6 colored lines, by adding a loop correction to the propagator of the fragmenting gluon, or by adding a loop correction to the propagator of the virtual heavy quark. There are additional loop diagrams in which the heavy quark line is attached to the eikonal line by both the fragmenting gluon that attaches to the end of the eikonal line and by a second gluon line, with the gluon that crosses the cut attached to either the fragmenting gluon or the eikonal line.
As in the LO cut diagrams, we denote the momenta of both the Q andQ by p and the momentum of the final-state gluon by q. We denote the loop momentum by l. The sum of the virtual one-loop cut diagrams at order α 3 s defines a function A virtual (p, q, l). The virtual NLO contribution to the fragmentation function can be expressed as where N CS is the Collins-Soper prefactor in Eq. (5) and dφ Born is the phase-space measure in Eq. (15). By means of standard tensor reduction techniques, the integral over the loop momentum l in Eq. (70) can be reduced to a sum of one-loop scalar integrals whose numerators are simply 1. Our procedure to apply this reduction is implemented with the use of the Mathematica package FeynCalc [20]. The denominators of the scalar integrals come from Feynman propagators with mass m, massless Feynman propagators, and eikonal propagators of the form i/[(l + P ).n + iǫ], where P is a linear combination of p and q. A product of eikonal propagators can be reduced algebraically to a linear combination of single eikonal propagators. A scalar integral with only Feynman propagators is a function of the invariant mass s = (2p + q) 2 . A scalar integral with one eikonal propagator is a function of s and the momentum fraction z = (2p).n/(2p + q).n. We need the Laurent expansion in ǫ = (4 − D)/2 for each scalar integral to order ǫ 0 . Our results for the scalar integrals with only Feynman propagators are in agreement with results available in the literature [21,22]. The Laurent expansions for the scalar integrals with a single eikonal propagator can be evaluated analytically, with the finite terms order of ǫ 0 expressed in terms of dilogarithms. For some of the integrals, the analytic expressions in terms of dilogarithms are very complicated, so the finite terms might as well be expressed in terms of finite integrals that can be evaluated numerically. Our results for the poles in ǫ in the scalar integrals with one eikonal propagator are given in Appendix B.
In all the poles in ǫ from the loop integral in Eq. (70), the Born squared amplitude A Born (p, q) appears as a multiplicative factor. The virtual NLO corrections to the fragmentation function can therefore be expressed as where f pole (p, q) has only poles in ǫ and A finite (p, q) is a finite function of s = (2p + q) 2 and z.
The terms in f pole (p, q) can be organized to make their cancellation against the poles from other contributions of the NLO correction more transparent: There are four terms in Eq. (72) with only UV poles: is the Casimir for the fundamental representation. In Feynman gauge, the terms U g , U Q (s), U gQQ , and U eikonal arise from virtual-gluon propagator corrections, virtual-quark propagator corrections, quark-gluon vertex corrections, and eikonal line corrections, respectively. There is one term in Eq. (72) with mixed UV and IR poles: In Feynman gauge, this term comes from loop correction to the gluon-eikonal vertex. There are two terms in Eq. (72) with only IR poles: The infrared poles originate from loop-momentum configurations in which partons becoming soft and/or collinear. The term S 1 is a soft pole that in Feynman gauge comes from oneloop diagrams obtained from the four LO cut diagrams by exchanging a gluon between the on-shell heavy quarks. All other infrared poles are included in the term S 2 (s, z). The virtual NLO corrections in Eq. (71) can be expressed as where NdφA Born (p, q) is the LO differential fragmentation function in Eq.

V. RENORMALIZATION
Beyond leading order in α s , the fragmentation function D g→η Q depends on a factorization scale µ F and the running coupling constant α s depends on a renormalization scale µ R . If we make those scales explicit, the expansion of the fragmentation function to NLO in Eq. (7) becomes The scales µ F and µ R are introduced through renormalization. The calculation of the fragmentation function is performed in terms of the renormalized fields Ψ r and A r , the renormalized coupling constant g, and the physical mass m of the heavy quark. Their relations with the corresponding bare quantities involve renormalization constants δ 2 , δ 3 , δ g , and δ m : The renormalization of the coupling constant is performed in the MS scheme, whereas the renormalization of the heavy quark mass is performed in the on-shell mass scheme. In the resulting expressions for the renormalization constants δ 2 , δ 3 , δ g , and δ m , it is convenient to pull out a common factor: The rescaled renormalization constantsδ i in the scheme specified above read where b 0 = (11N c − 4T F n f )/6 is the coefficient of −α 2 s /π in the beta function (d/dµ)α s (µ). In the counterterm for g in Eq. (80c), we have allowed for the renormalization scale µ R of α s to be different from the scale µ introduced through dimensional regularization.
The NLO contributions to the fragmentation function from the counterterms for the propagators and vertices in the LO cut diagrams are The terms with the coefficients 2C gQQ , C eikonal , C g , and C Q (s) are associated with the quarkgluon vertices, the eikonal-gluon vertex, the gluon propagator, and the quark propagator in the LO cut diagrams, respectively. The expressions of these coefficients in terms of the rescaled renormalization constantsδ i 's are The final expression for the counterterm contributions to the NLO fragmentation function are The field renormalization constantsδ 2 andδ 3 include single IR poles. They cancel the single IR poles that remain after adding together the real NLO corrections and the virtual NLO corrections. The real NLO corrections have single IR poles proportional to D  81), the linear combination of field renormalization constants is 2δ 2 +δ 3 . The IR pole in 2δ 2 cancels the IR pole in S 1 introduced in Eq. (72). The IR pole inδ 3 cancels the yet-to-be-cancelled single IR poles in the sum of the contributions from S 2 (s, z) introduced in Eq. (72) and from the terms V gg and V qq in Eqs. (67). This completes the cancellation of the IR poles.
The renormalization of the operator defining the fragmentation function also introduces a counterterm. Its expression in the MS scheme reads where P gg (y) is the Altarelli-Parisi splitting function including both real and virtual contributions: We have allowed for the factorization scale µ F to be different from the scale µ introduced through dimensional regularization. When the contribution from the counterterms D

VI. NUMERICAL RESULTS
Once all the poles have been cancelled in the NLO corrections to the fragmentation function, the function d NLO (z, µ R , µ F ) in Eq. (77) can be obtained by adding all the finite parts: • the subtracted real NLO corrections, which are integrated over the phase space of the two final-state partons in 4 dimensions, • the finite parts of the integrated subtraction terms in Eqs. (49), (65), and (69), • the finite parts of the virtual NLO corrections in Eq. (76), which are integrated over the Born phase space, • the finite parts from the renormalization counterterms in Eqs. (83) and (84).
The numerical integrations are performed with the use of the adaptive Monte Carlo integrator Vegas [23]. The NLO fragmentation function proportional to α 2 renormalization and factorization scales, we choose twice the mass of the heavy quark: µ R = µ F = 2m b . The LO term α 2 s d LO (z) increases monotonically from 0 to α 2 s /(36m 3 b ) as z increases from 0 to 1. The NLO term α 3 s d NLO (z) increases from −∞ as z → 0 to a broad maximum at an intermediate value of z, and then decreases to −∞ as z → 1. For µ R = µ F = 2m b , its maximum is 2.9 α 3 s /(36m 3 b ) at z = 0.45. The NLO term is negative near both endpoints, but it cannot be compared to the LO term in these regions, because the d NLO (z) is actually a distribution in z with delta-function contributions at the endpoints z = 0 and z = 1. In the integral of the product of d NLO (z) and a smooth function of z, the endpoint contribution cancels a divergence in the integral up to the endpoint, so that the integral over the endpoint region is well behaved. The NLO term can be compared to the LO term in the intermediate region of z, and it is larger than the LO term in this region. At z = 0.5, the NLO fragmentation function is larger than the LO fragmentation function by a factor of 2.7. The total fragmentation probability obtained by integrating the NLO fragmentation function over z from 0 to 1 is larger than the LO fragmentation probability by a factor 1.89 for the choice of scales µ R = µ F = 2m b . The mean value z of the longitudinal momentum fraction is 2/3 at LO, and it decreases to 0.54 at NLO.
The sensitivity of the LO and NLO fragmentation functions to the renormalization scale µ R is illustrated in Figure 3. The bands are obtained by varying µ R up or down by a factor of 2 around the central value 2m b (with µ F = 2m b ). The NLO band in Figure 3  wider than the LO band except near the endpoint at z = 1. One might have expected the sensitivity to µ R to be decreased by adding NLO corrections, but this is not the case simply because the NLO term in the fragmentation function is larger than the LO term in the central region of z. The ratio of the fragmentation functions at NLO and LO is less sensitive to the renormalization scale. At z = 0.5, the ratio changes from 2.79 to 2.71 to 2.67 as µ R varies from m b to 2m b to 4m b . The ratio of the fragmentation probabilities at NLO and LO changes from 1.75 to 1.89 to 1.99. The mean momentum fraction z increases from 0.49 to 0.54 to 0.57.
The sensitivity of the NLO fragmentation function to the factorization scale µ F is illustrated in Figure 4. The band is obtained by varying µ F up or down by a factor of 2 around the central value 2m b (with µ R = 2m b ). In the central region of z, the width of the band from varying µ F is much narrower than that from varying µ R in Figure 3. The width increases near the endpoints of z at 0 and 1, but the fragmentation function also has canceling endpoint contributions at z = 0 and z = 1. Therefore the increased sensitivity to µ F near the endpoints will not result in a large increase in sensitivity for the integral of the product of the fragmentation function and a smooth function of z. The ratio of the fragmentation probabilities at NLO and LO is more sensitive to µ F than to µ R , ranging from 2.58 to 1.89 to 1.21 as µ F varies from m b to 2m b to 4m b .

VII. SUMMARY
In this paper, we have presented the NLO calculation of the fragmentation function for a gluon into a spin-singlet S-wave quarkonium state η Q at leading order in v. This calculation represents the first NLO result for a fragmentation function into quarkonium that is a nontrivial function of z at LO. We have found that the real NLO correction can be organized in an efficient way by constructing a set of subtraction terms matching the poles in each phase-space boundary leading to singularities. This strategy allows for a transparent organization of both UV and IR poles and their cancellation among the different components of the calculation (real correction, virtual corrections, and counterterms). It also paves the way to automation of the NLO calculation of the fragmentation functions in other NRQCD channels.
We found that the NLO QCD corrections have a dramatic effect on the fragmentation function in the MS renormalization and factorization schemes. The effect on the shape of the fragmentation function is particularly dramatic. Instead of increasing monotonically with z as at LO, the NLO fragmentation function has a broad maximum in the central region of z. In this region, it is about a factor of 3 larger than at LO. As a consequence, the NLO fragmentation function displays strong sensitivity to the renormalization scale. These results suggest that QCD corrections to fragmentation functions could have a significant impact on the production of quarkonium states at large transverse momentum.
The phase-space measure for a single massless parton of momentum q can be expressed as d D−1 q (2π) D−1 2q 0 = 1 2(2π) 3−2ǫ q 1−2ǫ 0 dq 0 | sin θ| 1−2ǫ dθ dΩ ⊥ , where θ is the polar angle, and dΩ ⊥ is the measure for integration over angles in the transverse plane. The total transverse solid angle is Ω ⊥ = 2π 1−ǫ /Γ(1 − ǫ). The light-like 4-vector n defines a spatial direction that can be used as the polar axis that determines the polar angle θ. However it is sometimes convenient to choose the polar axis in the spatial direction of a different light-like 4-vector k. A convenient alternative set of variables from q 0 and θ is the longitudinal momentum q.n and a variable λ defined by λ = 2k.q/k.n.
The phase-space measure can be expressed as This phase-space measure is independent of the overall scales of n and k. In a Lorentz frame where the spatial parts of n and k are back-to-back, q.n = q 0 n 0 (1 − cos θ) and λ = q 0 (1 + cos θ)/n 0 , so Eq. (A4) reduces to the measure in Eq. (A2). A convenient parameterization of the two-parton phase-space measure in Eq. (A1) can be obtained by introducing two light-like 4-vectors that define the polar axes of q 1 and q 2 : a 4-vector k 1 that may depend on p and a 4-vector k 2 that may depend on p and q 1 . We also introduce a momentum fraction u defined by u = q 2 .n (q 1 + q 2 ).n .
After integrating over the longitudinal component and transverse angles of q 1 , the phasespace measure reduces to dφ real (p, q 1 , Phase space for UV and mixed poles. To obtain the phase-space parameterization used to integrate the UV and mixed poles in Section III C, we choose the light-like vectors k 1 and k 2 that specify the polar axes for q 1 and q 2 to be k µ 1 = 2p µ − m 2 p.n n µ , k µ 2 = (2p + q 1 ) µ − (2p + q 1 ) 2 2(2p + q).n n µ .
Phase space for IR poles. To obtain the phase-space parameterization used to integrate the IR poles in Section III E, we first introduce additional integrals over a light-like 4-vectorq that has the same longitudinal component as q 1 + q 2 and over the invariant mass (q 1 + q 2 ) 2 . We do this by multiplying the phase-space measure in Eq. (A1) by 1 in the form 1 = ∞ 0 dt d Dq δ D (q − (q 1 + q 2 ) + [t/2(q 1 + q 2 ).n]n) δ(t − (q 1 + q 2 ) 2 ). (A11) ǫ. The finite terms of order ǫ 0 can also be evaluated analytically in terms of logarithms and dilogarithms of functions of z and r. For some of the loop integrals, the evaluation of the integrals using a computer algebra program, such as Mathematica, gives dilogarithms with many different arguments. The number of different arguments can be reduced by using functional identities for dilogarithms. However, for many of the integrals, the expressions for the finite terms are still sufficiently complicated that we do not present them here.

Integrals with two Feynman propagators
There are 12 independent integrals with two Feynman propagators and an eikonal propagator. These integrals have a UV divergence that yields a single pole in ǫ. They may also have an IR divergence that yields a second pole in ǫ. If there is a double pole in ǫ, the nature of the subleading single pole is ambiguous. It is convenient to express these integrals as the product of a Laurent expansion in ǫ and an overall factor of It is convenient to introduce a compact notation for the Feynman propagators. The propagators for momentum k and masses 0 and m are denoted by D(k, m) = 1 k 2 − m 2 + iε . (B7b)