Next-to-leading power two-loop soft functions for the Drell-Yan process at threshold

We calculate the generalized soft functions at $\mathcal{O}(\alpha_s^2)$ at next-to-leading power accuracy for the Drell-Yan process at threshold. The operator definitions of these objects contain explicit insertions of soft gauge and matter fields, giving rise to a dependence on additional convolution variables with respect to the leading power result. These soft functions constitute the last missing ingredient for the validation of the bare factorization theorem to NNLO accuracy. We carry out the calculations by reducing the soft squared amplitudes into a set of canonical master integrals and we employ the method of differential equations to evaluate them. We retain the exact $d$-dimensional dependence of the convolution variables at the integration boundaries in order to regulate the fixed-order convolution integrals. After combining our soft functions with the relevant collinear functions, we perform checks of the results at the cross-section level against the literature and expansion-by-regions calculations, at NNLO and partly at N$^3$LO, finding agreement.

1 Introduction contributions calculated in [12] and collected in App. C. In the last step, we carry out the remaining UV renormalization and remove the initial state collinear singularities at cross section level.
To the best of our knowledge, this is the first time that soft functions at NLP are evaluated to O(α 2 s ). However, in the present work, we do not analyze the UV renormalization and the RG evolution of the soft functions since any such procedure requires an expansion around d → 4 which leads to the appearance of divergent convolutions integrals preventing naive renormalization attempts. Our intention is to provide more information about the higher order structure of these soft functions which could give a hint towards the solution of the divergent convolution problem, at least in the present case. However, it should be noted that at NLP accuracy, calculations and studies of the renormalization and evolution properties of the soft function needed for the h → γγ decay process were recently presented at O(α s ) in [34,35]. In the context of the q Tsubtraction method, calculations at fixed-order accuracy, which required the evaluation of several new integrals, have been recently carried out at NLP in QCD [36] without separating the different regions.
The paper is organized as follows. In Sec. 2 we review the structure of the factorized cross section which is one of the main results of [12]. In Sec. 3 we describe in detail the calculation of the two-loop soft functions. In particular, we discuss the evaluation of the canonical master integrals using the differential equation method. In Sec. 4 we calculate the convolution integrals of the soft functions with the corresponding collinear functions [12] which allows us to carry out a series of checks at the cross-section level against the literature, [30,31,33], and against expansion-by-regions calculations, both at NNLO and available results at N 3 LO [32]. In App. A we list the analytic results for the collinear functions calculated in [12], and in App. B we provide expressions for the relevant two-parton matrix elements used in this calculation. Useful cross-section level results from [12] are collected in App. C. Finally, App. D contains the expressions for the relevant Altarelli-Parisi splitting functions.

