Gluon fragmentation into quarkonium at next-to-leading order using FKS subtraction

We present the calculation at next-to-leading order (NLO) in αs of the fragmentation function of a gluon into heavy quarkonium in the color-octet spin-singlet S-wave channel. To calculate the real NLO corrections, we adapt a subtraction scheme introduced by Frixione, Kunszt, and Signer. Ultraviolet and infrared divergences in the real NLO corrections are calculated analytically by evaluating the phase-space integrals of the subtraction terms using dimensional regularization. The subtracted phase-space integrals are then evaluated in 4 space-time dimensions. The divergences in the virtual NLO corrections are also calculated analytically. After renormalization, all the divergences cancel. The NLO corrections significantly increase the fragmentation probability for a gluon into the spin-singlet quarkonium states ηc and ηb.


Introduction
The production of a hadron with large transverse momentum p T can be simplified by using a factorization theorem for inclusive hadron production at large p T [1]. It states that the leading power in the expansion of the inclusive cross section in powers of 1/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 i (gluons, quarks, and antiquarks). The symbol "⊗" in eq. (1.1) represents a convolution integral over the longitudinal momentum fraction z of the hadron H relative to the parton. The pQCD cross section dσ for producing the parton i can be expanded in powers of α s (p T ), and it includes convolutions with parton distributions if the colliding particles are hadrons. The nonperturbative factor D i→H (z) is a fragmentation function that gives the probability distribution for z. We refer to eq. (1.1) as the leading-power (LP) factorization formula. It was derived by Collins and Soper in 1981 for the case of a light hadron H at a transverse momentum satisfying p T Λ QCD [1]. The LP factorization formula in eq. (1.1) applies equally well to heavy quarkonium at a transverse momentum satisfying p T m, where m is the mass of the heavy quark. A proof of this factorization theorem that deals with issues specific to heavy quarkonium production was first sketched by Nayak, Qiu, and Sterman in 2005 [2]. The LP factorization formula gives the leading power in the expansion in powers of m/p T . In the case of cross sections summed over quarkonium spins, the corrections are suppressed by m 2 /p 2 T . 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 for cross sections for heavy quarkonium production [3]. It 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, where v is the typical relative velocity of the Q orQ in the quarkonium. The theoretical status of the NRQCD factorization conjecture is discussed in ref. [4]. The predictive power of the LP factorization formula for heavy quarkonium in eq. (1.1) can be increased by applying the NRQCD factorization formula to the fragmentation functions, reducing these nonperturbative functions of z to multiplicative constants. The fragmentation function for the parton i to produce the quarkonium H is expressed as a sum of pQCD fragmentation functions multiplied by NRQCD matrix elements: (1. 2) The sum extends over the color and angular-momentum channels n of a nonrelativistic QQ pair. The pQCD fragmentation functionD i→QQ[n] (z) for producing the QQ pair in the channel n can be expanded in powers of α s (m). The NRQCD matrix element O H n is proportional to the probability for a QQ pair created in the channel n to evolve into a

JHEP01(2019)227
final state that includes the quarkonium H. The nonperturbative constants O H n scale as definite powers of v [3].
The factorization theorem for inclusive production of heavy quarkonium at large p T has been extended to the next-to-leading power (NLP) of m 2 /p 2 T . The NLP factorization theorem was proven diagrammatically by Kang, Qiu, and Sterman [5][6][7], and it was derived using soft collinear effective theory by Fleming, Leibovich, Mehen, and Rothstein [8,9]. In addition to corrections to the terms in the LP factorization formula, the NLP factorization formula has additional terms suppressed by m 2 /p 2 T that are expressed as a sum of pQCD cross sections for producing a pair of collinear partons convolved with double-parton fragmentation functions. In the case of cross sections summed over quarkonium spins, the corrections to the NLP fragmentation formula are suppressed by a power of m 4 /p 4 T . The predictive power of the NLP fragmentation formula can be dramatically increased by applying the NRQCD factorization formula to the double-parton fragmentation functions, reducing these nonperturbative functions to multiplicative constants.
The NLP/NRQCD factorization formula opens the door to dramatic improvements in the accuracy of theoretical predictions for quarkonium production at 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 calculating all the pQCD factors to next-to-leading order (NLO) in α s . The pQCD factors are the cross sections for producing single partons, the cross sections for producing collinear parton pairs, the single-parton fragmentation functions for producing QQ pairs, the double-parton fragmentation functions for producing QQ pairs, and the evolution kernels for both sets of fragmentation functions.
The first fragmentation function for quarkonium production to be calculated to nextto-leading order (NLO) in α s was for gluon fragmentation into QQ in the color-octet 3 S 1 channel [10,11]. This calculation is particularly simple, because the LO fragmentation function is proportional to δ(1 − z). The first NLO calculation of a fragmentation function that at LO is a nontrivial function of z was that for gluon fragmentation into QQ in the color-singlet 1 S 0 channel [12]. It would be useful to have all the phenomenologically relevant fragmentation functions calculated to NLO in α s .
In NLO QCD calculations, the most challenging step is often the calculation of the phase-space integrals from real-gluon emission. A strategy that is often effective is to design subtractions that cancel the infrared divergences in the phase-space integrals, calculate the integrals of the subtraction terms analytically, and then calculate the subtracted phase-space integrals numerically. In the case of fragmentation functions, the phase-space integrals also have ultraviolet divergences. It is therefore necessary to design subtractions that also cancel these ultraviolet divergences. The NLO calculation of ref. [12] was carried out using a subtraction procedure for fragmentation functions that was adapted from the dipole subtraction procedure for parton cross sections introduced by Catani and Seymour [13]. An alternative subtraction scheme for parton cross sections that has some advantages was introduced by Frixione, Kunszt, and Signer (FKS) [14]. The FKS subtraction method has been used in the automation of next-to-leading order computations of parton cross sections in QCD [15].

JHEP01(2019)227
In this paper, we adapt the FKS subtraction method to the NLO calculation of fragmentation functions. We illustrate the method by applying it to gluon fragmentation into a QQ pair in the color-octet 1 S 0 channel. This fragmentation function is of phenomenological importance for the production of J P C = 0 −+ quarkonium states, such as the η b or η c . The color-singlet 1 S 0 channel is leading order in v. The color-octet 1 S 0 channel is one of three color-octet channels suppressed by only v 4 . In the case of production of the η c at large p T at the Large Hadron Collider, the color-octet 1 S 0 channel is numerically the most important of the three [16,17].
The outline of our paper is as follows. In section 2, we present the LO fragmentation function for gluon fragmentation into QQ in the color-octet 1 S 0 channel and define some quantities that are useful in the NLO calculation. In section 3, we introduce the FKS subtraction terms that cancel all the ultraviolet and infrared divergences in the real NLO corrections. The subtracted phase-space integrals are calculated in 4 dimensions, and their insensitivity to the cut parameters in the FKS subtractions is verified. In section 4, we present analytic results for the ultraviolet and infrared divergences from loop integrals in the virtual NLO corrections. In section 5, we verify that all the divergences from phase-space integrals and from loop integrals are canceled by renormalization of the parameters of QCD and by renormalization of the operator whose matrix element defines the fragmentation function. Some numerical illustrations of our results are presented in section 6. Our results are summarized in section 7. In appendix A, we derive parametrizations of massless two-parton phase-space integrals that are used to integrate the subtractions terms for the real NLO corrections. In appendix B, we calculate the pole terms in the dimensionally regularized phase-space integrals of the subtraction terms analytically. In appendix C, we show how MadGraph5 [18] can be used to generate helicity amplitudes for cut diagrams with a heavy-quark pair and two light partons in the final state.

