Covariant calculation of a two-loop test of NRQCD factorization

We test the nonrelativistic QCD factorization conjecture for inclusive quarkonium production at two loops by carrying out a covariant calculation of the NRQCD long-distance matrix element (LDME) for a heavy-quark pair in an $S$-wave, color-octet state to fragment into a heavy-quark pair in a color-singlet state of arbitrary orbital angular momentum. The NRQCD factorization conjecture for the universality of the LDME requires that infrared (IR) divergences that it contains be independent of the direction of the Wilson line that appears in its definition. We find this to be the case in our calculation. The results of our calculation differ in some respects from those of a previous calculation that was carried out by Nayak, Qiu, and Sterman using light-cone methods. We have identified the sources of some of these differences. The results of both calculations are consistent with the NRQCD factorization conjecture. However, the general principle that underlies this confirmation of NRQCD factorization at two-loop order has yet to be revealed.


I. INTRODUCTION
The nonrelativistic quantum chromodynamics (NRQCD) factorization conjecture 1 for quarkonium production in hard-scattering QCD processes [3] postulates that the production rates for those processes at large momentum transfer can be written as a sum of products of short-distance coefficients (SDCs) and NRQCD long-distance matrix elements (LDMEs).
The SDCs are process dependent, but they contain no infrared (IR) divergences, and can be calculated in perturbation theory. The LDMEs contain all of the IR sensitivity of the process. They are generally nonperturbative in nature, but they are conjectured to be universal (process independent)-a property that gives NRQCD factorization much of its predictive power.
The NRQCD factorization conjecture has enjoyed considerable phenomenological success.
(See, for example, Refs. [4][5][6].) However, there are also processes for which theory and experiment are in considerable tension. (See, for example, Refs. [5,6].) Although the NRQCD factorization conjecture has been in existence for many years, there is still no demonstration that it is correct to all orders in perturbation theory or that it fails in perturbation theory. Important progress toward an all-orders proof of NRQCD factorization was given in Ref. [7]. There, an all-orders argument was given, for the leading and first subleading powers of m H /p T , that quarkonium production rates can be written as a sum of SDCs convolved with fragmentation functions for one or two partons to fragment into a quarkonium. Here, m H and p T are the mass and transverse momentum of the quarkonium.
A proof of NRQCD factorization would require the further factorization of the fragmentation functions into sums of products of SDCs with NRQCD LDMEs. One difficulty in carrying out this factorization is that the fragmentation functions at the scale of the heavyquark mass m Q contain processes in which gluons are radiated with energies of order m Q in the quarkonium center-of-momentum (CM) frame. Such gluons cannot be absorbed into the NRQCD LDMEs, which contain radiation only at the small scale m Q v, where v is the typical velocity of the heavy quark Q or antiquarkQ in the quarkonium CM frame. On the 1 In this paper, when we refer to "NRQCD factorization," we mean NRQCD factorization for inclusive quarkonium production. A proof of NRQCD factorization for exclusive quarkonium production at leading power in the inverse of the hard-scattering momentum transfer is given in Refs. [1,2]. It is generally conceded that NRQCD factorization for quarkonium decays can be established along the lines of the argument that is given in Ref. [3], although no complete proof has been published.
other hand, soft gluons with arbitrarily small momenta can connect the order-m Q gluons to the Q orQ, and these soft gluons must be absorbed into the NRQCD LDME if NRQCD factorization is to hold. A resolution of this dilemma was suggested in Ref. [7]. It is based on modifying the original definition of an LDME to include a Wilson line, which acts as a proxy for an order-m Q gluon in that its interactions with additional soft gluons are identical to those of the order-m Q gluon. This so-called "gauge completion" of the LDMEs is also required in order to make the LDMEs manifestly gauge invariant [7]. (In this paper, we will always refer to the gauge-completion Wilson lines as "Wilson lines" in order to distinguish them from soft approximations to the heavy-quark lines, which we refer to as "eikonal lines.") There are no apparent obstacles to establishing that all of the soft singularities in the fragmentation functions can be absorbed into the gauge-completed NRQCD LDMEs, given that the effective field theory NRQCD is a valid approximation to QCD in the soft limit and that the gauge-completed LDMEs account for the soft interactions with the order-m Q gluons. The central problem in proving NRQCD factorization is then to demonstrate that the soft divergences in the LDMEs are independent of the direction(s) of the gauge-completion Wilson line(s). Without such a demonstration, the LDMEs would depend on the directions of the order-m Q gluons, and so universality of the LDMEs would be lost, even for a single type of quarkonium production process.
Perturbative tests at one and two loops of the proposition that the soft divergences in the LDMEs are independent of the direction of the gauge-completion Wilson line have been provided by Nayak, Qiu, and Sterman [7][8][9]. In Refs. [7,8], one-and two-loop contributions to the LDMEs were given at the leading nontrivial order in v (order v 2 ), with the details of the calculations being given in Ref. [7]. In Ref. [9], the one-and two-loop contributions were calculated to all orders in v, although explicit expressions were given only for the non-Abelian diagrams, which are the most complicated to evaluate. These calculations confirmed that the soft divergences in the LDMEs are independent of the Wilson-line direction through two loops. In the remainder of this paper, we refer collectively to Refs. [7,9] as "NQS".
The calculations in Refs. [7][8][9]  The calculations are quite complicated, especially for the non-Abelian diagrams, and the independence of the soft divergences from the Wilson-line direction emerges only after many complicated intermediate expressions are summed. Therefore, it seems worthwhile to check these calculations using manifestly covariant methods. One could hope that such covariant calculations might give insights into the independence of the soft divergences from the Wilson-line direction.
In this paper, we present such a manifestly covariant calculation of the soft divergences in an LDME at two-loop order. The calculation relies on a new UV regulator for the phase-space integrations that avoids the use of a hard cutoff. As we will see, our covariant calculation is actually considerably more complicated than the calculations in NQS. This is probably because the light-cone methods in NQS allow one to cancel contributions of realreal diagrams against contributions of real-virtual diagrams and to cancel certain real-virtual contributions by symmetry, without having to calculate them completely. In contrast, in our calculations, we compute all of the IR contributions of the real-real diagrams and real-virtual diagrams explicitly and implement cancellations only at the end. Consequently, we must deal with soft poles in dimensional regularization in D = 4 − 2 dimensions of order 1/ 4 , 1/ 3 , 1/ 2 , and 1/ in order to arrive at a final result of order 1/ . Although we can eliminate the 1/ 4 contribution trivially at the start by making a Ward-identity re-arrangement of the numerator structure, we still must work out the cancellation of the 1/ 3 and 1/ 2 poles, while the NQS calculation contains only 1/ 2 imaginary poles, which cancel in the absolute square of the amplitude, and 1/ real poles.
In NQS, it was found that IR poles in the LDME are independent of the Wilson-line direction, and we also find that to be the case. Our result for the non-Abelian diagrams agrees with that in NQS, aside from an overall sign. Our result for the Abelian diagrams differs from that in NQS in two respects. First, we find new contributions that have the form of a one-loop contribution to an SDC times the one-loop IR-divergent contribution to the LDME. For these contributions, the IR pole in the LDME is independent of the direction of the Wilson line. We conclude that some of these contributions are absent in NQS because a mismatch between virtual-gluon and real-gluon contributions in the UV region was neglected. We also find a new contribution that has the form of a two-loop IR-divergent contribution to the LDME. Again, the IR pole in the LDME is independent of the direction of the Wilson line. This new contribution can be reconciled with the result in NQS if we re-interpret a UV contribution in NQS as an IR contribution.
The remainder of this paper is organized as follows. In Sec. II we describe the LDME and its Feynman rules, the classes of diagrams that we calculate, and the kinematics. Section III contains a description of the covariant phase-space regulator that we use and also contains convenient formulas for the phase-space integration in dimensional regularization. The main body of our work is in Sec. IV, in which we give a description of the calculations for the various diagrams at a level of detail that would allow the reader, with some considerable effort, to reproduce our results. Our results are summarized in Sec. V, and we compare them with the results in NQS in Sec.VI. Finally, we give our conclusions in Sec. VII.