Factorization near threshold
In this section, we review the structure of the NLP factorization theorem for the DY process in the threshold region [12], and we remind the reader of the operatorial definitions of the NLP soft functions which are relevant for this work.
We consider the parton diagonal channel of the DY process, qq → γ * [→ ¯ ] + X in the kinematic region z = Q 2 /ŝ → 1, whereŝ = x a x b s is the partonic centre-of-mass energy squared, x a , x b are the momentum fractions of the partons inside the incoming hadrons and Q 2 is the invariant mass squared of the lepton pair. Up to NLP in the threshold expansion, the cross-section differential in Q 2 assumes the following form where f a/A (x a ) and f b/B (x b ) are the usual parton distribution functions (PDFs), and σ(z) with superscripts LP and NLP are the leading power and the next-to-leading power partonic cross sections, respectively. Since we only focus on the qq-channel, we omit the indices a, b in the following. The LP partonic cross section factorizes into a product of two functionsσ the hard function H(Q 2 ) and the soft function 3) The LP position-space soft function is a vacuum matrix element of soft Wilson lines 4 [38] S 0 ( We now turn our attention to the NLP part of the factorization formula, which is understood to be formally valid only in d-dimensions, before renormalization. To facilitate comparison with literature, we define the quantity ∆ related to the partonic cross section as follows (2.6) The NLP partonic cross section receives contributions from power corrections to the phase-space, the so-called kinematic corrections, and from insertions of subleading power Lagrangian terms in time-ordered product operators, referred to as dynamical corrections. We have The ∆ kin NLP (z) term has been presented in Eq. (5.11) of [12] at NNLO. In this work, we focus on the calculation of the generalized soft functions which appear in the factorization formula in the ∆ dyn NLP (z) contribution. The result for ∆ dyn NLP (z) takes the following form [12] ∆ dyn where Ω = Q(1 − z). In the equation above, C A0,A0 is the hard matching coefficient of the LP SCET current for the DY process. The J i are the collinear functions and the S i represent the generalized soft functions in momentum space, defined as 9) in terms of the position-space multi-local soft functions, S i (x 0 ; {z j− }). At NLP these are given by In the above definitions, µ, ν are the Lorentz indices, σ, λ are the Dirac indices, and A, B and a, b, f, g, h are adjoint and fundamental colour indices, respectively. The B ± (q + ) field is a soft building block formed by a soft covariant derivative (soft quark field) and soft Wilson lines The soft functions in Eqs. (2.10) -(2.14) are the fundamental objects of interest in this work. Thus far, only partial results for these objects have been reported in the literature.
The O(α s ) result for S 1 , expanded in , was given in [6], and with complete d-dimensional dependence in [12]. At O(α 2 s ), results for virtual-real soft diagrams have been presented in [12]. It was found that only the S 1 soft function receives such contributions. In this article, we complete the calculations of the bare soft functions at O(α 2 s ) by considering all the diagrams with two-real soft parton emissions.
A quick inspection of the soft functions in (2.10) -(2.14) reveals that S 1 and S 3 are conveniently defined as scalar objects. The remaining three functions, S 2 , S 4 and S 5 , contain instead a non-trivial dependence on Lorentz, Dirac, and adjoint and fundamental colour indices. These indices are contracted with the corresponding indices carried by the respective collinear functions, which are reported for convenience in App. A. We prefer to work with scalar objects and, for this reason, we absorb the colour, spin and Lorentz structures of the multi-local collinear functions into their corresponding soft functions S 4 and S 5 . For example, making use of (A.8), the part of the factorization formula at O(α 2 s ) which depends on S 4 is given by, 5 Focusing on the collinear and soft functions, we now redefine 4 (Ω; ω 1 , ω 2 ), (2.17) such that the new J are scalar functions. The S 5;bf gh,σλ soft function in Eq. (2.14) is redefined in an analogous way by factoring out the same scalar collinear function such that J The S 2;µν soft function is anti-symmetric under the exchange of the Lorentz indices µ, ν. Since this is a vacuum matrix element, it must be proportional to the epsilon tensor µν ⊥ which is the only anti-symmetric structure which can carry two transverse Lorentz indices. However, its appearance is forbidden due to parity conservation of QCD. Indeed, we checked by direct calculation that S 2;µν vanishes at O(α 2 s ).

Two-loop soft functions
We proceed with the main focus of this work by providing the calculation details of the double real emission corrections to the soft functions defined in (2.10) -(2.14). Techniques for solving integrals which appear in calculations of LP soft functions at NNLO have been developed over the years and several examples exist in the SCET literature.
In particular, results for the exclusive soft function relevant for small transverse momentum resummation in DY was obtained in [39], the soft function for the production of an electroweak boson at large transverse momentum was computed in [40], and the calculation of the soft function relevant for boosted top-quark pair production was presented in [41]. However, these methods are insufficient for the NNLO calculation of the NLP soft functions containing dependence on additional convolution variables. Therefore, we apply more advanced techniques developed for fixed order calculations. Similar methods were used in [42] to calculate the NNLO soft function for top quark pair production at threshold. The strategy is straightforward: we first obtain the squared amplitudes at O(α 2 s ). Subsequently we use LiteRed [27,28] to reduce such expressions to a linear combination of master integrals (MIs), and finally, we compute the necessary MIs by employing the differential equation method. Each of these steps is expanded upon in the following sections.

Reduction to master integrals
The two-loop expressions for the soft functions S 1 and S 3 are directly obtained from their matrix element definitions. The soft functions S 4 and S 5 are computed after the redefinition made in (2.17) by employing the NLP Feynman rules given in [20]. The expressions for the squared matrix elements of all the soft functions are collected in the ancillary.pdf file and they correspond to the diagrams in Figures 1 -4. We use k 1 and k 2 to label the momenta of the partons crossing the cut.