Leading-order fragmentation function
In this section, we present the perturbative fragmentation function for g → QQ, with the QQ pair in a color-octet 1 S 0 state, at leading order in α s . We also introduce some related expressions that are useful in the calculation of the real radiative corrections at next-to-leading order in α s .

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 with 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 on 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 JHEP01(2019)227 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]: they are summarized in ref. [12]. 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 gluon line attached to the operator and the remainder flows through the eikonal line. The particle into which the gluon is fragmenting has a specified 4-momentum. In the case of fragmentation of a gluon into a QQ pair with zero relative momentum, it is convenient to express the 4momentum of the QQ pair as 2p. The longitudinal momentum fraction z of the QQ pair is: The operator at the left end of the eikonal line is labelled by a Lorentz index µ and a color index a. The Feynman rule for the operator vertex is where q, λ and c are the 4-momentum, Lorentz index, and color index of the gluon line attached to the vertex. The operator at the right end of the eikonal line is labelled by a Lorentz index ν and a color index b. The fragmentation function is the sum of all cut diagrams contracted with −g µν and δ ab and multiplied by the Collins-Soper prefactor [1]: where N c = 3 is the number of colors of a quark and D = 4−2 is the space-time dimension.

LO fragmentation function
The NRQCD factorization formula in eq. (1.2) for the fragmentation function D i→H (z) for producing the quarkonium state H expresses it as a sum of pQCD fragmentation functions multiplied by NRQCD matrix elements. The NRQCD matrix elements O H n scale as definite powers of the relative velocity v of the heavy quark in the quarkonium. The pQCD fragmentation functionsD i→QQ[n] (z) can be calculated as power series in α s (m), where m is the heavy-quark mass. The fragmentation function for a gluon into a 0 −+ quarkonium state η Q , including the color-singlet 1 S 0 and color-octet 1 S 0 terms explicitly, has the form NLO (z) + . . .
The color-singlet NRQCD matrix element O 1 ( 1 S 0 ) η Q is leading order in v, and it can be expressed in terms of the wavefunction at the origin for the η Q . The matrix element O 8 ( 1 S 0 ) η Q is one of three color-octet matrix elements that are suppressed by only v 4 . It is related by heavy-quark spin symmetry to the NRQCD matrix element O 8 ( 3 S 1 ) ψ Q for a 1 −− quarkonium state ψ Q , which can be determined phenomenologically by fitting cross sections for production of ψ Q . In eq. (2.4), the PQCD fragmentation functions multiplying have been expanded to next-to-leading order (NLO) in α s . The color-singlet 1 S 0 fragmentation function was calculated at leading order (LO) in α s by Braaten and Yuan in 1993 [19]: (2.5) The NLO term D NLO (z) in this fragmentation function was calculated in ref. [12]. Our goal is to calculate the color-octet 1 S 0 fragmentation function to NLO.
The pQCD fragmentation functions in eq. (2.4) can be determined from perturbative QCD calculations of the fragmentation function for producing a QQ pair. Fragmentation functions for producing QQ in a 1 S 0 state can be determined most easily by taking the QQ pair to be in a spin-singlet state with zero relative momentum. To determine the colorsinglet fragmentation function, the Q andQ are projected onto the color-singlet state QQ 1 by contracting their color indices i and j with δ ij / √ N c . To determine the color-octet fragmentation function, the Q andQ are projected onto the color-octet state QQ 8 with color index a by contracting their color indices i and j with √ 2 T a ij . Given the normalizations of the NRQCD operators defined in ref. [3], the perturbative approximations to the NRQCD matrix elements are where N c = 3. On the left side of eq. (2.6b), there is an implied sum over the N 2 c −1 colors of QQ 8 . If dimensional regularization is used to regularize ultraviolet and infrared divergences, JHEP01(2019)227 these matrix elements have no NLO corrections. By dividing the perturbatively calculated fragmentation function D g→QQ 8 (z) by the perturbative matrix element in eq. (2.6b), we obtain the fragmentation function multiplying O 8 ( 1 S 0 ) η Q in eq. (2.4).

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 obtained by interchanging the two gluon vertices on the left side of the cut and interchanging the two gluon vertices on the right side of the cut. 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 in figure 1 is given in eq. (2.10) of ref. [12]. The QQ pair is projected onto the colorsinglet state QQ 1 in eq. (2.11) of ref. [12]. If the amplitude is instead projected onto the color-octet state QQ 8 with color index b, the net effect is the replacement where a is the color index of the gluon-field-strength operator and c is the color index of the final-state gluon. After squaring the amplitudes, contracting the color indices for the operators, and summing over the colors of the gluon and the QQ pair, the net effect is the substitution This procedure projects the QQ pair onto the color-octet states QQ 8 . Thus the LO fragmentation function D g→QQ 8 (z) differs from the LO fragmentation function D g→QQ 1 (z) just by the multiplicative color factor (N 2 c − 4)/2. The LO fragmentation function D (2.9) In the calculation of the NLO fragmentation function for g → QQ 8 , it is useful to have the LO fragmentation function for g → QQ 8 expressed as an integral over the gluon phase space in D = 4 − 2 dimensions: 10) where N CS is the Collins-Soper prefactor in eq. (2.3). The Born phase-space measure dφ Born is the product of the differential phase space for the final-state gluon of momentum q and JHEP01(2019)227 a factor 2πδ(K.n − (2p + q).n) from the cut through the eikonal line. It can be reduced to a single differential in the invariant mass s of the QQg system: where s and z expressed as functions of p and q are There is an implied Heavyside theta function that imposes the constraint s > 4m 2 /z. The Born squared amplitude A Born is obtained by multiplying the right side of eq. (2.15) in ref. [12] by the color factor (N 2 c − 4)/2: (2.14) Dividing by the perturbative NRQCD matrix element in eq. (2.6b), we obtain the LO fragmentation function D LO (z) in eq. (2.9) multiplied by α 2 s .

Born tensors
To facilitate the calculation of the real NLO corrections to the fragmentation function, it is convenient to generalize the integration measure dφ Born (p, q) for the LO fragmentation function in eq. (2.11) by allowing q to be a more general light-like 4-vector. It could be the momentum q 1 or q 2 of a massless final-state parton, or it could be another light-like 4-vector constructed from q 1 and q 2 . The variables s = (2p + q) 2 and z = (2p.n)/(2p + q).n defined in eqs. (2.12) can be regarded as functions of this more general light-like 4-vector q.