II. PRELIMINARIES
The gauge-completed LDME is given by [9] O H n (0) = 0|χ † (0)K n,e ψ(0)Φ where ψ † and χ are two-component Pauli fields that create a heavy quark and a heavy antiquark, respectively, K n,e and K n,a are local combinations of spin and color matrices and polynomials in covariant derivatives, and a † H is the operator that creates the quarkonium H. The operators Φ (A) (0) are Wilson lines, which are defined by where P stands for path ordering, g s = √ 4πα s is the QCD coupling constant, A (A)µ = A µ a t a is the gauge field in the adjoint representation, (t a ) bc = −if abc , and µ is the momentum of the Wilson line.
Following NQS, we consider the LDME for a QQ pair in an S-wave, color-octet state to evolve into a QQ pair in a color-singlet state of arbitrary orbital angular momentum. That is, we consider the matrix elements [9] M (8→I) (P 1 , P 2 , ) = where the superscript I = 1, 8 denotes the color of the QQ pair, the T (q) i are the generators of color SU(3) in the fundamental representation, and P 1 and P 2 are momenta of Q and Q , respectively. We do not explicitly consider the spin state of the QQ pair, as the soft approximation for the Q andQ lines is independent of the spin.
Since we wish to test for a dependence on the direction of the LDME Wilson line, we need to consider only those diagrams in which at least one gluon attaches to the Wilson line.
Furthermore, it is clear that a diagram has a nonzero color factor only if one or more gluons crosses the final-state cut. The specific classes of diagrams that we compute are illustrated in Figs. 1-6. We use the notation X P j P k i for the contributions of these diagrams. Here, X i denotes the class of the diagram. For each class, the subscript i labels the position of the final-state cut, and the superscripts P j P k specify the gluon connections to the Q andQ lines, as we will describe below. The classes are defined as follows.
• A P j P k i ( Fig. 1) designates normal ladder diagrams. P j indicates the left-most gluon connection to the heavy-quark lines; P k indicates the right-most gluon connection to the heavy-quark lines.
• B P j P k i ( Fig. 2) designates crossed ladder diagrams. P j indicates the left-most gluon connection to the heavy-quark lines; P k indicates the right-most gluon connection to the heavy-quark lines.
• C P j P k i • F P j P k i ( Fig. 6) designates the non-Abelian diagrams. P j indicates the left-most connection of a gluon to the heavy-quark lines; P k indicates the right-most connection of a gluon to the heavy-quark lines.
The classes of diagrams are constructed such that one can obtain all of the contributions that we calculate by summing over the indices i, j, and k and including hermitian-conjugate diagrams, except for the classes A and B, for which the sums over the indices i, j, and k already include the hermitian-conjugate contributions. We note that the assignments of our loop momenta k 1 and k 2 for the diagrams C i , D i , E i , and F i are different from those of NQS.
We define the total momentum of the QQ pair to be 2p and the relative momentum to QQ rest frame to be q, and so In the rest frame of the QQ pair, p and q are given by where E Q = m 2 Q + q 2 . The relative velocity of the Q andQ is In comparing with the light-cone calculation in NQS, we need to make use of light-cone momentum coordinates for a D-dimensional vector V = (V 0 , V 1 , · · · , V D−1 ), which we take to be Then, a scalar product of two vectors V and W is given by We note that, in NQS, the momentum of the Wilson line is specified to be along the minus light-cone direction: However, our covariant calculation does not make use of this assignment.
We define the following Lorentz-invariant quantities, which appear throughout our calculations: Since we are interested in analyzing the soft singularities, we take soft approximations for the interactions of the gluons with the Q andQ lines. As we have mentioned, we refer to the Q andQ lines in the soft approximation as "eikonal lines," and we refer to the gaugecompletion Wilson lines as "Wilson lines," in order to distinguish them from each other.
The Feynman rules for the eikonal lines and Wilson lines are given in Refs. [7,9]. They are summarized below.
• The interaction between a gluon and a Q orQ line to the left of the final-state cut is shown in Fig. 7. The propagator for a Q orQ line to the left of the cut is given by and the vertex is given by where α is the Lorentz index of the vertex, a is the color index of the vertex, b and c are the color indices of the Q, µ is introduced to account for the dimensionality of the coupling constant in D = 4 − 2 dimensions, and the diagrammatic momentum of thē Q line is opposite to the physical momentum. The rules for propagators and vertices to the right of the final-state cut are obtained by taking the hermitian conjugates of the rules that are shown. Note that the product of a propagator and a vertex is invariant with respect to a change of the scale of P i .
• The interaction of a gluon with a Wilson line is illustrated in Fig. 8. On the left side of the final-state cut, the Wilson-line propagator is given by and the vertex is given by The Feynman rules for the non-eikonal and non-Wilson-line parts of the diagram are the usual ones. We follow the conventions in [10], which are consistent with the conventions that we have chosen for the eikonal and Wilson lines. We work in the Feynman gauge throughout.
We find that our expressions for the diagrams C i , D i , E i , and F i , which contain one gluon connection to the Wilson line, differ by an overall sign from the expressions for the corresponding diagrams in NQS.

III. PHASE SPACE
In this section we discuss the covariant phase-space regulators that we employ and also present convenient formulas for carrying out the phase-space integrations in dimensional regularization.
A In the fragmentation function in which the LDME is embedded, there are restrictions on the final-state phase space that follow from Dirac δ-functions that express conservation of the light-cone energy: , when the gluon with momentum k 2 is real, , when the gluons with momenta k 1 and k 2 are real.
Here k is the momentum of the parton (gluon) that fragments into the QQ pair. In the fragmentation function at a scale of order m Q , the light-cone energy that is available to the final-state gluons, · k = · (k − P 1 − P 2 ), is also of order m Q .
Motivated by this constraint on the fragmentation-function phase space, we impose a UV cutoff on the phase space of the LDME. In NQS, there is also a UV cutoff on the phase space of the LDME. It is imposed as a hard cutoff on · k (and on the remaining components of k ). We find it calculationally more convenient to provide a cutoff of order m Q on · k by applying a weight function to the available light-cone energy · k : where Λ is a cutoff parameter of order m Q . Then, integrating over the constraints in Eq. (15), when the gluon with momentum k 2 is real, and when the gluons with momenta k 1 and k 2 are real. We call the factors on the right side of Eq. (17) our standard phase-space regulators. For some parts of the calculation, as we will explain later, we will find it necessary to temporarily impose additional UV regulators in order to ascertain IR/UV nature of poles in .

B. Phase-space integration
In computing the phase-space integration for the final-state gluons, it is convenient to extend the range of integration to infinity, relying on UV regulators to remove UV poles.
Then, we obtain the following phase-space integration formulas: where p 0 is the temporal component of p, p is the vector of spatial components of p, and we define the measure of the phase-space integration as In all of the calculations in this paper, we have arranged the parameter integrations so that

IV. ANALYSES OF THE DIAGRAMS
In this section, we outline our calculations for the various classes of diagrams.