Topologies
The calculation of the double real emission corrections to the soft functions includes two types of phase space constraints. The first constrains the total energy radiated into the final state, enforced by the δ(Ω − 2E X ) condition, where E X is the total radiation energy. The second type instead constrains specific light-cone components of the soft parton momenta. Indeed, the NLP soft functions are also differential in ω, or ω 1 and ω 2 , which are the convolution variables that connect the soft functions to their corresponding collinear functions, as prescribed by (2.8). We find that up to O(α 2 s ), three different constraints of the second type are possible. Namely, the integrands of the soft functions depend on δ(ω − n − k 1 ), 6 . These constraints, along with the on-shell cut propagators conditions δ(k 2 1 ) and δ(k 2 2 ), are set in the LiteRed program. We now define the auxiliary topologies, A, B, and C, which implement the δ(ω − n − k 1 ) constraint and only differ among themselves by the choice of one propagator. Topology A is defined by the following set of seven propagators where the last four propagators are cut propagators. This means that and equivalent relations hold for P 5 , P 6 and P 7 . Topology B is obtained starting from the list of propagators in (3.1) and replacing the single propagator Similarly, topology C requires the substitution P 3 → n + (k 1 + k 2 ). The integrals which appear in our calculations are written aŝ where the index T indicates the specific topology T ∈ {A, B, C, . . .}, and can be expressed as a linear combination of the independent MIs. It turns out that all the MIs for the topologies A and C are a subset of the MIs for the topology B. This is not surprising since the three topologies share most of the propagators. In particular, we find the following five MIs for topology B (1, 0, 0, 1, 1, 1, 1),Î 4 (Ω, ω) ≡Î B (1, 1, 0, 1, 1, 1, 1), where the integralÎ 1 (Ω, ω) represents the phase space integral. The squared matrix elements with the δ(ω − n − k 1 − n − k 2 ) constraint require four additional topologies to be reduced. The topology D is defined by the list of propagators where P 4 to P 7 are cut. Topology E is obtained from (3.5) by replacing P 3 → n − k 1 , topology F by substituting P 2 → n + k 1 , and the topology G by exchanging both P 2 → n + k 1 and P 3 → n − k 1 . After reduction, we find that only two additional MIs appear for the set of topologies which implement the constraint δ(ω − n − k 1 − n − k 2 ): Two additional topologies are needed to reduce the integral expressions with a double constraint given by δ(ω 1 − n − k 1 ) δ(ω 2 − n − k 2 ). We define the H topology as where only the first two propagators of the list remain uncut. The topology I is related to H by a replacement of the second propagator P 2 → n + k 1 . Only one MI is found for these last topologiesÎ In total we find eight new MIs which need to be computed to evaluate the O(α 2 s ) corrections to the NLP soft functions.

Results after reduction
The integrals belonging to each of the topologies defined above can be reduced by employing the program LiteRed. In this subsection we present the results for the soft functions in terms of linear combinations of MIs. The soft function S 1 carries an additional 2r0v superscript since it is the only one that also receives a virtual-real (superscript 1r1v) contribution. In the following expressions we omit for simplicity the Ω and ω (ω 1 ,ω 2 ) dependence in the MIs and we find where = (4 − d)/2. The S 3 soft function has the following form in terms of MIs and S