JHEP01(2019)227
The Born phase-space measure in eq. (2.11), with the factor of 1/(2p + q).n replaced by 1/(2p + q 1 + q 2 ).n and with its coefficient expressed as a function of s and z, will be denoted by dφ Born (s, z): Similarly, the Born squared amplitude A Born (p, q) in eq. (2.13), with the factor of [(2p + q).n] 2 replaced by [(2p + q 1 + q 2 ).n] 2 and with its coefficient expressed as a function of s and z, will be denoted by A Born (s, z): (2.17) The product of N Born , dφ Born , and A Born in eqs. (2.15), (2.16), and (2.17) depends only on s and z. We introduce a more concise notation for this product: 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. If the weight function is simply 1, the integral is the LO fragmentation function for g → QQ 8 in D dimensions defined in eq. (2.10): 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 real NLO corrections point-by-point in the phase space. 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 obtained by multiplying the right side of eq. (2.24) in ref. [12] by the color factor (N 2 c − 4)/2:

JHEP01(2019)227
where l µ and T µν are The tensor A µν eikonal is orthogonal to n µ and n ν . Its contraction with −g µν is where A Born is given in eq. (2.17).
The Born tensor with Lorentz indices associated with the final-state gluon is obtained by multiplying the right side of eq. (2.27) in ref. [12] by the color factor (N 2 c − 4)/2: where the tensors are Their coefficients are where z is the momentum fraction in eq. (2.12b). The tensor A µν gluon is orthogonal to q µ and q ν . Its contraction with −g µν is where A Born is given in eq. (2.17).

Real NLO corrections
The real NLO corrections to the perturbative fragmentation function for g → QQ, with the QQ pair in a color-octet 1 S 0 state, come from cut diagrams with two real partons in the final state. The two partons can be two gluons (gg) 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, in which case the fragmenting gluon is 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. Each of the cut diagrams involves an integral over the phase space of the two real partons in the final state. We denote the equal momenta of the Q andQ by p and the momenta of the final-state partons (which can be two gluons or a light quark and antiquark) by q 1 and q 2 . The real NLO contribution to the fragmentation function can be expressed as where N CS is the Collins-Soper prefactor in eq. (2.3) and 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. The longitudinal momentum fraction of the QQ pair is In the integrand of eq. (3.1), A real is the squared amplitude from cut diagrams with two gluons crossing the cut, and A (qq) real is the squared amplitude from cut diagrams with a light quark and antiquark crossing the cut. The factor 1 2 multiplying A (gg) real in eq. (3.1) compensates for overcounting the states by integrating over the entire phase space of the two identical gluons.

Anatomy of the poles
The phase-space integrals in eq. (3.1) diverge in several regions, yielding poles in = (4 − D)/2 of both infrared (IR) and ultraviolet (UV) nature. The nature of the IR poles can be soft or collinear.

Soft infrared limits
Soft singularities only show up in the squared amplitude A (gg) real from final-state gluons since a soft (anti)quark does not yield any pole. First let us consider the limit q 2 → 0 (soft limit for the gluon 2). The tree-level amplitude is proportional to the color factor d acb in eq. (2.7), where a, b, and c are the color indices of the eikonal line, the QQ pair, and JHEP01(2019)227 the final-state gluon, respectively. In the eikonal approximation associated with the soft limit q 2 → 0, the amplitude for the emission of the soft-gluon is obtained from the tree amplitude by replacing d acb with where µ and d are the Lorentz index and color index of the gluon with soft momentum q 2 . The first term inside the parentheses is an eikonal factor built upon the momentum of the heavy-quark pair. Such a term arises only when the heavy-quark pair is projected onto a color-octet state. Squaring the amplitude and summing over color and Lorentz indices, one readily obtains the following approximation for the squared amplitude in the soft limit q 2 → 0: where A Born is the Born squared amplitude defined in eq. (2.13). The analogous result in the soft limit q 1 → 0 reads

Collinear infrared limits
Collinear singularities show up in both the squared amplitudes A (gg) real and A (qq) real from the kinematic region in which the light partons crossing the cut are collinear. In order to describe the collinear limit, it is convenient to define the light-like four-vector which satisfiesq 2 = 0 andq.n = (q 1 + q 2 ).n. In the collinear limit, the four-vectorq coincides with q 1 + q 2 , and the expression for the squared amplitude factorizes over the Born amplitude up to spin correlations. In order to account for these spin correlation effects in collinear splittings, it is convenient to define a 4-vectorq by This 4-vector is orthogonal toq:q.q = 0. The 4-momenta of the two partons can be expressed as In the case of A (gg) real , the factorisation formula when the momenta of the two gluons crossing the cut are nearly collinear reads

JHEP01(2019)227
where the Born tensor A µν gluon is defined in eq. (2.23). The tensor P (gg) µν is (3.10) In the case of A (qq) real , the factorisation formula when the momenta of the light quark and antiquark crossing the cut are nearly collinear reads The tensor P where T F = 1 2 is the trace of the square of a generator for the fundamental representation and n f is the number of light flavours.
In the squared amplitude A (gg) real , additional collinear singularities arise from kinematic regions in which one gluon crossing the cut is collinear to the eikonal line. When the gluon momentum q i is collinear to the four-vector n, one has where Q = 2p + q 1 + q 2 is the sum of the momenta of the particles crossing the cut.

Ultraviolet limits
Ultraviolet singularities in the squared amplitude A (gg) real arise from the kinemetic region in which the invariant mass s = (2p + q 1 + q 2 ) 2 of all the final-state particles goes to ∞. The factorisation formula in each of the two UV limits s (2p + q j ) 2 for j = 1, 2 also holds up to spin correlations. In order to account for these spin correlation effects, it is convenient to introduce the 4-vectors l j defined by These 4-vectors are orthogonal to n: l j .n = 0. The 4-vector l µ 1 is the component of q µ 2 orthogonal to n. It is also convenient to define the longitudinal momentum fraction y j of the system consisting of the QQ pair and the parton of momentum q j : The factorisation formula includes a factor of A µν eikonal (p, q j ), where A µν eikonal is the Born tensor defined in eq. (2.20) 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

JHEP01(2019)227
with longitudinal momentum y j K.n into a QQ pair with longitudinal momentum zK.n via the radiation of a gluon of momentum q j . In the limit s (2p + q j ) 2 , one has The tensor P eik µν is defined by