A. Method of calculation
Our general method of calculation is as follows.
For the Abelian diagrams, we carry out the k 1 integration first, holding k 2 fixed. We control UV divergences by imposing our standard phase-space regulators, and, in some cases, temporarily impose additional UV regulators in order to determine the IR/UV nature of the poles in . We first combine the denominators by using Feynman parameters that run from 0 to 1 to combine denominators that have a common term involving k 1 and by using Feynman parameters that run from 0 to ∞ to combine the remaining denominators. We then carry out the integration over k 1 , using standard formulas of dimensional regularization if k 1 is virtual and using Eq. (18) if k 1 is real. We identify the singular regions of the parameter integrals and isolate these as integrations over a single parameter, either by rescaling the parameters that run from 0 to ∞ or, in a few cases, by using sector decomposition. We then expand the integrals around the singularities by making use of formulas such as which applies when the domain of integration is 0 ≤ x ≤ 1. Here, the plus distribution for a regular function f (x) in x ∈ [0, 1]. However, we avoid expanding factors involving k 2 in powers of , as the complete dependence of these factors is needed to find the coefficient of the soft pole from the k 2 integration.
We find for the Abelian diagrams that the k 1 integrations for the sum over a class plus hermitian conjugate are IR finite when k 2 is fixed. 2 In some cases, there are UV poles from the k 1 integration. The IR finiteness of the k 1 integration when k 2 is fixed is expected on general principles because the sum over cuts in a class is sufficient to effect a unitarity cancellation of singularities that occur when k 1 is soft relative to k 2 or k 1 is collinear to .
Some of the contributions from the k 1 integration remain finite as k 2 goes to zero, and these can be considered to be SDCs that multiply possible one-loop IR poles from the k 2 integration. In the cases of diagrams C and E, there are contributions that become singular as k 2 goes to zero, and these yield two-loop IR-divergent contributions to the LDME.
Once we have evaluated the k 1 integration, we can carry out the k 2 integration straightforwardly by combining denominators with Feynman parameters and using Eq. (18) to carry out the phase-space integration. The remaining parameter integrals are easily reduced, by using expansions of the type in Eq. (20), to elementary integrals.
In the case of the non-Abelian diagrams, the k 1 integration for the sum over the class plus hermitian conjugate is no longer IR finite when k 2 is held fixed. In this case, as we will explain, the unitarity cancellation requires that one carry out the k 2 integration, as well as the k 1 integration. Therefore, for the non-Abelian diagrams, we combine denominators involving both k 1 and k 2 using Feynman parameters that range from 0 to 1 or from 0 to ∞, as outlined above. We then carry out the k 1 and k 2 integrations, isolate the IR singularities in single parameter integrals by using re-scaling or sector decomposition, and evaluate the IR poles by using Eq. (20).
In most cases, we have checked the results from direct integration of the parameter integrals by using the Mellin-Barnes representation to decompose a denominator so as to obtain parameter integrations that can be carried out in terms of beta functions. We then evaluate the Mellin-Barnes integrations by isolating poles in using the methods of Tausk [11] or Smirnov [12] and computing the remaining finite integrals as a sum over residues in the complex plane of Mellin-Barnes integration variable.

B. A i diagrams
We do not need to evaluate the momentum integrations for A i diagrams, as the color factors vanish: The color factors of the B i are given by The Feynman amplitudes for the diagrams B P 1 P 2 i with our standard phase-space regulators are given by , where c and d are defined in Eq. (10) and we have suppressed unnecessary iε's in the denominators.
1. IR finiteness of the k 1 or k 2 integration We wish to establish that, for a certain combination of the B i , the integration over k 1 with k 2 fixed is IR finite, while for the remaining part of the B i the integration over k 2 with k 1 fixed is IR finite. For this purpose, we temporarily impose additional regulators to control all of the UV divergences.
It can be seen by power counting that our standard UV regulators render the integrations over · k 1 and · k 2 UV finite. However, there are still potential UV divergences that could arise from the integrations over the other components of k 1 and k 2 . In order to control these divergences, we introduce the following additional UV-regulator factor where the cutoff Λ is of order m Q . We denote the B i into which we have inserted this UV-regulator factor by B i Reg .
The B i contain rapidity (collinear) divergences that are associated with the vanishing of the denominators of Wilson-line propagators. 3 We expect these divergences to cancel, because of unitarity, in the sum over cuts of the B i . 4 We can separate the parts of B 1 that contribute to the cancellations of the rapidity divergences in B 2 and B 3 by making use of the following partial-fraction identity in B P 1 P 2 1 Reg : Then we can write 3 With appropriate deformations of the k 1 integration contours into the complex plane, these can be considered to be collinear-to-divergences. See, for example, p. 291 of Ref. [13]. 4 In general, one must sum over all cuts of a diagram in order to effect a unitarity cancellation. However, in the case of divergences that arise from the vanishing of specific propagator denominators, it is necessary only to sum over all cuts of the singular denominators. Other propagators are relatively far off shell and can be contracted to a point for the purpose of evaluating divergences. where , In Appendix C 1, we show that the k 1 integration with k 2 fixed in B P 1 P 2 Reg is IR finite. We also show that the integrand of B P 1 P 2 1a Reg + B P 1 P 2 2 Reg is the hermitian conjugate of the integrand of B P 2 P 1 1b Reg + B P 2 P 1 3 Reg with k 1 ↔ k 2 . It then follows that the k 2 integration with k 1 fixed in B P 1 P 2 1b Reg + B P 1 P 2

3
Reg is also IR finite. We note, for use below, that it also follows that B P 2 P 1 1b

Computation of the diagrams
Having established that k 1 integration with k 2 fixed in B P 1 P 2 1a Reg + B P 1 P 2 2 Reg is IR finite and that the k 2 integration with k 1 fixed in B P 1 P 2 1b Reg + B P 1 P 2

3
Reg is IR finite, we can remove the extra UV-cutoff factor (25) without introducing any ambiguities in dimensional regularization.
Carrying out the k 1 integrations, we obtain where we defineμ and γ E is the Euler-Mascheroni constant. Then, combining B P 1 P 2 1a and B P 1 P 2 2 in Eq. (29) and carrying out the y integration, we find that Carrying out the k 2 integrations and expanding the resulting expression through O( −1 ), we have where ζ(z) ≡ ∞ n=1 n −z is the Riemann zeta function. It is easily seen, from the Feynman rules, that Then, we have Since, as we have noted at the end of Sec. IV C 1, That is, there are no poles in the sum over the B i diagrams.

D. C i diagrams
The color factors of the C i diagrams are given by The expressions for the diagrams C i , with our standard phase-space regulators in Eq. (17), are given by where a and c are defined in Eq. (10).
1. IR finiteness of the k 1 integration in C 1 and C 2 It can be seen by power counting that our standard UV regulators render the integration over · k 1 UV finite. However, there are still potential UV divergences that could arise from the integrations over the other components of k 1 . In order to control these divergences, we introduce the following additional UV-regulator factor We denote the C i into which we have inserted this UV-regulator factor by C i Reg .
The k 1 integrations of the individual C i diagrams contain rapidity (collinear) divergences that arise when the denominators of the Wilson-line propagators vanish. We expect these divergences to cancel, by unitarity, in the sum over cuts in the C i .
In Appendix C 2, we have carried out the k 1 integrations of C P 1 P 2 1 Reg and C P 1 P 2 2 Reg . The results are where the first term comes from C P 1 P 2 1 Reg and the second term comes from C P 1 P 2 2 Reg . We see the explicit cancellation of the poles in that arise from the rapidity (collinear) divergences.