Master integrals
We begin by describing the calculation of the MIs for the topology B given in (3.4). Starting from those expressions, it is convenient to make the variable change ω → r Ω, and redefine the MIs by factoring out their mass dimensions in Ω where the prime integrals only depend on the variable r. We use Canonica [43] to guide us in finding the canonical basis of MIs [29]. In particular, this is achieved by the following transformations where the canonical integrals are the ones without the prime. The system of differential equations for the vector of integrals I(r) ≡ I 1 (r), I 2 (r), I 3 (r), I 4 (r), I 5 (r) is given by where (3.16) We notice from the structure of A(r) that the integral I 2 (r) only couples to I 1 (r). The alphabet simply reads {1 − r, r}. The A(r) matrix in Eq. (3.16) is lower triangular and can be solved iteratively. The integral I 1 (r), which is the starting integral of our system of equations, is obtained by direct integration and we find We notice that I 2 (r) and I 3 (r) satisfy the same differential equation. Hence, they will lead to identical results. Starting from the result for I 1 (r), it is possible to compute I 3 (r) (and I 2 (r) = I 3 (r)) by solving the differential equation. The result is The -dependent constant C 3 ( ) can be fixed by requiring that the integral of I 3 (r) in the range r ∈ [0, 1] is equal to the parent integral where the δ(ω − n − k 1 ) constraint is removed. This latter integral is easily evaluated by direct integration. Specifically, we require that which fixes the -dependent constant of I 3 (r) to be We now focus on the integral I 4 (r) that satisfies a differential equation which involves both I 1 (r) and I 3 (r), as dictated by Eq. (3.15). Its solution reads (3.21) The most complicated part of the integration concerns the hypergeometric function which we rewrite using its integral representation Then, we make the variable transformation r → 1 + R (r − 1) and integrate over the range R ∈ [0, 1]. Finally, we carry out the integration over t and arrive at the result which retains exact d-dimensional dependence. In the above equation the integration constant has already been fixed following the same procedure as used for the integral I 3 (r). We checked that the resulting expression for I 4 (r) satisfies its initial differential equation. Finally, we consider the last and most difficult integral of topology B: I 5 (r). The differential equation for I 5 (r) reads 24) and the solution has the following structure We do not write explicitly the function f I 5 (r , ) since it is too lengthy, but we know its exact expression in d-dimensions. Unfortunately, we are not able to directly integrate f I 5 (r , ) in d-dimensions since it involves a 3 F 2 hypergeometric function. Nevertheless, we devise a technique which allows us to retain the dependence on r −2 terms, after r integration, which is the relevant information that is needed to regularize the convolution integrals. Indeed these are the potential problematic contributions since the division by r in the transformation to the non-canonical basis (see the last term of (3.14)) will generate delta terms and plus distributions, after expansion, for the non-canonical integral I 5 (r). Therefore, we need to treat these terms with care. We follow the strategy of expanding the function f I 5 (r , ) in the limit r → 0 (up to finite order in r ) and add and subtract this term in the following way dr lim where the -limit in the above equation means that one needs to perform the → 0 expansion up to the relevant order. In Eq. (3.26) we split the integral in two terms, the first which we are able to integrate in d-dimensions exactly and the second which we need to -expand before integration. For the first term, we find The second term in (3.26) contains terms that are non-singular in the r → 0 limit and can be expressed in terms of standard HPLs [44] as follows We now have all the ingredients to construct the non-canonical integral I 5 (r) by using the last equation of (3.14) which combines the canonical integrals I 4 (r) and I 5 (r). The expression retains the exact d-dependence for terms of the type r −1−2 and (1 − r) −1−4 which makes it suitable to be convoluted with O(α s ) collinear functions avoiding undefined convolutions at fixed-order accuracy. In total we find 12 π 2 − 6 + 30 12 6 ζ 3 + 5 + 1 6 (r − 1)r 6 (6 − 7r)Li 2 (r) + 48 (r − 1) ln 2 (1 − r) + r 3 ln 2 r + π 2 + 6 ln r − 12 ln(1 − r) (r + 1) ln r + 2(r − 1) θ(r)θ(1 − r) + O( ) .

(3.29)
One can also expand the boundary singular terms in → 0 by using the relation (3.30) to find In order to reproduce the cross section we need to integrate our results over r in the range r ∈ [0, 1] and we obtain The integration constants are fixed similarly to I 3 and I 4 . It turns out that they are zero up to the finite order in . We still need to discuss the calculation of the two MIs in Eq. (3.6) which implement the constraint δ(ω − n − k 1 − n − k 2 ) and the last MI in Eq. (3.8) with the constraint δ(ω 1 − n − k 1 ) δ(ω 2 − n − k 2 ). These integrals are carried out by direct integration in a straightforward way. For completeness we report the results beloŵ where r 1 = ω 1 /Ω and r 2 = ω 2 /Ω.

Results
In this subsection we collect the final expressions for the bare soft functions which constitute the main results of our work. We retain the relevant d-dimensional dependence at the integration boundaries to avoid divergent convolutions when combining the soft functions with collinear functions [12] to fixed-order accuracy. We refrain from expanding our results in d → 4 since a consistent procedure for the renormalization of the soft functions beyond LL accuracy is not yet available in the literature. Starting from the reduced result for S (Ω, ω) here due to its length, but it is possible to easily reconstruct it from the information provided in the two subsections above. With an analogous procedure we obtain the results for the S 3 , S 4 and S 5 soft functions In the above expressions the strong coupling constant is understood to be the renormalized coupling α s ≡ α s (µ) in the MS scheme obtained via the relation Z α α s µ 2 = (4πe −γ E ) α 0 s , where α 0 s is the bare coupling constant and Z α = 1 − β 0 α s /(4π ) with β 0 = 11 3 C A − 4 6 n f . The discussion of the renormalization procedure for the soft functions is outside the scope of this paper due to the aforementioned divergent convolution problem. However, one should keep in mind that, according to coupling renormalization, contributions proportional to Z α S (1) 1 will appear at O(α 2 s ) and must be taken into account before operator renormalization.