Subtraction scheme and kinematics
In order to extract the poles resulting from phase-space integrals in the singular regions outlined in the previous section, we adopt in this work a subtraction scheme that is inspired by the formalism introduced by Frixione, Kunszt, and Signer (FKS) [14]. In the FKS formalism for parton cross sections [14], the phase space is partitioned into kinematic regions that in any singular limit give either (a) a finite contribution or (b) a soft IR singularity or (c) a collinear IR singularity or (d) the product of a soft IR singularity and a collinear IR singularity. In the extension of the FKS formalism to fragmentation functions, a kinematic region in a singular limit may also give (e) a UV singularity or (f) the product of a UV singularity and a collinear IR singularity. For the (gg) term in eq. (3.1), the partition of the phase space can be implemented by multiplying real by a partition of unity with four terms: where the weight functions S i,j are scalar functions of p, q 1 , q 2 , and n. Each S i,j acts as a damping factor for some singular phase-space regions. The second index j is associated with the parton of momentum q j . For the (qq) term in eq. (3.1), no partition of unity is required since there is only one kinematic region yielding a pole, namely the singular limit of a collinear qq pair. Our strategy to extract the poles in the expression in eq. (3.1) is to design a subtraction term T (gg) i,j for each of the four terms 1 2 A (gg) real S i,j in the integrand and a subtraction term T (qq) for the term A (qq) real . Hence the fragmentation function in eq. (3.1) can be decomposed as where the dependence on the momenta p, q 1 , q 2 is implicit in each function in the integrands. The subtraction terms T (gg) i,j and T (qq) are designed so that the first and second integrals on the right side of eq. (3.19) are finite and can be evaluated in D = 4 dimensions. The third and fourth integrals are evaluated in D = 4 − 2 dimensions, so the UV and IR divergences appear as poles in . The first and second integrals in eq. (3.19) are decomposed into sums of terms such that the subtracted singularity in each term is associated with a very simple JHEP01(2019)227 kinematic region. Such a decomposition is particularly suited for multichannel Monte Carlo integration, as it naturally separates the different channels by making use of the partition of unity in eq. (3.18).
Unlike the construction in the FKS formalism, Lorentz invariance is explicitly manifest in our construction of the subtraction terms, i.e. we do not specify any specific frame. Singular regions are characterised by Lorentz-invariant quantities. We define a dimensionless variable λ that vanishes at the boundary that gives IR singularities associated with two collinear partons and/or one soft parton: We define a dimensionless variable u 1 (u 2 ) that vanishes at the boundary that gives IR singularities when gluon 2 (gluon 1) is soft and/or collinear to the eikonal line: They satisfy u 1 + u 2 = 1. (The mismatch between the indices of the variable u j that approaches 0 and the momentum q 3−j that becomes soft or collinear simplifies the expressions for subtraction terms.) The total invariant mass s of the final-state particles and the invariant mass s j of the heavy-quark pair and the gluon j are We define dimensionless variables ζ j that vanish at the boundaries that give UV singularities when s is much larger than s j : where the momentum fraction y j is defined by eq. (3.15). The boundaries of the singular regions are represented in figure 2.
The subtraction terms T (gg) i,j and T (qq) introduced in eq. (3.19) will be defined using cutoff variables, so that the subtraction is applied only in the vicinity of the singular phasespace boundaries. The variables on which these cutoffs apply are selected in such a way to make the integration of the poles in in the subtraction integrals tractable.

Partition of unity
To define the weight functions S i,j in the partition of unity in eq. (3.18), we first introduce functions w i,j : Representation of the singular regions for the integration of the real emission amplitude. Each line represents a specific limit, which is specified in terms of Lorentz invariants in the inner part of the figure. A line represents a phase-space boundary that gives a single pole in . A dot connecting two lines represents a phase-space region leading to a double pole in . For each dot, there is a specific subtraction term T (gg) i,j appearing in eq. (3.19) whose integral includes the double pole. The subtraction of the single pole associated with a line connecting two dots is shared among the two subtraction terms associated with these dots.
The sum of these four functions is The weight functions in the partition of unity are defined by At the singular boundaries of phase space, one or two of the variables λ, u 1 , u 2 , ζ 1 , and ζ 2 vanishes. In the singular limits, some of the weight functions S i,j vanish and the others simplify. In the collinear IR limit λ → 0, only S 2,1 and S 2,2 are nonzero (and they add up to 1): (3.27) In the soft IR limit u j → 0, only S 1,j and S 2,j are nonzero (and they add up to 1):

JHEP01(2019)227
In the UV limit ζ j → 0, only S 1,j is nonzero: In the construction of the subtraction terms T (gg) i,j in eq. (3.19), we use Heavyside theta functions in the variables λ, u j , and ζ j associated with the singular boundaries: (3.30) In the construction of the subtraction terms T (gg) 2,j in eq. (3.19), we also use Heavyside theta functions in variables δ and δ j defined by , (3.31a) The theta functions are In the construction of the subtraction terms T (qq) in eq. (3.19), we use the Heavyside theta function θ (λ) in the variable λ.

Subtraction terms for
The phase-space boundaries that give singularities in the integral of 1 2 A (gg) real S 1,j are u j = 0 and ζ j = 0. The singular piece can be expressed as (3.33) The remainder obtained by subtracting T (gg) 1,j from 1 2 A (gg) real S 1,j can be integrated numerically in D = 4 dimensions, because the θ (u j ) term subtracts the singularity when u j → 0, the θ (ζ j ) term subtracts the singularity when ζ j → 0, and the θ (ζ j ) θ (u j ) term adds back the double singularity that is over-subtracted by the other two terms.
The weight function S (u j ) 1,j in the first subtraction term in eq. (3.33) is where λ j is defined by (3.35) In the soft limit q 3−j → 0, this weight function approaches S 1,j , whose limiting behavior is given in eq. (3.28a). The subtraction term D (u j ) 1 is constructed using the Born squared amplitude A Born in eq. (2.13) with s replaced by s j /y j : The subtraction term D (ζ j ) in eq. (3.33) is constructed using the Born tensor A µν eikonal in eq. (2.20) with Lorentz indices associated with the eikonal line: where the 4-vector l j is defined in eqs. (3.14). The over-subtraction term D (ζ j ,u j ) in eq. (3.33) is constructed using the Born squared amplitude A Born (s j /y j , z): The poles in the integral of D (ζ j ) are evaluated analytically in section B.2 of appendix B. The over-subtraction integral summed over j is given in eqs. (B.11) and (B.12).