Computation of C 1 and C 2
Having established that the k 1 integration with k 2 fixed in C P 1 P 2 1 Reg + C P 1 P 2 2 Reg plus hermitian conjugate is IR finite, we can remove the extra UV-regulator factor in Eq. (38) without introducing any ambiguities in dimensional regularization. The expressions for the diagrams are given by where the k 1 integrations of C P 1 P 2 1 and C P 1 P 2 2 are defined by Introducing Feynman parameters and performing the k 1 integrations, we obtain where we define the k 2 -dependent denominators K 1 and K 2 as Then, carrying out the λ 1 integration in Eq. (42) by using the formula in Eq. (A8) and inserting our results for C P 1 P 2 1 and C P 1 P 2 2 into Eq. (40), we find that In Eq. (44), the factor (2P 1 · k 2 ) 2 in the denominator of the first integrand arises from the k 1 integration. It has an IR sensitivity, in that it becomes singular when k 2 goes to zero. It affects the strength of the IR pole that will appear in the k 2 integration. Consequently, all of the contributions from the first set of brackets in Eq. (44) should be regarded as IR in nature, except for the pure UV pole, which is obtained by setting the factor (2P 1 · k 2 ) 2 to unity. However, this pure UV pole has an imaginary coefficient, and cancels when we add the hermitian-conjugate contribution. The factor (2 · k 2 + Λ 2 ) 1+2 in the denominator of the second integrand in Eq. (44) also arises from the k 1 integration. However, in this case there is no IR sensitivity because this factor is finite as k 2 goes to zero. Consequently, all of the contributions from the second set of brackets in Eq. (44) are UV in nature.
Combining denominators in Eq. (44) by using Feynman parameters and carrying out the k 2 phase-space integration by using Eq. (18), we find that and where we have made change of variables λ 3 → ξλ 2 , carried out λ 2 integrations in terms of the beta function, and introduced the definitions It can be seen from Eq. (10) that a ≥ 1. Then, the ξ integrations in Eq. (45) can be carried out by making use of the formulas in Eq. (A9). The result is where we have omitted imaginary contributions, which cancel when we add the hermitian conjugate. The subscripts "IR" and "UV" are reminders of the origins of the contributions.
The factor labeled "IR" in the first term is a two-loop contribution to the LDME. The factor labeled "IR" in the second term is the one-loop contribution to the LDME (absent its color factor), which is computed in Appendix D. The category UV includes all IR-finite contributions, as well as UV-divergent contributions.
We can find all the other C P j P k i contributions from the relations Taking into account the color factors in Eq. (36), we obtain where we have kept only the real part because the imaginary contributions cancel when we add the hermitian-conjugate contribution.
The color factors of D i diagrams are given by The expressions for the diagrams D i , with our standard phase-space regulators in Eq. (17), are given by , 1. IR finiteness of the k 1 integration in D 1 and D 2 As in the case of the C i diagrams, it can be seen by power counting that our standard UV regulators render the integration over k 1 UV finite. However, there are still potential UV divergences that could arise from the integrations over the other components of k 1 . In order to control these divergences, we introduce the following additional UV-regulator factor We denote the D i into which we have inserted this UV regulator factor by D i Reg .
The k 1 integrations of the individual D i diagrams contain rapidity (collinear) divergences that arise when the denominators of the Wilson-line propagators vanish and also contain soft divergences. We expect these divergences to cancel, by unitarity, in the sum over cuts in the D i .
In Appendix C 3, we have carried out the k 1 integrations in D 1 Reg and D 2 Reg . The results are Combining the results for the k 1 integrations of D P 1 P 2 1 and D P 1 P 2 2 , we find that We see that the real double and single poles cancel in the sum over cuts, leaving only an imaginary pole that cancels when we add the hermitian conjugate. Hence, we find that the plus hermitian conjugate is IR finite.