Next-to-next-to-leading order: soft contributions
After presenting the main results of this work in the section above, we now proceed with their validation through comparisons against cross-section level results which are available in the literature. The soft functions calculated in this work carry a dependence on the convolution variables ω or ω 1 , ω 2 . In order to evaluate their contributions to the cross section, one needs to carry out the convolution integrals of the ω-dependent soft functions with their respective collinear functions.
We begin with the contribution of S 1 by considering the relevant part of the factorization formula in (2.8) with the variable transformation r = ω/Ω. We have (Ω, r) . (Ω, r) soft function and the tree-level collinear function in (A.5), then integrating over the convolution variable r one finds to all orders in : 7 Setting the soft scale to Ω = Q(1 − z), the scale µ = Q, and finally expanding in we arrive at the following expression We recall that the complete contribution to S 1 proportional to C F C A comprises an additional term, which stems from diagrams involving virtual-real corrections. Such contribution can be found in Eq. (5.10) of [12], and reads It is interesting to notice that the leading logarithmic contribution proportional to C F C A cancels at cross section level when summing the double real, Eq.
In addition to S 1 , we need to take into account the contributions due to the other soft functions S 3 , S 4 , and S 5 , which can be obtained integrating over the convolution variables the corresponding term in (2.8), similarly to what written in (4.1). S 3 and S 4 read Expanding in with Ω = Q(1 − z) and µ = Q we obtain After convolution with the corresponding collinear function, S 3 and S 4 gives opposite contributions to the partonic cross section, such that they effectively cancel each other at this order. The last contribution to the partonic cross section is given by the term involving S 5 , and reads Setting Ω = Q(1 − z), µ = Q and expanding in we find We point out that the contributions to the cross section due to the soft functions starting at O(α 2 s ), namely S 3 , S 4 , and S 5 , in (4.7) and (4.9), do not contain leading logarithmic terms. This confirms an assumption made in [6], where it was claimed that a logarithmically enhanced off diagonal mixing of these soft functions with the single gluon soft function is not possible.
Results concerning the calculation of soft gluon contribution to the partonic Drell-Yan cross section within a diagrammatic approach have been presented in [31]. The contribution due to double real soft radiation is provided in Eq. (5.2) of [31] and contains the contribution due to ∆ dyn (2)2r0v NLP−soft presented here in Eq. (4.3). However, a direct comparison is not straightforward, because the expression in Eq. (5.2) of [31] contains also NLP corrections due to the expansion of the phase space from the integration of the LP matrix element squared. In the present approach, these are a part of the kinematic correction in Eq. (2.7) discussed in [12]. Similarly, the contribution due to virtual-real soft radiation given in Eq. (4.4) of this paper is included in Eq. (4.6) of [31]. However, the two contributions cannot be compared directly, as Eq. (4.6) of [31] contains also the contribution due to hard and collinear loops. The problem can be overcome by comparing with the individual terms giving rise to Eq. (4.6) and (5.2) of [31], provided by one of us, and we confirm that the whole contribution to the cross section due to the soft function S 1 in (4.5) agrees with the cross-section level diagrammatic calculation of [31]. Moreover, the contribution due to S 5 has not been considered in [31], and we validate its expression against an in-house calculation performed with the expansion-by-regions method [45].