Subtraction terms for T
The phase-space boundaries that give singularities in the integral of 1 2 A (gg) real S 2,j are u j = 0 and λ = 0. The singular piece can be expressed as (3.40) The remainder obtained by subtracting T (gg) 2,j from 1 2 A (gg) real S 2,j can be integrated numerically in D = 4 dimensions, because the θ (u j ) term subtracts the singularity when u j → 0, the θ (λ) term subtracts the singularity when λ → 0, the θ (u j ) θ (δ j ) term subtracts additional singularities when both u j → 0 and δ j → 0, and the last two terms add back singularities that are over-subtracted by the other terms.
The weight function S (u j ) 2,j in the first subtraction term in eq. (3.40) is wheres = (2p +q) 2 . In the soft limit q 3−j → 0, this weight function approaches S 2,j , whose limiting behavior is given in eq. (3.28b). The subtraction term D (u j ) 2 is constructed using the Born squared amplitude A Born in eq. (2.13) with s replaced bys = (2p +q) 2 : z). (3.43) Since u 1 + u 2 = 1, they satisfy S In the collinear IR limit λ → 0, this weight function approaches S 2,j , whose limiting behavior is given in eq. (3.27). The subtraction term D (λ) is constructed from the Born tensor A µν gluon in eq. (2.23) with Lorentz indices associated with the final-state gluon: where the tensor P The subtraction term D (u j ,δ j ) is constructed using the Born squared amplitude A Born (s j , z), where the variables j is defined by with s j = (2p + q j ) 2 and u j defined in eq. (3.21). The subtraction term is where w j is defined by The over-subtraction terms D (u j ,δ) and D (λ,u j ) in eq. (3.40) are the same. They are constructed using the Born squared amplitude A Born (s, z):

Subtraction term for T (qq)
The only phase-space boundary that gives singularities in the integral of A (qq) real is λ = 0. The singular piece can be expressed as (3.50) The remainder obtained by subtracting T (qq) from A (qq) real can be integrated numerically in D = 4 dimensions, because the θ (λ) term subtracts the singularity when λ → 0.
The subtraction term D (qq) is constructed from the Born tensor A µν gluon in eq. (2.23) with Lorentz indices associated with the final-state gluon: The tensor P

Insensitivity to cut parameters
The finite parts of the real NLO corrections to the fragmentation function for g → QQ 8 can be obtained by adding two contributions: • the subtracted real NLO corrections, which are given by the first two integrals on the right side of eq. (3.19). The integrals over the phase space of the two final-state partons in 4 dimensions are evaluated numerically.
• the finite parts of the subtractions for the real NLO corrections, which are given by the last two integrals on the right side of eq. Both contributions depend on the cut parameters u cut , λ cut , δ cut , and ζ cut that define the phase-space regions where the subtractions are applied. The dependence on the cut parameters must cancel between the two contributions. The fact that the overall sum of finite pieces should be independent of the cut parameters can be used as a sanity check for the subtraction procedure. The numerical integrations are performed with the use of the adaptive Monte Carlo integrator Vegas [20]. To validate the subtraction procedure, we use two sets of cut parameters:

Virtual NLO corrections
The virtual NLO corrections to the perturbative fragmentation function for g → QQ 8 , where the QQ pair is in a color-octet 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, like the one in figure 1, 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 operator vertex 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 equal momenta of the Q andQ by p and the momentum of the final-state gluon by q. The sum of the virtual one-loop cut diagrams at order α 3 s defines a scalar function A virtual (p, q, l, n) that is integrated over the loop momentum l. The virtual NLO contribution to the fragmentation function can be expressed as where N CS is the Collins-Soper prefactor in eq. (2.3) and dφ Born is the phase-space measure in eq. (2.11). The integral of A virtual over the loop momentum l is a homogeneous function of n of degree 2. It can be expressed as the product of [(2p + q).n] 2 /s 2 , where s = (2p + q) 2 , and a dimensionless function of s/m 2 and the momentum fraction z defined in eq. (2.12b). By means of standard tensor reduction techniques, the integral of A virtual over l can be reduced to [(2p + q).n] 2 /s 2 multiplied by a linear combination of one-loop scalar master integrals whose numerators are simply 1 and whose coefficients are functions of s, z, the JHEP01(2019)227 number of dimensions D = 4 − 2 , and the number of colors N c . The denominators of the master integrals come from at most three Feynman propagators with mass m or 0 and at most one eikonal propagator of the form 1/[(l + P ).n + i ], where P is a linear combination of p and q. The set of master integrals is the same as in the NLO calculation of the fragmentation function for g → QQ 1 in ref. [12]. Each master integral can be expanded as a Laurent expansion in to order 0 . The coefficients of the poles in can be evaluated analytically. Analytic expressions for the poles in for the master integrals with an eikonal propagator are given in appendix B of ref. [12]. The poles and the finite parts of the virtual NLO correction can be readily obtained from the calculation of the fragmentation function for g → QQ 1 in ref. [12] by changing the coefficients of the master integrals.
In all the poles in from the loop integral in eq. (4.1), the Born squared amplitude A Born (p, q) in eq. (2.13) 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. (4.3) with only UV poles: where C F = (N 2 c −1)/(2N c ) 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. (4.3) with mixed UV and IR poles: In Feynman gauge, this term comes from loop correction to the operator vertex. There are two terms in eq. (4.3) with only IR poles:

JHEP01(2019)227
The infrared poles originate from loop-momentum configurations with partons that can be soft and/or collinear. The term S 1 is a soft pole that in Feynman gauge comes from one-loop 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. (4.2) can be expressed as where D 1 (z) is the LO fragmentation function in 4 − 2 dimensions in eq. (2.19) and N dφA Born is the LO differential fragmentation function in eq. (2.18). Many of the IR poles cancel against terms in the real NLO corrections, which are given by the integrals of the subtraction terms in eq.

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

JHEP01(2019)227
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 schemes specified above read s /π in the beta function (µd/dµ)α s (µ). In the counterterm for g in eq. (5.3c), we have allowed for the renormalization scale µ R of α s to be different from the scale µ introduced through dimensional regularization.
The contributions to the NLO fragmentation function from inserting propagator counterterms and vertex counterterms into the LO cut diagrams are The terms with the coefficients C g , C Q (s), C gQQ , and C eikonal are associated with the virtualgluon propagator, the virtual-quark propagator, the quark-gluon vertices, and the eikonalgluon vertex in the LO cut diagrams, respectively. The expressions for these coefficients in terms of the rescaled renormalization constantsδ i 's are The counterterm contributions to the NLO fragmentation function can be reduced to The field renormalization constantsδ 2 andδ 3 have single IR poles that must cancel the single IR poles that remain after adding the real NLO corrections and the virtual NLO

JHEP01(2019)227
corrections. The linear combination 2δ 2 +δ 3 in eq. (5.6) has infrared poles proportional to C F , n f T F , and N c . The IR pole proportional to C F inδ 2 cancels the C F term in the IR pole in S 1 in eq. (4.6a). The IR pole proportional to n f T F inδ 3 cancels the IR pole in eq. (B.39). The IR pole proportional to N c inδ 3 cancels single IR poles in S 1 and S 2 (s, z) and single IR poles in eqs. (B.20b), (B.26), and (B.27b). This completes the verification of the cancellation of the IR poles.
The renormalization of the operator defining the fragmentation function introduces an additional counterterm. Its expression in the MS scheme reads 1 where P gg (y) is the Altarelli-Parisi splitting function: We have allowed for the factorization scale µ F to be different from the scale µ introduced through dimensional regularization. The UV poles in the sum of the real NLO corrections in eq. (3.19) and the virtual NLO corrections in eq. (4.7) must be canceled by the counterterms in eqs. (5.4) and (5.7). The UV poles in the NLO virtual corrections from the coefficients U g , U Q (s), and U gQQ in eqs. (4.4) are canceled by the corresponding counterterms C g , C Q (s), and C gQQ in eqs. (5.5). The UV pole from real NLO corrections in eq. (B.7b) is a convolution integral over y. It is canceled by the contribution to the operator counterterm in eq. (5.7) from the region y < 1. The contribution to the virtual NLO corrections from the coefficient U eikonal in eq. (4.4d) is canceled by the sum of the eikonal counterterm C eikonal in eq. (5.5d) and the contribution to the operator counterterm in eq. (5.7) from the endpoint y = 1. This completes the verification of the cancellation of the UV poles.
The complete NLO term D (NLO) g→QQ 8 (z) in the fragmentation function for g → QQ 8 is obtained by adding the real NLO corrections in eq. (3.19), the virtual NLO corrections in eq. (4.7), and the counterterms in eqs. (5.4) and (5.7), and then taking the limit → 0. After dividing by the perturbative NRQCD matrix element in eq. (2.6b), we obtain the NLO fragmentation function D

Numerical results
Beyond leading order in α s , the fragmentation function D g→η Q (z) depends on a factorization scale µ F and on a renormalization scale µ R . If we make those scales explicit, the expression for the color-octet 1 S 0 term in the NLO fragmentation function for g → η Q in eq. (2.4) is

JHEP01(2019)227
The scales µ F and µ R were introduced through renormalization. The Altarelli-Parisi evolution equation for the fragmentation function can be used to sum up large logarithms of µ F /m to all orders in α s . The solution to the evolution equation is sensitive to the behavior of the NLO fragmentation function near the upper endpoint z → 1. The fragmentation function includes terms with endpoint singularities of the form log(1 − z) and log 2 (1 − z). It is worthwhile to make those singularities explicit. Our final result for the NLO term in the fragmentation function has the form LO (z/y) where D LO (z) is the LO fragmentation function in eq. (2.9), D finite (z) is a function with a smooth limit as z → 1, and D (8) sing (z) is a function with logarithmic singularities in that limit: 3) The coefficients c finite (z) are calculated by subtracting the singular term D (8) sing (z) defined by eq. (6.3) from the fragmentation function D (8) NLO (z) with µ R = µ F = 2m: The resulting curves are shown in figure 4, both as a function of z with a linear scale (left panel) and as a function of 1−z with a logarithmic scale (right panel). By construction, the curves for D (8) finite (z) with n f = 0 and for the contribution to D finite (z) from one additional light-quark flavor approach plateaus when 1 − z gets very close to 1.
The NLO fragmentation function is compared with the LO fragmentation function in figure 5 for the case of bottomonium. The LO fragmentation function is proportional to α 2 s (µ R )D LO is given in eq. (2.9). The NLO fragmentation function is proportional to the sum of α 2 s (µ R )D LO (z) increases monotonically from 0 to 5α 2 s /(96m 3 b ) as z increases from 0 to 1. The NLO term α 3 s D NLO (z) increases from −∞ as z → 0 to a maximum at z ≈ 0.8, and then decreases to −∞ as z → 1. For µ R = µ F = 2m b , its maximum is 1.1 × 10 −3 /m 3 b at z = 0.80. At z = 0.5, the NLO fragmentation function is larger than the LO fragmentation function by a factor of 2.67. The NLO term is large and negative in the z → 1 region because of the log 2 (1 − z) term. To determine the fragmentation function accurately as a function of z in this region, it would be necessary to sum the leading logarithms of 1 − z to all orders. This is not essential because the logarithmic singularities as z → 1 are integrable. The NLO term is also large and negative in the z → 0 region. This may arise from a singular term of the form log z. The large NLO corrections in this region may not be a problem, because kinematics generally excludes contributions from the small-z region of the fragmentation function.
The sensitivity of the LO and NLO fragmentation functions to the renormalization scale µ R and to the factorization scale µ F is also illustrated in figure 5. The bands are obtained by varying µ R or µ F up or down by a factor of 2 around the central value 2m b (with the other scale held fixed). The left panel of figure 5 shows the NLO band from varying µ R (with µ F = 2m b ). The NLO band is significantly wider than the LO band in the central region of z. 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 right panel of figure 5 shows the band from varying µ F (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 , partly because the function multiplying log(µ 2 F /4m 2 ) in eq. (6.2) vanishes near z ≈ 0.75. The width of the band increases at larger z and near the z → 0 endpoint.

Summary
In this paper, we have presented the calculation of the NLO QCD corrections to the fragmentation function for a gluon into a color-octet 1 S 0 QQ pair at leading order in v. To calculate the real NLO correction, we developed a generalization of the FKS subtraction scheme that can be applied to fragmentation functions. The UV and IR poles in the real NLO corrections arise from boundaries in the phase space for the final-state partons. We constructed subtraction terms for the integrand of the real NLO correction that cancel the singular contributions from each of these boundaries and for which the poles in the subtraction integrals can be calculated analytically. Our development of an FKS subtraction scheme for fragmentation functions paves the way toward the automated calculation of the NLO fragmentation functions in other NRQCD channels.
We found that the NLO QCD corrections have a dramatic effect on the fragmentation function. Instead of increasing monotonically with z as at LO, the NLO fragmentation function has a maximum in the central region of z. For µ F = µ R = 2m, the NLO fragmentation function is larger than the LO fragmentation function by about a factor of 3. As a consequence, the NLO fragmentation function displays strong sensitivity to the renormalization scale. These results suggest that QCD corrections to the color-octet 1 S 0 fragmentation function could have a significant impact on the production of the spin-singlet quarkonium states η c and η b at large transverse momentum.

JHEP01(2019)227
A Two-parton phase-space measures The real NLO corrections to the fragmentation function involve integrals over the phase space for two massless partons whose longitudinal momenta are constrained to add up to (K − 2p).n. The dimensionally regularized phase-space measure is where the single-particle phase-space measure in D dimensions is Explicit parametrizations of this phase-space measure are required in order to extract the poles in in the integrals of the subtraction terms for the real NLO corrections.
In appendix A of ref. [12], two parametrizations of the two-parton phase-space measure were derived. The parametrization in eq. (A.10) of ref. [12]  The subtraction term D (u j ,δ j ) 1 in eq. (3.47) is the product of the Born squared amplitude A Born (s j , z), wheres j is defined in eq. (3.46), and a simple function of two dimensionless variables: To obtain an expression for the phase-space measure in eq. (A.1) that is differential in w j and u j , we begin with the expression for the phase-space measure in eq. (A.6) of ref. [12], which depends on two light-like 4-vectors k 1 and k 2 that define polar axes for the parton momenta q 1 and q 2 : where the scalar variables ρ j are defined by The first light-like 4-vector k 1 should be a linear combination of p µ and n µ with coefficients that are scalar functions of p and n. The second light-like 4-vector k 2 should be a linear

JHEP01(2019)227
combination of p µ , n µ , and the first parton momentum q µ j with coefficients that are scalar functions of p, n, and q j . The transverse angular measure dΩ 2⊥ is that for the second parton momentum.
The appropriate choices for the two light-like 4-vectors in this case are The first scalar variable ρ j can be expressed as By using the identity q 3−j = (q 1 +q 2 )−q j , the second scalar variable ρ 3−j can be expressed as After changing variables from ρ j and ρ 3−j tos j and w j , we obtain our final result for the two-parton phase-space measure: where dφ Born is the Born phase-space measure defined in eq. (2.16) and dφ is the parton-emission measure (A.10) The range of w j is from (1 − z)u j /z to ∞. The range of u j is from 0 to 1.

A.2 Phase-space measure for D (λ)
After averaging over transverse angles, the subtraction term D (λ) in eq. (3.45) can be expressed as the product of the Born squared amplitude A Born (s, z), wheres = (2p +q) 2 andq is the 4-vector defined in eq. (3.6), and the sum of two functions of dimensionless variables. In the second function, the two dimensionless variables are .
To obtain an expression for the phase-space measure in eq. (A.1) that is differential in λ and v, we begin by expressing the phase-space measure in an iterated form involving the product of the phase-space measure for Q = q 1 + q 2 and the differential phase space for the decay of a particle with 4-momentum Q into particles with 4-momenta q 1 and q 2 : where the differential 2-particle phase space is

JHEP01(2019)227
The product of the phase-space measure for Q and the delta function in eq. (A.12) is most easily simplified in the rest frame of P . The phase-space measure for Q is where θ is the angle between Q and n in that frame. The energy Q 0 and the magnitude |Q| of the 3-momentum can be expressed in terms of Lorentz invariants: where λ(a, b, c) = a 2 + b 2 + c 2 − 2(ab + bc + ca). After dividing each term by 2p.n, the constraint provided by the delta function in eq. (A.12) can be expressed as The delta function can be used to integrate over cos θ. It is convenient to change variables from s = (2p + Q) 2 tos = (2p +q) 2 : The resulting expression for the product of the phase-space measure and the delta function is The differential 2-body phase-space measure dΦ(Q, q 1 , q 2 ) is most easily simplified in the rest frame of Q. It can be reduced to where θ is the angle between P and q 2 in that frame and dΩ ⊥ is the transverse angular measure. The dimensionless variable v defined in eq. (A.11) can be expressed in terms of cos θ and Lorentz invariants: The differential phase space then reduces to Inserting eqs. (A.18) and (A.21) into eq. (A.12), we obtain our final result for the two-parton phase-space measure: dφ real (p, q 1 , q 2 ) = dφ Born (s, z) dφ

JHEP01(2019)227
where dφ Born is the Born phase-space measure defined in eq. (2.16) and dφ 2 is the partonemission measure The variables v + and v − depend on λ ands:

B Subtraction integrals
Explicit expressions for the poles in the integrals of the subtraction terms can be obtained by carrying out the integration over a (3 − 2 )-dimensional slice of phase-space. The choice of this slice depends on the subtraction term under consideration, which is labelled by a superscript (A) that indicates the variables that vanish in the singular region. The general decomposition of the phase-space measure reads The Collins-Soper prefactor N CS is defined in eq. (2.3). The Born phase-space measure dφ Born , which is differential in s (A) , is defined in eq. (2.16). The parton-emission measure dφ (A) can be reduced to a differential in two variables multiplied by a transverse angular measure dΩ ⊥ , whose integral is For each subtraction term (A), the choices of the arguments s (A) and z (A) of dφ Born and the choice of the measure dφ (A) are designed to ease the extraction of the poles in .
B.1 Subtraction term proportional to θ (ζ j ) In the subtraction term in T 1,j that is proportional to θ (ζ j ) , the function D (ζ j ) is given in eqs. (3.37), where the tensor V (ζ j ) µν is given in eq. (3.38). To extract the poles in the integral of the subtraction term, we use the decomposition of the phase-space measure in eqs. (3.7) and (3.8) of ref. [12]. The factor N (p, q j ) is N CS y −2+2 j , where y j is the momentum fraction defined in eq. (3.15). The factor dφ Born (p, q j ) can be expressed in terms of the measure defined in eq. (2.16) as (1/y j )dφ Born (s j , z/y j ), where s j = (2p + q j ) 2 . In the decomposition of the phase-space measure in eq. (B.1), the Born phase-space measure is dφ Born (s j , z/y j ) and the parton-emission measure is

JHEP01(2019)227
The extraction of the poles follows the same approach as in section 3.3 of ref. [12]. The only difference is that the phase-space boundary s > s j /y j is replaced by the cut s > (s j /y j )/ζ cut . The tensor V (ζ j ) µν can be averaged over transverse angles using results in section 3.3 of ref. [12]. The angular average of D (ζ j ) is (B.4) The integral over s gives a UV pole. Up to terms of relative order 3 , it can be expressed as (B.5) The integral over y j gives an IR pole. After renaming the differential variable in dφ Born as s j → s, the subtraction integral is the same for j = 1 and j = 2. Their sum is The functions I where P (real) gg (y) is the real-gluon contribution to the Altarelli-Parisi splitting function for g → g: The function D 1 is defined in eq. (2.19). The functions D log and D log 2 are defined by where the measure N dφA Born (s, z) is given in eq. (2.18). Note that D log (z, A cut ) depends on A cut through the term log(A cut )D 1 (z).

JHEP01(2019)227
dφ Born (s, z) but different parton-emission measures. The V (λ) 1 term is integrated using the parton-emission measure dφ term, the integrals over λ and u both give IR poles. For the V (λ) 2 term, the integrals over λ and v give a double IR pole. The complete subtraction integral is The functions I The function D 1 is defined in eq. (2.19), the functions D log and D log 2 are defined in eqs. (B.9), and the function F (λ) is Its limiting behavior as λ → 0 is v − → λ/r 2 .

B.6 Subtraction terms proportional to θ (u j )
In the subtraction term proportional to θ (u j ) in T below by considering appropriate choices for s (A) and dφ (A) for each subtraction term. It is therefore convenient to consider these two subtraction terms together.
To extract the poles in the integral of the D (u j ) 1 subtraction term given in eq. (3.36), we use the decomposition of the phase-space measure in eq. (B.1) with the Born phase-space measure dφ Born (s j /y j , z), where s j = (2p + q j ) 2 and y j is defined in eq. (3.15), and with the parton-emission measure where u j is defined in eq. (3.21) and λ j is defined in eq. (3.35). The parton-emission measure dφ in eq. (B.30) can be derived from dφ (ζ j ,u j ) in eq. (B.10) by changing variables from s to λ j . To extract the poles in the integral of the D (u j ) 2 subtraction term given in eq. (3.42), we use the decomposition of the phase-space measure in section B.3 with the Born phasespace measure dφ Born (s, z) and with the parton-emission measure dφ The integral of the subtraction term S over the Born phase space and over the variable λ is The integral of the subtraction term S over the Born phase space and over the variable λ j can be expressed in a similar form by renaming integration variables s j /y j →s and λ j → λ: 1,j in eq. (3.34), with the substitutions s j /y j →s and λ j → λ, and S (u j ) 2,j in eq. (3.41) are weight functions that add up to 1. The sum of the integrals in eqs. (B.31) and (B.32) therefore factors into a Born phase-space integral and an integral over λ. With dimensional regularization, the integral over λ is 0, but it can be expressed as the difference between a UV pole and an IR pole. The integral over u j gives an IR pole. The sum of the two subtraction integrals is the same for j = 1 and j = 2. The sum over j of the subtraction integrals is This constant vanishes if UV = IR .

JHEP01(2019)227
B.7 Over-subtraction term proportional to θ (λ) θ (u j ) In the over-subtraction term proportional to θ (λ) θ (u j ) in T (gg) 2,j , the function D (λ,u j ) is given in eq. (3.49b). To extract the poles in the integral of the over-subtraction term, we use the same decomposition of the phase-space measure as in section B.3, with the Born phase-space measure dφ Born (s, z) and the parton-emission measure dφ (λ,u j ) = dφ (u j ,δ) in eq. (B.13).
The integrals over λ and over u j give IR poles. The over-subtraction integral is the same for j = 1 and j = 2. Their sum is where the constant V (λ,u) is (B.36)

B.8 Subtraction term involving light quarks
In the subtraction term T (qq) in eq. (3.50), which is proportional to θ (λ) , the function D (qq) is given in eqs. (3.51). It is proportional to the contraction of the Born tensor A µν gluon in eq. (2.23) and the tensor P (qq) µν in eq. (3.12). To extract the poles in the integral of the subtraction term, we use the same decomposition of the phase-space measure as in section B.3, with the Born phase-space measure dφ Born (s, z) and the parton-emission measure dφ (qq) = dφ (u 1 ,δ) defined in eq. (B.13).
The tensor P (qq) µν in eq. (3.12) can be averaged over the transverse angles using results in section 3.5 of ref. [12]. The angular average of D (qq) is where the constant V (qq) is

C Helicity amplitudes with MadGraph
We use MadGraph5 [18] to generate the amplitude associated with the sum of cut diagrams with an S-wave color-octet heavy-quark pair and two light partons in the final state. The squared amplitudes generated by MadGraph5 are expressed in the helicity basis, and they are summed over the sets of helicities that give non-zero contributions. Such a computation in terms of helicity amplitudes is particularly well suited for numerical purposes.

JHEP01(2019)227
The MadGraph5 generator is primarily aimed at generating matrix elements and events for scattering or decay processes. The input is a UFO model [21] that encodes all the Feynman rules associated with the process. The Feynman rules associated with fragmentation processes can be translated into a UFO model, so the MadGraph5 can also be used to generate the matrix elements associated with fragmentation cut diagrams. A convenient precedure is to append the following additional (fictitious) particles and interactions to the UFO model associated with the Standard Model: Particles: the new states have the following properties: • the state eik is a zero-mass, scalar, color-octet state; • the state source is a zero-mass, vector, color-singlet state; • the state co is a duplicate of a heavy-quark pair state (either charm or bottom) that is projected onto a color-octet state.
Interactions: the new interactions specify the Lorentz structure of the coupling of the gluon field with the source and/or the eikonal lines. Two vertices must be defined: • the source-eik-g vertex, • the eik-g-eik vextex.
The additional Feynman rules for fragmentation functions include a delta function for the cut through the eikonal line. That delta function is included in the phase-space measure, so it does not appear in the UFO model. The relabelling of the heavy-quark+antiquark state co simply acts as a trigger to force the heavy quark and antiquark in the final state into a color-octet configuration, a projection that is applied by means of an ad-hoc modification of the color treatment in MadGraph5. Given that we are only interested in a final-state heavy-quark pair projected onto a scalar state in this work, we also add a filter to remove diagrams where a single gluon is attached to the co fermion line.
With the above UFO model at hand, MadGraph5 can be used to generate (i) the diagrams associated with a fragmentation process and (ii) a Fortran code for the numerical evaluation of the amplitude associated with the cut diagrams. The left-hand side of the cut of each generated diagram for gluon fragmentation with two gluons crossing the cut in addition to the 1 S 0 QQ 8 state is displayed in figure 6.
The projection of the heavy-quark pair onto a spin-singlet state has not yet been applied. Moreover, with the above extension of the UFO model for SM processes, Mad-Graph5 has treated the eik and source states as if they were regular particles. This calls for several modifications in the default Fortran code generated by MadGraph5: • The projection onto a spin-singlet state can be easily achieved by combining amplitudes with different heavy-quark helicities in the appropriate configurations.
• MadGraph5 automatically decomposes the calculation of helicity amplitudes into generic building blocks, and writes the routines associated with these blocks using JHEP01(2019)227

JHEP01(2019)227
a module called Aloha [22]. In the two routines associated with the two vertices involving the eikonal line, MadGraph5 identifies the four-vector n µ associated with the eikonal line with the momentum flowing along the eikonal line. Thus the two routines must be modified so that the four-vector n µ is properly set to the (fixed) eikonal four-vector.
• A cut diagram with Lorentz index µ at the operator on the left side of the cut and Lorentz index ν at the operator on the right side of the cut must be contracted with −g µν . Instead, MadGraph5 contracts each index with the polarization vector µ (λ) [or * ν (λ)] built upon the helicity state λ of the source. Summing over the helicity states of the source will not yield the correct result in general, unless we define the momentum K associated with the source with some care. For a fixed value of K.n, we define light-like four-vectors K + and K − whose spacial components are equal and opposite and such that K + .n = K.n and K µ − is proportional to the eikonal fourvector n µ . The tensor that is used by MadGraph5 to contract indices µ and ν at the sources on both sides of the cut can be decomposed as follows: However by gauge invariance, any amplitude associated with the sum of cut diagrams with Lorentz index µ and ν at the sources on both sides of the cut vanishes if it is contracted with either the eikonal four-vector n µ or the eikonal four-vector n ν . Since K ν − is porportional to n µ , the second term in eq. (C.1) will not contribute after contraction of the Lorentz indices.
After implementing the above modifications, the code generated by MadGraph5 can be used to evaluate the matrix element of processes involving gluon fragmentation into an S-wave color-octet heavy-quark pair plus other light partons. Gauge invariance and Lorentz invariance have been checked numerically. We also checked that the above procedure adapted to the case of fragmentation into a color-singlet heavy-quark pair reproduces numerically the squared amplitudes calculated analytically in ref. [12].
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.