Calculation of D 1 and D 2
Having established that the k 1 integration with k 2 fixed in D P 1 P 2 plus hermitian conjugate is IR finite, we can remove the extra UV-cutoff factor (25) without introducing any ambiguities in dimensional regularization. That is, we can assume that any poles that we encounter in the final result for the sum over diagrams are UV in origin. We note that, as can be seen from Eq. (51), the k 1 integration in D P 1 P 2 2 , without the additional UV regulator in Eq. (52), is scaleless and vanishes in dimensional regularization. That is, the IR and UV poles cancel. Therefore, where the k 1 integration D 1 is given by The λ 1 integration can be carried out exactly to obtain The k 2 -dependent denominator factor that arises from the k 1 integration, (2 · k 2 + Λ 2 ) 1+2 , is non-singular as k 2 goes to zero. Therefore, we conclude that the entire contribution from the k 1 integration is UV in nature.
Making use of the formula for the k 2 integration in Eq. (45), we obtain Here, as we have mentioned, the subscripts "IR" and "UV" on the brackets indicate the origins of the contributions, and the factor labeled "IR" is the one-loop contribution to the LDME, absent its color factor. At this point we should, in principle, discard any imaginary contributions, since our proof of the IR finiteness of the k 1 integration with k 2 fixed was valid only for D P 1 P 2 1 + D P 1 P 2 2 plus hermitian conjugate. However, the expression in Eq. (58) contains no imaginary parts. (The imaginary IR pole in D P 1 P 2 2 Reg in Eq. (53) is also present in D P 1 P 2 2 , but it cancels against an imaginary UV pole.) The relations in Eq. (48) also hold for the D P j P k i . Taking these relations into account, along with the color factors in Eq. (50), we obtain The color factors of E P j P k i are given by The expressions for the diagrams E i with our standard phase-space regulators in Eq. (17) are given by where It can be seen by power counting that the k 1 integrations of the individual diagrams are rendered UV finite by our standard phase-space regulators. There is no need in this case to introduce additional UV regulators.
The k 1 integrations of the individual E i diagrams contain rapidity (collinear) divergences that arise when the denominators of the Wilson-line propagators vanish. We expect these divergences to cancel, by unitarity, in the sum over cuts in the E i .
Applying Feynman parametrization to Eq. (62), we obtain Then, carrying out the k 1 integration and the parameter integrations, except for λ 1 integration in E 1 , we obtain where K 1 and K 2 are given in Eq. (43). The first λ 1 integration yields The second λ 1 integration is given in Eqs. (A8). Then, we obtain The denominator factors K 1+2 1 = (2P 1 ·k 2 ) 1+2 become singular as k 2 goes to zero, and so all of the terms in this expression yield IR contributions. Inserting these results into Eq. (61), we find that .
Making use of the formula for the k 2 integration in Eq. (45), we obtain where we have dropped imaginary contributions that cancel when we add the hermitianconjugate contributions.
The relations in Eq. (48) also hold for the E P j P k i . Taking into account these relations and the color factors in Eq. (50) we obtain Now we turn to the calculation of the non-Abelian diagrams F The color factor of the F P j P k i is given by We note that our expression for the color factor of the non-Abelian diagrams differs by a factor of (1/ √ N c ) 2 from the color factor that is given in Eq. (16) of Ref. [9] because the factor 1/ √ N c from the color-singlet projector was omitted in Eq. (16) of Ref. [9].
The expressions for the F P 1 P 2 i diagrams, with our standard UV phase-space regulator in Eq. (17) are given by The numerator factor N F is given by where V µ 1 µ 2 µ 3 is the triple-gluon vertex with all momenta q i 's flowing into the vertex.
The numerator factor can be written as The first two terms on the right side of Eq. (74) cancel the eikonal-propagator denominators P 2 · k 2 and P 1 · (k 1 + k 2 ), respectively, resulting in expressions that are independent of P 1 or P 2 . These expressions cancel when we sum over gluon connections to the quark and antiquark lines. (This is a manifestation of the graphical Ward identities.) Therefore, we drop these terms in the numerator and use a modified numerator factor It can be seen from power counting that the individual diagrams F P 1 P 2 i are UV finite with respect to the k 1 and k 2 integrations. UV divergences that might occur when the component of k 1 or k 2 that is parallel to becomes large cancel because the terms in the modified The k 1 and k 2 integrations contain IR divergences that arise from several sources: (1) the rapidity (collinear) divergence that occurs when the denominator · k 1 vanishes, (2) the divergence that occurs when k 1 is collinear to and k 2 is collinear to k 1 , (3) the divergence that occurs when k 1 is collinear to k 2 , (4) the divergence that occurs when k 1 is soft relative to k 2 , and (5)  We expect all of the remaining divergences to cancel by unitarity in the sum over diagrams, except for divergence (5), which is the object of our calculation. In comparison with the unitarity cancellations that we found in the Abelian case, the unitarity cancellations in the case of the F i diagrams are rather involved. In particular, when k 1 is soft or collinear to k 2 , the cancellations involve the hermitian conjugate of the original diagram with k 1 ↔ k 2 and P 1 ↔ P 2 . Rather than attempting to identify and implement all of the individual unitarity cancellations, we simply carry out a straightforward evaluation of the diagrams and cancel poles in the sum over diagrams. This approach leads to rather complicated intermediate expressions that contain poles of orders 1/ 3 and 1/ 2 that cancel to leave a final result of order 1/ .
Using Eq. (71a) with the modified numerator factor in Eq. (75), we obtain where Note that any terms in F α 11 and F α 12 that are proportional to α vanish upon contraction with the associated factors in Eq. (76). It is this property that eliminates the contributions of order 1/ 4 . Applying Feynman parametrization to Eq. (77), we obtain We first carry out the k 1 integration by making use of the phase space integration formula in Eq. (18). In the result, we introduce a Feynman parameter x to combine the two k 2dependent denominators and then carry out the k 2 integration by using Eq. (18) once again.
The results are Here, we have dropped the terms that are proportional to α because they cancel on contraction with the other factors in Eq. (76), and we have introduced the definitions We make a change of variables, replacing λ 2 , λ 3 , and λ 4 with Then, we carry out the λ 1 integration by making use of the integral formula in Eq. (A3) and Next, we carry out the parameter integrations, except for the ξ integration, by making use of Eq. (20). We find that the 1/ 3 , 1/ 2 , and 1/ contributions to F P 1 P 2 1 are given by the following expressions: where Li 2 (z) is the dilogarithm function Using Eq. (71b) with the modified numerator factor in Eq. (75), we have where Applying Feynman parametrization to Eq. (86), we obtain where we have made the translation We carry out the k 1 and k 2 integrations in Eq. (87) by using the standard formulas for the virtual-gluon loop integration and the phase-space integration formula in Eq. (18). The result is where and Now we make a change of variables, replacing λ 2 , λ 3 , λ 4 , and x with Then, we can express the λ 1 integrations in terms of the beta function by making use of the integral formula in Eq. (A3) to obtain are given by the following expressions: 4iπ + 12 log 2c − 7 log A + 10 log B AB , (93b) 3. Result for the IR poles in F 1 + F 2 Now we compute the IR poles of various orders in in F 1 + F 2 .
First, let us consider the 1/ 3 poles. From Eqs. (83a) and (93a), we find that where we have made a change of variable ξ → 1/ξ for F P 1 P 2 1 1/ 3 and used the definitions of A and B that are given in Eq. (46). Since the integrand in Eq. (94) is antisymmetric under P 1 ↔ P 2 , or c ↔ d, we find that the triple poles cancel after symmetrization under P 1 ↔ P 2 : The triple poles in F P 1 P 1 1 + F P 1 P 1 2 and F P 2 P 2 1 + F P 2 P 2 2 cancel in a similar fashion.
Next, let us consider the 1/ 2 poles. From Eqs. (83b) and (93b), we find that where We compute the combination F P 1 P 2 1 + F P 1 P 2 2 1/ 2 AB by making use of the same method that we used to compute the 1/ 3 contributions. The result is Also, one can show that The ξ integrations can be carried out through a straightforward, but tedious process, by partial-fractioning the denominators and factoring the arguments of the logarithms. We obtain a lengthy result, which we do not reproduce here, that contains logarithms and dilogarithms. This result can be greatly simplified by making use of the polylogarithm identities in Eq. (A10) to obtain This expression contains no real IR poles.
Next, let us consider the 1/ contribution. From Eqs. (83c) and (93c), we find that where Taking only the real parts, we obtain For the contribution that contains a dilogarithm in the integrand, we eliminate the dilogarithm by integrating by parts. Then, the ξ integrations can again be carried out by partialfractioning the denominators and factoring the arguments of the logarithms. This leads to an expression, which we do not reproduce here, that is hundreds of terms long and contains logarithms, dilogarithms, and trilogarithms. This result can be reduced, by making use of the polylogarithm identities in Eq. (A10) and symmetrizing under P 1 ↔ P 2 to obtain a remarkably simple form: Since this result shows no dependence on the UV cutoff, we conclude that it is entirely IR in nature. We have confirmed this conclusion by carrying out a calculation of the k 1 integration for F using light-cone methods [14].
We can find all the other F P j P k i contributions from the relations Taking into account the color factors in Eq. (70), we obtain

V. SUMMARY OF RESULTS
Now let us summarize the results of our calculations.

From Eq. (22) we have
and from Eqs. (23) and (35) we have That is, there are no IR divergences in the ladder diagrams.
From Eqs. (49) (59), (69), and (106), we have where we have used g 2 s = 4πα s . We remind the reader that the subscripts "IR" and "UV" on the brackets indicate the origins of the contributions and that the factor labeled "IR" is the one-loop contribution to the LDME, which is computed in Appendix D. The expression for C contains two contributions: (1) a UV factor times an IR factor, which is a one-loop contribution to an SDC times the one-loop IR-divergent contribution to the LDME; (2) a two-loop IR-divergent contribution to the LDME. The expression for D contains a UV factor times an IR factor, which is a one-loop contribution to an SDC times the one-loop IR-divergent contribution to the LDME. The expressions for E and for the non-Abelian contribution F are both two-loop IR-divergent contributions to the LDME. Note that all of the IR-divergent contributions to the LDME are independent of the direction of the Wilson line, and, so are consistent with the NRQCD factorization conjecture. All of the two-loop IR-divergent contributions to the LDME are of the same form, and their sum is nonzero.

A. The NQS calculation and its correspondence to our calculation
In NQS, the integration over the minus component of the loop momentum of a virtual gluon is carried out by closing the integration contour in the complex plane and using Cauchy's theorem. Some of the pole residues cancel against contributions in which the virtual gluon is replaced with a real gluon, and those contributions are not calculated. As we will see below, this cancellation is not exact because of sensitivity of the real-gluon contribution to the UV regulator.
In the calculations in Ref. [7], a velocity expansion is made for the interactions of the gluons with the Q andQ lines. This expansion mixes the contributions of some of the diagrams that appear in our calculations. Consequently, there is not a one-to-one correspondence between the diagrams in NQS and our diagrams. However, the following correspondences hold.
Here, the roman numerals on the left sides of the correspondences are the notations for the diagrams in Ref. [7].

B. Comparison
Our result for the sum of the ladder and crossed ladder diagrams A + B contains no IR poles, which is in agreement with the order-v 2 results for I + II in Ref. [7].
Our result for the non-Abelian diagrams F is in agreement with the corresponding result in Ref. [9], except for an overall sign. The result for the non-Abelian diagrams in Ref. [9], when expanded to order-v 2 , is identical to the result for III in Ref. [7].
Our result for the Abelian diagrams that involve one interaction on the Wilson line is contained in C + D + E. In order to compare this result with the calculation of IV + V + V I in Ref. [7], we need to complete some of the calculations in NQS, which were left in the form of integrals in Ref. [7]. We have carried out these calculations in Appendix E, correcting some typographical errors in the expressions in Ref. [7], inserting missing color factors, and correcting the overall signs. Our results, from Eqs. (E2), (E8), and (E17), are The integral over k + 1 yields an infrared pole: Note that the sum of diagrams IV and VI is rotationally invariant: The real IR pole comes only from the diagram V: We note that, in NQS, this contribution is dropped because it is regarded as a one-loop correction to an SDC, times the one-loop correction to the LDME. However, in Appendix E 2, we have checked this assignment of the contribution of diagram V and conclude that it is a two-loop IR-divergent contribution to the LDME.
In order to compare our results with those in NQS, we expand our results in Sec. V through order v 2 (order q 2 ) to obtain where, following NQS, we normalized the quark mass m Q to unity.
There are several differences between our results and those in NQS.  117)]. As we have said, this coefficient is regarded in NQS as being UV in origin, but we believe it to be IR in origin. It is plausible that these discrepancies might also be removed once the UVregulator mismatch has been taken into account fully in the NQS calculation, but we have not checked this. In any case, these terms, some of which are dependent, are consistent with NRQCD factorization because, owing to their UV origins, they can be interpreted as one-loop contributions to an SDC.
UV contributions from the UV-regulator mismatch could, in principle, occur in individual contributions from the non-Abelian diagrams in the calculation in NQS. There could be a residual contribution from the mismatch between the double-real-gluon contribution and real-virtual-gluon contribution that comes from the residue of the pole in the virtual-gluon propagator that is labeled as k 1[(k 1 −k 2 ) 2 ] in NQS. There could also be a residual contribution from the residue of the pole in the virtual-gluon propagator that is labeled as k 1[k 2 1 ] in NQS. That residue is discarded in NQS because of the antisymmetry of the integrand under the interchanges P 1 ↔ P 2 , k + 1 ↔ k + 2 , and k 1⊥ ↔ k 2⊥ . However, that symmetry is violated by the UV phase-space regulator. Evidently, such residual contributions, if they are present, cancel in the final result, since we find no UV contributions in the coefficient of the IR pole in our covariant calculations.
Finally, let us consider the two-loop IR-divergent contributions, which, in the Abelian case, come from diagrams C and E in our calculation. The sum of these contributions is equal to the contribution in V in the NQS calculation, and so our result for the two-loop IR-divergent contribution agrees with that in NQS, provided that we interpret V as being a two-loop IR-divergent contribution.

VII. CONCLUSIONS
The central issue in establishing the NRQCD factorization conjecture for inclusive quarkonium production is the question of the universality of the NRQCD long-distance matrix elements (LDMEs). Universality of the LDMEs requires that any IR divergences that they contain be independent of the direction of the Wilson line which makes the LDMEs gauge invariant.
In order to test for a possible dependence on the Wilson-line direction, we have used covariant methods to carry out a two-loop calculation of the NRQCD LDME for a heavy QQ pair in an S-wave, color-octet state to evolve into a color-singlet state. Our calculation provides a check of a previous two-loop calculation by Nayak, Qiu, and Sterman in Refs. [7][8][9], who used light-cone methods to carry out their calculation. Although our results differ from those of Nayak, Qiu, and Sterman in several respects, we find, as did they, that the LDME is independent of the direction of the Wilson line.
One might hope that a covariant calculation would reveal simplifications in comparison with the rather complicated light-cone calculation. That did not prove to be the case in our calculation, probably because, for non-Abelian diagrams, we did not implement the unitarity cancellations between real-and virtual-gluon corrections at the integrand level.
Consequently, we had to deal with poles in the dimensional-regularization parameter of orders 1/ 3 , 1/ 2 , which canceled, in order to obtain a final result of order 1/ . Furthermore, in the case of the order-1/ contribution, our intermediate expressions contained hundreds of terms involving dilogarithms and trilogarithms, which ultimately canceled to yield a very simple expression.
Our results for the non-Abelian diagrams F agree with those in Refs. [7][8][9], up to an overall sign, and are independent of the Wilson line direction. Our results for the Abelian diagrams A and B that involve two gluon connections on the Wilson line do not contribute an IR pole, in agreement with the calculations in Ref. [7]. However, our results for the Abelian diagrams C, D, and E that involve one gluon connection on the Wilson line differ in some respects from those in Ref. [7], although both our results and those in Ref. [7] are consistent with the NRQCD factorization conjecture. In our case, we find dependences on the Wilsonline direction in some contributions. However, these contributions can be factored into a one-loop contribution to a short-distance coefficient times the one-loop contribution to the LDME, the latter of which is independent of the Wilson-line direction.
We have identified one source of the differences between our calculation and that in Ref. [7]. In Ref. [7], a double-real-gluon contribution is canceled against a particular realvirtual-gluon contribution that comes from the residue of the pole in the virtual-gluon propagator. That cancellation is exact in the IR limit, but leaves residual UV contributions from the virtual-gluon loop because, unlike the corresponding real-gluon loop, the virtual-gluon loop is not constrained by phase space in the UV. These residual contributions account for UV double poles that are present in our result, but not in the result in Ref. [7]. They may also account for UV single poles and UV constant terms that are present in our results, but not in the results in Ref. [7], although we have not checked this explicitly.
In our result, we find a two-loop IR-divergent contribution to the LDME that arises from the Abelian diagrams C and E. This contribution is not present in the result for the Abelian diagrams in Ref. [7]. However, if we reinterpret the contribution from diagram V in Ref. [7] as a two-loop IR-divergent contribution to the LDME, rather than as a one-loop contribution to the SDC times the one-loop contribution to the LDME, as is done in Ref. [7], then we find agreement with Ref. [7] for the total two-loop IR-divergent contribution to the LDME.
While our result confirms the NRQCD factorization conjecture through two loops, the principle that underlies that result has not emerged, and its elucidation remains a challenge for future work. The U.S. Government retains for itself, and others acting on its behalf, a paid-up nonexclusive, irrevocable worldwide license in said article to reproduce, prepare derivative works, distribute copies to the public, and perform publicly and display publicly, by or on behalf of the Government.

Appendix A: Useful formulas
In this appendix, we compile some formulas that are useful in our calculation.
The Feynman parametrization that we use for parameters that run from 0 to 1 is given The Feynman (Schwinger) parametrization that we use for parameters that run from 0 to ∞ is given by We derive the following integral formula by calculating first with real values of the parameters a, b, α, and β and then constructing an analytic continuation in those parameters that is consistent with the original integral: where a, b, α, and β are complex numbers. This formula is used for the phase-space integrations and the parameter integrations in our calculations.
In computing the k 1 integrations of the Abelian diagrams, we encounter a parameter integration of the following form: where K 1 and K 2 are k 2 -dependent Lorentz products. We perform the λ 1 integration in Eq. (A4) by making use of the Mellin-Barnes representation with The initial contour for the z integration runs parallel to the imaginary axis with Re(z) = z 0 and separates the left poles and right poles of the Γ functions, where the left (right) poles have a positive (negative) sign in the argument of the Γ function. We analytically continue the function in Eq. (A5) in 1 , 2 , and 3 to a common value near zero, deforming the contour for the z integration so that the poles in the Γ functions never cross the contour. The result is a curved contour, in the style of the method of Smirnov [12], that still separates the left and right poles. If we close the contour at ∞ in the right half of the z plane, we pick up only the right poles. The result is an asymptotic expansion in powers of K 1 /K 2 : In our applications of this parameter integral, K 1 /K 2 is small, and we can retain the contribution of leading order in K 1 /K 2 , which comes from the n = 0 terms. This leading contribution is given below.
We also make use of the following elementary integrals, which are valid for a ≥ 0: We simplify a number of expressions in our calculations by applying the following polylogarithm identities: and where Li 3 (x) is the trilogarithm, which is defined by In this appendix, we derive the phase-space integration formulas that are given in Eq. (18).
First, let us consider where p is an arbitrary real four-vector and M 2 is a real number.
In order to simplify the calculation, we initially choose a coordinate frame in which p ⊥ = 0. Then, re-write the result of this calculation in a form that is manifestly rotationally invariant. In the frame in which p ⊥ = 0, I 1 is given, in light-cone coordinates, by Note that, because of the δ-function and the θ-function, we have k + ≥ 0 and k − ≥ 0.
Using the δ-function to integrate over k − , we find that In the last line we have introduced the definitions and λ ≡ k + |p + | , where p + = sign(p + )|p + |. Next, let us perform the k ⊥ integration. Carrying out the angular integrations and making the change of variables z = |k ⊥ | 2 , we find that where in the last equality, we have carried out the z integration by making use of the integral formula in Eq. (A3) and have made use of the n-dimensional solid angle .

(B6)
Substituting this result into Eq. (B3), we have We carry out the λ integration by using Eq. (A3) once again. The result is Now we restore rotational invariance in our result by making the replacements p + → Clearly, this expression is rotationally invariant. It is also invariant under boosts because the magnitude of the first two denominator factors is |p 2 | D 2 −1 and the signs of the arguments in the first two denominator factors are also boost invariant.
Next, let us consider where q is an arbitrary four-vector.
Again, we initially choose a coordinate frame in which p ⊥ = 0 and, then, re-write the result of the calculation in a form that is manifestly rotationally invariant. In the frame in which p ⊥ = 0, q α I α 2 is given, in light-cone coordinates, by Using the δ-function to carry out the k − integration and performing the change of variables that is below Eq. (B3), we obtain where we have dropped terms that are odd in k ⊥ . The k ⊥ integration yields Substituting this result into Eq. (B12) and carrying out the remaining λ integration by making use of Eq. (A3), we find that Now, we restore rotational invariance in our result by making the replacements p + → (B15) Appendix C: IR finiteness of integrations in B i , C i , and D i In this appendix, we demonstrate the IR finiteness of the k 1 integration with k 2 fixed in B 1a + B 2 , C plus hermitian conjugate, and D plus hermitian conjugate.
1. IR finiteness of the k 1 integration in B 1a + B 2 and the k 2 integration in B 1a + B 3 We extract the k 1 integration in B 1a and B 2 and extract the k 2 integration in B 1b and B 3 , We can obtain B P 1 P 2 1b Reg and B P 1 P 2

3
Reg from the expressions for B P 1 P 2 1a Reg and B P 1 P 2 2 Reg , respectively, by taking their hermitian conjugates and making the substitutions k 2 ↔ k 1 and P 1 ↔ P 2 in the integrands. Therefore, it follows that B P 2 P 1 1b Reg + B P 2 P 1 3 Reg is the hermitian conjugate of B P 1 P 2 1a Reg + B P 1 P 2 2 Reg .
Applying Feynman parameters to Eq. (C2), we find that Let us compute B 1a Reg first. Performing the k 1 phase-space integration and the z integration, we obtain Separating the contributions that correspond to the first and second terms in braces in Eq. (C4), we have where Let us consider the λ 1 integration in B P 1 P 2 1aa . We can rotate the contour of integration counterclockwise by an angle of π without encountering any singularities. Therefore, by Cauchy's theorem, we have Next, let us consider B 2 Reg in Eq. (C3). Carrying out the k 1 virtual integration and λ 1 integration, we obtain where, in the result, we have made a change of variables λ 2 → 1/λ 1 . Then, Therefore, we find that Carrying out the x and y integrations, we obtain Note that the first two terms in the numerator cancel as λ 1 → 0, and so they result in a single IR pole. The λ 1 integration of the remaining numerator term can be carried out by making use of the integration formula in Eq. (A3), and it results in a real double pole. Since we can conclude that the double IR pole from the k 1 integration in B P 1 P 2 2 Reg is canceled by the IR pole from the k 1 integration in B P 1 P 2 1aa Reg . The residual imaginary pole will be canceled by the corresponding pole in the hermitian-conjugate contribution B P 2 P 1 1b Reg + B P 2 P 1 3 Reg . The remaining contribution to B P 1 P 2 is B P 1 P 2 1ab Reg [Eq. (C6)]. Let us consider the parameter integrations of B P 1 P 2 1ab . Carrying out the y integration, we obtain In the last line, we have rotated the λ 1 contour of integration counterclockwise by an angle of π, using the fact that one does not encounter any singularities in carrying out that rotation.
Then, we find that Carrying out the x integration, we obtain The λ 1 integration does not produce a pole in . Therefore, we can expand the integrand in powers of . The order-0 term vanishes and the order-1 term cancels the 1/ in Eq. (C15).
Hence, we conclude that the k 1 integration with k 2 fixed in B P 1 P 2 1ab Reg is IR finite.
Therefore, we conclude that the k 1 integration with k 2 fixed in B 1a + B 2 is IR finite, from which it follows that the k 2 integration with k 1 fixed in B 1b + B 3 is IR finite.
2. IR finiteness of the k 1 integration in C 1 + C 2 We extract the k 1 integrations in C P 1 P 2 1 and C P 1 P 2 2 , writing , . (C17) Applying Feynman parameters to Eq. (C17), we find that Let us compute C 1 Reg first. Performing the k 1 phase-space integration and λ 2 and x parameter integrations, we obtain Since the λ 1 integration converges at both 0 and ∞, we can expand the integrand as a series in and carry out the λ 1 integration term by term to obtain Next, let us consider C 2 Reg . Carrying out the virtual k 1 integration and the λ 1 integration, using Eq. (A3), we obtain where, in the result, we have made a change of variables λ 2 → 1/λ 1 . Then, carrying out the remaining λ 1 and x integrations, we obtain Inserting the results in Eqs. (C20) and (C22) into Eq. (C16), we obtain Hence, the k 1 integration in C 1 + C 2 is IR finite.
3. IR finiteness of the k 1 integration in the D 1 + D 2 We extract the k 1 integrations in D P 1 P 2 1 Reg and D P 1 P 2 2 Reg , writing , . (C25) Applying Feynman parameters to Eq. (C25), we find that Let us compute D 1 Reg first. Carrying out the k 1 phase space integration and λ 2 and x parameter integrations, we obtain The two terms in brackets cancel as λ 1 → 0, and so there is no divergence in the λ 1 integration from that limit. However, a divergence remains in the first term from the limit λ 1 → ∞. This divergence will eventually be canceled by the corresponding divergence in D 2 Reg . We separate the divergences in the first term by inserting where the first term on the right side of Eq. (C28) yields an IR pole, and the second term gives a UV pole. Then, we can write Eq. (C27) as follows: Now the λ 1 integration of the first term in Eq. (C29) converges in the UV and the IR, and so, for it, we can expand the integrand in a series in . The integral of the second term in Eq. (C29) can be expressed as a beta function and gives an IR pole. Then, expanding D 1 Reg in a series in , we obtain Next, let us consider D 2 Reg . Carrying out the virtual k 1 integration and the λ 1 integration using Eq. (A3), we obtain where, in the result, we have made a change of variables λ 2 → 1/λ 1 . Then, carrying out the x integration, we find that As in the case of D 1 Reg , the divergence in Eq. (C32) that appears in the integration over λ 1 as λ 1 → 0 cancels between the two terms in brackets. However, there remains a divergence that appears as λ 1 → ∞. We separate the UV and IR divergences in the first term by inserting 1 = 1 to obtain In Eq. (C34), the integration of the first term over λ 1 is finite, and so, for this term, we can expand the integrand in a series in . The integration of the second term λ 1 can be expressed as a beta function and yields an IR pole. Then, expanding D 2 Reg as a series in , we obtain Inserting the results in Eqs. (C30) and (C35) into Eq. (C24), we find that We see that the real double and single poles cancel between D 1 and D 2 , while the single imaginary pole in D 2 cancels when we add the hermitian-conjugate contribution. Hence, the k 1 integration with k 2 fixed in D 1 + D 2 plus hermitian-conjugate is IR finite.
Appendix D: The one-loop IR contribution to the LDME The only possible one-loop contribution to the matrix element in Eq. (3) is shown in Fig. 9. The superscripts P j P k specify the gluon connections to the Q orQ lines on either side of the cut.
The color factor of the O diagram is given by The expression for the diagram O with our standard phase-space regulators in Eq. (17) is given by .
Combining the denominators by using Feynman parameters and carrying out the k 2 phasespace integration by using Eq. (18), we find that where we have made change of variables λ 3 → ξλ 2 , carried out λ 2 integration in terms of the beta function, used the definitions of A and B that are given in Eq. (46), and carried out the ξ integration by using Eq. (A9).
Then, taking into account the color factor in Eq. (D1) and summing over the gluon attachments, we obtain If we expand Eq. (D4) to leading order in q 2 , then it agrees with the expression in Eq. (38) of Ref. [7], aside from the color factor, which has been dropped in Ref. [7].

Appendix E: NQS Calculations
In this appendix, we correct some signs and typographical errors and supply color factors in the expressions for diagrams IV, V, and VI in Ref. [7]. We mark these changes relative to the expressions in Ref. [7] with double brackets. We also complete the calculations of the integrals in the expression for diagram V.

Diagram IV
The contribution of diagram IV is given by which agrees with Eq. (51) of Ref. [7] up to the missing color factor f 2 abc 4Nc and an overall sign. Performing the remaining integrations, except for the k + 1 integration, we obtain The contribution of diagram V is given by which agrees with Eq. (54) of Ref. [7] up to the corrections that are marked with the double brackets. Using the rescaling of integration variables in Ref. [7] k + 2 = (k + 1 )y, k i⊥ = ( we can rewrite Eq. (E3) as follows: which agrees with Eq. (55) of Ref. [7], up to the corrections that are marked with double brackets. J V (κ 1 ) is defined by Performing the κ 2 integration, we can write Eq. (E6) as which agrees with Eq. (57) of Ref. [7]. Performing the remaining y and κ 1 integrations, we In Ref. [7], the complete factor that multiplies the k + 1 integral in Eq. (E8) was considered to be UV in nature and, therefore, to be a contribution to an SDC. Consequently, the entire contribution in Eq. (E8) was taken to be a one-loop contribution to an SDC times the one-loop contribution to the LDME and was dropped in Ref. [7]. As we will now show, the real contribution that comes from the product of the UV pole and the term −π 2 in Eq. (E8) should actually be considered to be IR in nature. We show this by completing the k 2 integration in V before carrying out the k 1 integration.
First, we re-write Eq. (E3) as where δ + (k 2 1 ) ≡ δ(k 2 1 )θ(k + 1 + k − 1 ), and K a 2 and K b 2 are Making the changes of variables k + 2 → k + 2 + k + 1 and k 2⊥ → k 2⊥ + k 1⊥ and carrying out the k 2⊥ integration, we obtain We can carry out the k + 2 integration by making use of the following formulas: ∞ −∞ dy 1 (y + iε) y(y + k + 1 + k − 1 ) + iε Then, we find that We can see that these expressions have an IR sensitivity because of the denominator factors (k + 1 + k − 1 ) 1+2 and (k + 1 + k − 1 ) 2 , which becomes singular when k 1 goes to zero, as happens for the IR divergence in the k 1 integral. These factors from the k 2 integration affect the strength of the IR pole in the k 1 integration. The only part of Eqs. (E12) that can be considered to be UV in nature is the pure UV pole in the expression for K b 2 , for which the denominator factor (k + 1 + k − 1 ) 2 is set to unity. However, this pure UV pole cancels when we add the hermitian-conjugate contribution. All of the other contributions in Eqs. (E12) are IR in nature. This is true, in particular, of the only real contribution, which comes from the product of the UV pole and the −π 2 term in K b 2 . Having demonstrated the IR nature of the real part of the result, we complete the calculation by carrying out the k 1 integration. Using our results from the k 2 integration, we We carry out the k 1 integration by making use of the δ + -function and make the change of variables k 1⊥ → √ 2k + 1 κ 1 and carry out the κ 1 integration. The result is which is in agreement with Eq. (E8).

Diagram VI
The contribution of diagram VI is given by In this section, we demonstrate how the UV double poles in C and D arise in an NQSstyle calculation. As we have mentioned, the source of the double poles is a mismatch between a real-gluon contribution, which is subject to a phase-space UV regulator, and the contribution from the residue of a pole in the propagator of the corresponding virtual gluon, which is unregulated in the UV.

Diagram IV
Let us first consider the contribution of diagram IV in Ref. [7]. In terms of light-cone variables the real-virtual contribution is [7] IV A = 4ig 4 Note that the last denominator can be rewritten as .

(F2)
We close the k − 2 contour in the upper half-plane and pick up the residue of the pole at (k 2 − k 1 ) 2 = 0. The location of the pole is and the residue is nonzero only when k + 2 > k + 1 . The residue is Carrying out the k + 2 , k 1⊥ , and k 2⊥ integrations, we obtain We note that the poles in this expression are purely IR in origin. That is, the integrations in IV A (k 2 −k 1 ) 2 are UV convergent. The real-real contribution IV B is precisely the negative of the real-virtual contribution IV A (k 2 −k 1 ) 2 , except that the k 2 integration in IV B is cut off by a UV regulator. Since the integrations in both IV A (k 2 −k 1 ) 2 and IV B are UV convergent, they cancel, up to terms that are suppressed by inverse powers of the UV-regulator scale Λ.
Hence, the diagrams IV do not yield any UV double poles.

Diagram V
Next, let us consider the contribution of diagram V in Ref. [7]. In terms of light-cone variables the real-virtual contribution is [7] Again, we close the k − 2 contour in the upper half-plane and pick up the residue of the pole at (k 2 − k 1 ) 2 = 0. The location of the pole is given in Eq. (F3), and the residue is nonzero only when k + 2 > k + 1 . The residue is After the changes of variables for i = 1 and 2, we can perform the κ 2 integration to obtain The contribution of the corresponding real diagram is Here, we have imposed a UV cutoff on the k + 2 integration. While this cutoff is different from our standard UV phase-space regulator, we expect that this will not affect the coefficient of

Diagram VI
Finally, let us consider the contribution of diagram VI in Ref. [7]. In terms of light-cone variables the real-virtual contribution is [7] V IA = i2 5/2 g 4 Again, we close the k − 2 contour in the upper half-plane and pick up the residue of the pole at (k 2 − k 1 ) 2 = 0. The location of the pole is given in Eq. (F3), and the residue is nonzero only when k + 2 > k + 1 . The residue is Making the changes of variables we obtain (F15) It can be seen easily that the κ 1 and κ 2 integrations cannot give UV poles, and so the diagram VI does not contribute any double UV poles.