Next-to-next-to-leading order: complete contribution
The result obtained in this paper, together with the results given in [12], can be compared with the reference [33] which gives the full NNLO contribution to the Drell-Yan process. To this end, we recall that within the present approach the full partonic cross section is given according to Eq. (2.7), that is, as the sum of a dynamic and a kinematic contribution. Writing explicitly all the terms contributing to the cross section, we have 8 (4.10) The first three terms have been calculated in [12], and are reported explicitly in Eqs. (C.1), (C.2) and (C.3) of App. C. The last term is given by the sum of the S 1 , S 3 , S 4 and S 5 contributions given in (4.5), (4.7) and (4.9) respectively. Explicitly, it reads This is the non-singlet contribution to the qq unrenormalized partonic cross section, expressed in terms of the bare coupling constant. To compare with [33] we need to write it in terms of the UV-renormalized coupling constant and remove the initial state collinear singularities at cross-section level. These procedures amount to adding the following counterterms to the full NNLO cross section: when both sides are evaluated in terms of the renormalized coupling constant, for µ f = µ r = Q. In this equation the symbol ⊗ indicates convolution 9 : (4.14) furthermore, by ∆ (1) (z)| 0 and ∆ (1) (z)| we indicate respectively the coefficient of the 0 and terms of the NLO qq Drell-Yan correction, and P 0 qq and P 1,NS qq represent the one-and two-loop Altarelli-Parisi splitting functions, that we provide for completeness in App. D. We note that Eq. (4.13) is valid to all powers in (1 − z). Expanding all terms in (1 − z) and selecting the NLP contribution we finally get the finite partonic cross section which agrees with the NLP content of Eq. (B.7) in [33], for µ f = µ r = Q, namely ∆ qq,AC expanded to NLP. These terms are provided explicitly in [33] in Eqs. (B.30) -(B.33).

Next-to-next-to-next-to-leading order
Recently, partial results for the NLP expansion of the C 3 F contribution to the N 3 LO Drell-Yan cross section were calculated in [32] using the expansion-by-region method. We are able to compare to this result by combining the result of the collinear function J obtained in this paper. To this end, we now focus on the following part of the factorization formula expanded to the third order in the coupling constant We perform the convolution according to (4.16) and arrive at the following d-dimensional result Setting Ω = Q(1 − z) and µ = Q, and expanding in we find Our expanded result in the equation above agrees with Eq. (45) of [32] up to the finite constant terms which are not reported there and a factor of two accounting for the anticollinear contribution.

Cusp anomalous dimension
In addition to the checks performed at the cross-section level in the two sections above, we use the leading pole of the two-loop soft function S 1 to extract the first diagonal entry of the anomalous dimension matrix defined in Eq. (3.50) of [6] finding agreement. Currently, the resummation beyond LL is hampered by the appearance of endpoint divergent convolutions [12]. However, once cured, the results obtained in this work will be useful to extract the soft anomalous dimension matrix beyond LL accuracy.

Conclusions
In this article, we calculated the real-real contributions to the NLP generalized soft functions, which enter the bare factorization theorem for the Drell-Yan process in the threshold region [12]. This allowed us to complete the comparison of the NNLO Drell-Yan cross-section up to NLP against existing fixed-order results. The generalized soft functions, listed in Eqs. (2.10) -(2.14), contain a dependence on additional convolution variables ω or ω 1 , ω 2 with respect to the LP soft function. We carried out the calculation by employing methods developed for fixed-order calculations such as the reduction to MIs and the use of the differential equations for the direct evaluation of the MIs. Our results retain the exact d-dimensional dependence on the convolution variables at the integration boundaries which allows us to perform the convolution integrals with collinear functions at fixed-order accuracy. Given the current issues stemming from the expansion in d → 4 of the soft and collinear functions before the convolution is performed, we leave the non-trivial study of the renormalization procedure of the generalized soft functions for future work.
We showed that combining the soft functions with their respective collinear functions, as prescribed by the NLP factorization theorem, and performing the d-dimensional integrals yields the correct NNLO cross-section expressions up to NLP in the threshold expansion. In addition, we reproduced partial N 3 LO results available in the literature. We also confirmed the result for the diagonal entry of the anomalous dimension computed in [6], and explicitly validated to NNLO the assumption made in [6]
It is useful to write the J 1 (n + p, x a n + p A ; ω) collinear function in terms of two scalar components in the following way J 1;γβ (n + p, x a n + p A ; ω) = δ γβ J 1,1 (x a n + p A ; ω) δ(n + p − x a n + p A ) + J 1,2 (x a n + p A ; ω) ∂ ∂(n + p) δ(n + p − x a n + p A ) .
1,2 (n + p; ω) = 2 . and We also require the LP amplitudes +g 2 s n η 1 + n η 2 + 1 n + k 2 1 n + (k 1 + k 2 ) T K 1 T K 2 + 1 n + k 1 1 n + (k 1 + k 2 ) T K 2 T K 1 η 1 (k 1 ) η 2 (k 2 ) +g 2 s −n η 1 − n η 2 (B.6) For S 3 we need C Terms contributing to the NNLO cross section We report in this appendix the terms contributing to (4.10), which have been calculated in [12]. First of all one has the contribution to the partonic cross section due to the kinematic correction at NNLO: