Nested soft-collinear subtractions in NNLO QCD computations

We discuss a modification of the next-to-next-to-leading order (NNLO) subtraction scheme based on the residue-improved sector decomposition that reduces the number of double-real emission sectors from five to four. In particular, a sector where energies and angles of unresolved particles vanish in a correlated fashion is redundant and can be discarded. This simple observation allows us to formulate a transparent iterative subtraction procedure for double-real emission contributions, to demonstrate the cancellation of soft and collinear singularities in an explicit and (almost) process-independent way and to write the result of a NNLO calculation in terms of quantities that can be computed in four space-time dimensions. We illustrate this procedure explicitly in the simple case of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal O(\alpha _\mathrm{s}^2)$$\end{document}O(αs2) gluonic corrections to the Drell–Yan process of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q \bar{q}$$\end{document}qq¯ annihilation into a lepton pair. We show that this framework leads to fast and numerically stable computation of QCD corrections.

However, in spite of these remarkable successes, it is important to recognize that existing implementations of subtraction schemes are complex, not transparent, and require significant CPU time to produce stable results. On the other hand, slicing methods, while conceptually simple, have to be carefully controlled to avoid dependence of the final result on the slicing parameter. Given these shortcomings, it is important to study whether improvements to existing methods are possible. In the context of the N -jettiness slicing method, there has been recent progress toward a better control of the soft-collinear region [65,66].
In this paper we study the residue-improved subtraction scheme introduced in Refs. [10][11][12]. This scheme is interesting because it is the only existing framework for NNLO QCD computations that is fully local in multi-particle phase space. As such, it should demonstrate exemplary numerical stability, at least in theory. Although this scheme is well understood and was applied to a large number of non-trivial problems, we will argue in this paper that certain aspects of it are redundant. Interestingly, once this redundancy is recognized and removed, the residue-improved subtraction scheme becomes very transparent and physical. In addition, the technical simplifications that occur become so significant that the cancellation of the divergent terms can be demonstrated independently of the hard matrix element and almost entirely analytically, and the final finite result for the NNLO contribution to (in principle) any process can be written in a compact form in terms of generic four-dimensional matrix elements. 2 Although the improvements that we just described hold true for an arbitrary complicated process, in this paper, for the sake of clarity, we restrict our discussion to the production of a colorless final state in qq annihilation. This allows us to discuss all the relevant conceptual and technical aspects of the computational framework, without cluttering the notation and limiting the bookkeeping to a minimum. The generalization of the framework described here to arbitrary processes is -at least conceptually -straightforward.
Admittedly, compared to NNLO QCD problems studied recently, the production of a colorless final state in qq annihilation is a very simple process, which has been discussed in the literature many times. However, we believe that the simplicity of our approach and the structures that emerge justify revisiting it one more time. Moreover, thanks to the simplicity of this process, we will be able to describe our approach in detail and demonstrate many intermediate steps of the calculation. Hopefully this will allow us to make the rather technical subject of NNLO subtractions accessible to a broader part of the particle physics community.
The paper is organized as follows. We begin with preliminary remarks in Sect. 2, where we also precisely define the problem that we plan to address. In Sect. 3 we discuss the next-to-leading order (NLO) QCD computation as a prototype of the following NNLO QCD construction. In Sect. 4, we describe how the NLO computation generalizes to the NNLO case. We elaborate on this in Sect. 5, where we discuss ultraviolet and collinear renormalization, and in Sects. 6, 7 and 8, where we study two-loop virtual corrections, one-loop corrections to single-real emission process, and the double-real emission contributions, respectively. In Sect. 9, we combine the different contributions and present the final result for the NNLO QCD corrections to color singlet production in qq annihilation. In Sect. 10, we show some numerical results and a comparison with earlier analytic calculations. We conclude in Sect. 11. A collection of useful formulas is provided in the appendices.

Preliminary remarks
We consider the production of a colorless final state V in the collision of two protons P + P → V + X. (2.1) We are interested in computing the differential cross section for the process in Eq. (2.1) where dσ i j is the finite partonic scattering cross section, f i, j are parton distribution functions and x 1,2 are momenta fractions of the incoming hadrons that are carried to a hard collision by partons i and j, respectively. The dependence on the renormalization and factorization scales and all other parameters of the process is understood. The finite partonic scattering cross section is obtained after the renormalization of the strong coupling constant removes all ultraviolet divergences, all soft and final state collinear divergences cancel in the sum of cross sections with different partonic multiplicities, and the initial state collinear divergences are subtracted by redefining parton distribution functions. Since the process under consideration is driven by a conserved current that is independent of α s , the ultraviolet renormalization reduces to the following (MS) redefinition of the strong coupling constant: where S = (4π) e − γ E , γ E ≈ 0.577216 is the Euler-Mascheroni constant, = (4 − d)/2 and is the leading-order (LO) QCD β-function.
Collinear divergences associated with initial state QCD radiation are removed by a redefinition of parton distributions. In the MS scheme, this amounts to the replacement (2.5) where ⊗ stands for the convolution andP (0, 1) i j are the Altarelli-Parisi splitting functions. As we already mentioned, we focus on gluonic corrections to the qq annihilation channel, q +q → V + ng. (2.7) This allows us to present all the features of the framework while limiting the bookkeeping to a minimum and, therefore, to keep the discussion relatively concise. All other partonic channels relevant for the Drell-Yan process can be obtained by a simple generalization of what we will describe. The collinear-renormalized partonic cross section for qq annihilation into a vector boson is expanded in series of α s . We write dσ ≡ dσ qq = dσ LO + dσ NLO + dσ NNLO , (2.8) where dσ NNLO = dσ VV + dσ RV + dσ RR + dσ ren + dσ CV . (2.9) Various contributions in Eq. (2.9) refer to virtual and real corrections, as well as to contributions to cross sections that arise because of the ultraviolet and collinear renormalizations. The latter are obtained with the procedure just described and read dσ ren = − α s (μ) 2π β 0 dσ NLO , 11) and the relevant splitting functions are provided in Appendix A. The cross section dσ (N)NLO is finite but all the individual contributions in Eq. (2.9) are divergent. The well-known problem is that these divergences are explicit in some of the terms and implicit in the others. Indeed, soft and collinear divergences appear as explicit 1/ poles in virtual corrections but they only become evident in real corrections once integration over gluon momenta is performed. However, since we would like to keep the kinematics of all the final state particles intact, we cannot integrate over momenta of any of the final state particles if it is resolved. It is this point that makes extraction of implicit singularities complicated and requires us to devise a procedure to do it.
Depending on how these implicit singularities are extracted, it may or may not be straightforward to recognize how they combine and cancel, once all contributions to the physical cross section are put together. At NNLO, this was done for the antenna subtraction scheme and, in a less transparent way, for the residue-improved sector decomposition. One thing we would like to do, therefore, is to combine the individual terms that contribute to partonic cross sections, and cancel all the 1/ divergences explicitly, without any reference to the matrix elements that contribute to the different terms in Eq. (2.9). In the next section we show how to do that at next-to-leading order in the perturbative expansion for the Drell-Yan process. This will allow us to set up the formalism and the notation that will be used for the NNLO analysis of Sects. 4-8.

The NLO calculation
We will illustrate our approach by studying the production of a lepton pair in quark-antiquark annihilation at next-toleading order in perturbative QCD. We note that, at this order, the method that we would like to describe is identical to the FKS subtraction scheme introduced in Refs. [67,68]. However, we formulate the FKS method in a way that makes its extension to next-to-next-to-leading order as straightforward as possible. 3 One point that we found helpful, especially for bookkeeping, was to introduce soft and collinear subtraction operators, and we show how to use them in the NLO computation below.
We are interested in the calculation of the finite partonic cross section dσ NLO defined in Eq. (2.9). It receives contributions from the virtual and real corrections and the collinear subtraction term. We will start the discussion with the realemission contribution. It refers to the process where V is a generic notation for all colorless particles in the final state. We write the cross section for the process in Eq. (3.1) as where s is the partonic center-of-mass energy, and In Eq. (3.4), dLips V is the Lorentz-invariant phase space for colorless particles, including the momentum-conserving is the matrix element for the process in Eq. (3.1) and F kin (1, 2, 4, V ) is an (infrared safe) observable that depends on kinematic variables of all particles in the process. Also, E max is an arbitrary auxiliary parameter that has to be large enough to accommodate all possible kinematic configurations for qq → V + g. The need to introduce such a parameter is a consequence of our construction, as explained in detail below. We would like to isolate and extract soft and collinear singularities that appear when the integration over [dg 4 ] in Eq. (3.2) is attempted. To this end, we introduce two operators that define soft and collinear projections 3 For earlier efforts, see Ref. [69].
where ρ i j = 1 − n i · n j and n i is a unit vector that describes the direction of the momentum of the ith particle in (d − 1)dimensional space. By definition, operators in Eq. (3.5) act on everything that appears to the right of them. The limit operations, on the right hand side of Eq. (3.5), are to be understood in the sense of extracting the most singular contribution provided that limits in the conventional sense do not exist. We will also use the averaging sign ... to represent integration over momenta of final state particles. This integration is supposed to be performed in the center-of-mass frame of incoming partons. We emphasize that this remark is important since our construction of the subtraction terms is frame-dependent and not Lorentz-invariant. We rewrite Eq. (3.2) in the following way: where I is the identity operator andÔ NLO is a short-hand notation for a combination of soft and collinear projection operatorŝ Note that in Eq. (3.6) soft and collinear projection operators act on F L M (1,2,4), which, according to Eq. (3.4), contains the energy-momentum conserving δ-function; we stress that the soft and collinear limits must be taken in that δ-function as well.
The reason for rewriting dσ R as in Eq. (3.6) is that the last term there is finite, thanks to the nested structure of subtraction terms. This term can, therefore, be integrated numerically in four dimensions. We emphasize again that the subtraction terms, as formulated here, are not Lorentz-invariant. This means that all the three terms in Eq. (3.6) should be computed in the same reference frame that, as already mentioned, is taken to be the center-of-mass reference frame of the colliding partons.
We now consider the remaining two terms in Eq. (3.6). Their common feature is either complete or partial decoupling of the gluon g 4 from the matrix element thanks to the fact that they contain either soft or collinear projection operators. Hence, those terms can be re-written in such a way that all singularities are extracted and canceled, without specifying the matrix elements for the hard process.
To see this explicitly, consider first two terms in Eq. (3.6) and write them as (1,2,4) . (3.8) It is easy to see that the first term in Eq. (3.8) vanishes. 4 Indeed, in the limit when the gluon g 4 becomes soft, we find where g s,b is the bare QCD coupling, C F = 4/3 is the QCD color factor, and F L M (1,2) is closely related to the LO cross section (3.10) The action of the collinear operators on the ρ gives Since for head-on collision ρ 12 = 2, ρ 24 = 2 − ρ 14 , we find Hence, the only term that we need to consider in Eq. (3.8) is the collinear subtraction (1,2,4) . (3.13) We will consider the action of the operator C 41 on F L M (1,2,4) and then infer the result for the operator C 42 . First, we find the collinear limit (3.14) We note that a new variable z = 1 − E 4 /E 1 is introduced in Eq. (3.14). The notation F L M (z · 1, 2) implies that in the computation of F L M (1, 2), cf. Eq. (3.10), the momentum p 1 is replaced with zp 1 everywhere, including the energymomentum conserving δ-function. We also used P qq (z) to denote the splitting function To simplify C 41 F L M (1,2,4) , we integrate over the emission angle of the gluon g 4 , rewrite the integration over its energy as an integral over z and express g s,b in terms of the renormalized coupling α s (μ). After straightforward manipulations we find where z min = 1− E max /E 1 and we introduced the short-hand notation We note that in Eq. (3.16) integration over z leads to divergences caused by the soft z → 1 singularity in the splitting functions. These singularities are regulated dimensionally in Eq. (3.16). On the other hand, this equation has the form of a convolution of a hard matrix element with a splitting function, so we expect that divergences present there will cancel against the collinear subtraction terms. However, collinear subtractions employ regularization of soft singularities that is based on the plus-prescription. Our goal, therefore, is to rewrite Eq. (3.16) in such a way that all soft singularities are regulated by the plus-prescription; once this is done, combining this contribution with virtual corrections and collinear subtractions becomes straightforward.
To simplify the notation, we denote F L M (z · 1, 2)/z = G(z) and split P qq (z) into a piece that is singular at z = 1 and a regular piece We also note that we can extend the integral over z in Eq. (3.16) to z = 0 since if E 4 > E max , F L M (z · 1, 2) will vanish because there is not enough energy to produce the final state. We will use this fact frequently in our NNLO analysis. We write . (3.19) 123 The expression in Eq. (3.19) can be expanded in a power series in to the required order and the plus-distributions can be used to write the result in a compact form. Indeed, the following equation holds: It is now straightforward to rewrite Eq. (3.16) in such a way that all soft, z → 1, singularities are regulated using the plus-prescription. We use the fact that we are in the center-of-mass frame of the incoming qq pair, so that 2E 1 = 2E 2 = √ s. We find The splitting function in Eq. (3.21) reads 5 (3.23) The result for C 42 F L M (1,2,4) is obtained by a simple replacement F L M (z · 1, 2) → F L M (1, z · 2) in Eq. (3.21).
Putting everything together, we find the following result for the real-emission cross section: (3.24) 5 The O( 2 ) contribution to P qq,R , relevant for NNLO contributions, is reported in Appendix A.
We note that in Eq. (3.24) all singularities of the real-emission contribution are explicit and a straightforward expansion in is, in principle, possible. However, such an expansion is inconvenient since it involves higher-order terms of lowermultiplicity amplitude. To avoid these contributions, it is useful to combine Eq. (3.24) with virtual corrections and collinear counterterms.
For the virtual corrections, all divergent parts can be separated using Catani's representation of renormalized one-loop scattering amplitudes [70]. We obtain where F fin LV (1, 2) is free of singularities and μ-independent. To arrive at the final result, we add virtual, real and collinear subtraction terms, cf. Eq. (2.9), and obtain 27) and the extra z terms in the denominator of the last line of Eq. (3.26) appear because of the z-dependent flux factor in the collinear counterterm cross section. Taking the limit → 0 in Eq. (3.26), we find the final formula for the NLO QCD contribution to the scattering cross section for q( p 1 ) +q( p 2 ) → V + X . It reads It follows that the NLO cross section is computed as a sum of low-multiplicity terms, including those where F L M (z · 1, 2) or F L M (1, z · 2) are convoluted with particular splitting functions, and the subtracted real-emission term described by O NLO F L M (1,2,4) . We note that terms that involve matrix elements of different multiplicities, as well as terms that involve different types of convolutions, are separately finite. We will use this observation at NNLO, to check for the cancellation of 1/ divergences in an efficient way.

The NNLO computation: general considerations
We would like to extend the above framework to NNLO in QCD. Apart from the UV and collinear renormalization discussed in Sect. 2, the NNLO cross section receives contributions from two-loop virtual corrections to qq → V (double virtual), from one-loop corrections to the process with an additional gluon in the final state qq → V + g (real-virtual), and from the tree-level process with two additional gluons in the final state qq → V + gg (double real). The double-virtual corrections can be dealt with in a straightforward way since all the divergences of the two-loop matrix elements are explicit, universal and well understood [70]. For our purposes, we only need to write them in a convenient form. The real-virtual corrections are more tricky, but do not require new conceptual developments. Indeed, the kinematic regions that lead to singularities in one-loop amplitudes with an additional gluon in the final state are identical to those appearing in the NLO computations and, furthermore, the limiting behavior of one-loop amplitudes with one additional parton is well understood [71][72][73]. Hence, we can deal with the real-virtual contribution by a simple generalization of what we did at NLO.
The qualitatively new element of the NNLO computation is the double-real emission process qq → V + gg. The methods that are applicable at next-to-leading order need to be adjusted to become useful in the NNLO case. However, somewhat surprisingly, these adjustments appear to be relatively minor although, of course, the bookkeeping becomes much more complex.
We begin by discussing the ultraviolet and PDF renormalizations at NNLO, as well as the double-virtual and the real-virtual contributions. We then move on to a more involved analysis of the double-real emission contribution to dσ NNLO .

The NNLO computation: ultraviolet and PDF renormalization
In this section we study the contributions to dσ NNLO coming from the ultraviolet and the collinear renormalization, beginning with the former. As mentioned previously, because the process qq → V is driven by a conserved current which is independent of α s , the ultraviolet renormalization contribution is very simple. Combining Eq. (2.10) and the first two lines of Eq. (3.26), it is straightforward to obtain We proceed to the collinear subtraction. Rewriting Eq. (2.10) to make the convolutions explicit, we obtain Terms that involve convolutions of the various splitting functions with F L M are, in principle, straightforward to deal with. These terms are fully regulated and can be expanded in powers of without further ado. In practice, we combine those terms with other contributions in order to cancel the singularities prior to integration over z.

123
It is less straightforward to rewriteP ⊗ dσ R+V and dσ R+V ⊗P in a form convenient for combining them with other contributions to dσ NNLO . We focus onP ⊗dσ R+V , and consider the effect of the convolution on the first two lines of Eq. (3.26).
First, we consider the term proportional to S(1, 2) F L M (1,2) in Eq. (3.26). It receives contributions from the divergent part of virtual corrections and from the soft regularization of collinear subtraction terms. These terms scale differently with z. The virtual correction depends on s − which becomes (zs) − once the momentum p 1 is changed to zp 1 . On the other hand, the soft remainders of the collinear subtracted terms scale as E −2 i , with i = 1, 2. Hence, in the calculation of dσ R+V (z · 1, 2), the corresponding contribution scales with z either as ∼ z −2 or as ∼ 1. Therefore, we have to compute For future convenience, we rewrite Eq. (5.3) as where the splitting function P qq,NLO CV (z) is defined in Eq. (A.6). The other two terms that we need involve convolutions of splitting function P qq,R with F L M (1, 2), cf. Eq. (3.26). The first term can be written as a double convolution where the splitting function P qq ⊗ P qq NLO CV is defined in Eq. (A.7). The second term is the left-right convolution, Combining all these terms we find Each term that appears on the right hand side of Eq. (5.8) is regularized and can be expanded in powers of independently of the other terms in that equation.

The NNLO computation: double-virtual corrections
The calculation of double-virtual corrections proceeds in the standard way. We start with the scattering amplitude for qq → V and write it as an expansion in the renormalized strong coupling constant where the dependence of the scattering amplitudes on the momenta of the external particles is suppressed. By analogy with what was done in Sect. 3, we write where in the second line we removed the renormalization contribution that is already accounted for in Eq. (5.1). F LV V can be directly expanded in a Laurent series in and integrated over the phase space of the final state particles since this integration does not introduce soft or collinear divergences. Before doing that, it is convenient to explicitly extract the 1/ poles from F LV V . Soft and collinear singularities of a generic scattering amplitude are given in Ref. [70]. Using these results, we rewrite F LV V as 3) The sum of the last two terms in Eq. (6.3) is a finite remainder of the O(α 2 s ) contribution to the virtual corrections once its divergent part is written in a form suggested in Ref. [70]. More specifically, F fin LV 2 is the finite remainder of the oneloop amplitude squared, while F fin LV V is the genuine two-loop finite remainder. Note that, contrary to F fin LV (1, 2) and F fin As follows from Eq. (6.3), the singularities of the doublevirtual corrections are proportional to the leading-order contribution F L M (1,2) and to the finite part of the virtual corrections to the NLO cross section F fin LV (1, 2), given in Sect. 3. Our goal is to rewrite the real-virtual and the double-real emission contributions in a way that allows explicit cancellation of the divergences in Eq. (6.3) without specifying hard matrix elements.

The NNLO computation: real-virtual corrections
The kinematics of the real-virtual corrections is identical to the NLO case described in Sect. 3. The procedure for making these corrections expandable in is, therefore, the same. We write We remind the reader that, according to our notation, softand collinear-projection operators in Eq. (7.1) do not act on the phase space of the gluon g 4 . It remains to compute the corresponding limits in Eq. (7.1) and to rewrite them, where appropriate, as convolutions of the hard matrix elements with splitting functions. The soft limit of a general one-loop amplitude is discussed in Refs. [71,72]. Adapting those results to our case, we find We need to integrate Eq. (7.2) over the phase space of the gluon g 4 . This can easily be done, with the result Note that in order to obtain a meaningful result, it is crucial that the integration over gluon energy is bounded from above; as we already explained in the context of the NLO computations, we use a parameter E max for this purpose, cf.
The second term that we need to consider is the softregulated collinear subtraction term We will only discuss the collinear projection operator C 41 ; the contribution of C 42 is obtained along the same lines. Collinear limits of loop amplitudes were studied in Refs. [71,73]. Using these results and adapting them to our case, we find where the splitting function P loop,i qq is given in Eq. (A.8). The soft-collinear limit is easily obtained by taking the collinear approximation in Eq. (7.2). We find Integrating over emission angles of the gluon g 4 and rewriting the result through plus-distributions, following the discussion in Sect. 3, we obtain a convenient representation for the softregulated collinear subtraction term. It reads where the two splitting functions are defined in Eqs. (A.10) and (A.11). We replace We are now in a position to present the final result for the real-virtual part. Collecting results shown in Eqs. (7.3) and (7.7), we obtain We stress that each term on the r.h.s. in Eq. (7.8) can be expanded in powers of ; we will make full use of this to cancel 1/ singularities when combining Eq. (7.8) with other contributions to dσ NNLO . To this end, it is useful to make all the 1/ singularities explicit in Eq. (7.8) by writing F LV (1,2) in terms of F fin LV (1,2) and We used s i j = 2 p i · p j and denoted a finite remainder which does not depend on the scale μ by F fin LV (1, 2, 4).

General considerations
In this section, we discuss the double-real emission contributions to dσ NNLO . Similar to the NLO case, we need to determine all kinematic configurations that may lead to singularities and understand the factorization of the matrix element that describes qq → V + gg in these regions. In the case of the two-gluon emission in qq annihilation into a colorless final state, the singular regions correspond to soft and/or collinear emissions, with collinear directions being the collision axis and the direction of either one of the two gluons. The difficulty of the NNLO case is that each of these kinematic limits can be approached in several different ways and all of them have to be identified and regularized separately. To do so, we introduce several soft and collinear projection operators. They are defined as follows. Consider a quantity A that depends on the four-momenta of some or all of the particles in the process. The action of operators SS, S 4,5 , CC 1,2 , C 14 , C 15 , C 24 , C 25 , C 45 on A is described by the following formulas: To make use of these projection operators, we need to rewrite the two-gluon phase space in a way that allows these limits to be taken. It is convenient to order gluon energies as the first step. We write where, as in Sect. 3, dLips is the phase space for the final state V , including momentum-conserving delta-function. The gluon phase space elements [dg 4,5 ] are defined as in Eq. (3.3) As we already saw when considering the real-virtual contribution, it is necessary to introduce the θ -function in order to define integrals over gluon energies in the soft limits. We work in the center-of-mass frame of the colliding quark and antiquark; it is in this frame that all the energies in the above formulas are defined. We recall that, similar to the NLO case, soft and collinear projection operators act on everything that appears to the right of them. However, in the NNLO case we will find it convenient, occasionally, to also simplify the phase space in the collinear limits. If so, we will explicitly show the corresponding part of the phase space to the right of the operator. For example We begin by extracting soft singularities from the doublereal process, largely repeating what we have done at next-toleading order. 6 We write In Eq. (8.5), the last term is soft-regulated, in the second term gluon g 5 is soft and soft singularities associated with g 4 are regulated, and in the first term both g 4 and g 5 are soft. All of these terms contain collinear singularities. Regulating these is more difficult since collinear singularities overlap. For this reason, we first need to split the phase space into mutually exclusive partitions that, ideally, select a single kinematic configuration that leads to singularities. We write where the label i runs through the elements of the set i ∈ {14, 15; 24, 25; 14, 25; 15, 24}. We refer to the first two elements of the set as triple-collinear and to the last two elements of the set as double-collinear partitions. We construct the weights w i in such a way that when they multiply the matrix element M(1, 2, 4, 5) squared, the resulting expression is only singular in a well-defined subset of limits. For example, in the partition 14, 15 collinear singularities in w 14,15 |M| 2 only occur when gluons 4 and/or 5 are emitted along the direction of the incoming quark q( p 1 ) or when their momenta are parallel to each other. Similarly, in the partition 24, 25, the singularities occur when momenta of g 4 and/or g 5 are parallel to p 2 or to each other. In the partitions 14, 25 and 15, 24, singularities only occur when momenta of g 4 and g 5 are collinear to p 1 and p 2 or p 2 and p 1 , respectively. Apart from these requirements, the specific form of w i is arbitrary. Weights used in this calculation are given in Appendix B. In the following, we assume that weights do not depend on gluon energies and, therefore, commute with soft operators. The triple-collinear partitions require further splitting to factorize all the relevant singularities. The purpose of this splitting is to establish a well-defined hierarchy for the parameters ρ 4i , ρ 5i , ρ 45 , since different orderings correspond to different limiting behavior. This splitting is not unique; a possible choice consistent with the phase space parametrization that we employ later (cf. Appendix B) reads where, as usual, We will refer to the four contributions shown in Eq. (8.7) as sectors a, b, c and d. We note that only two of the sectors are qualitatively different, since the other two are just obtained by the 4 ↔ 5 replacement. However, because of the energy ordering E 5 < E 4 , we no longer have the 4 ↔ 5 symmetry, and we have to consider all the four sectors separately.
A suitable parametrization of all angular variables that supports splitting of the angular phase space as shown in Eq. (8.7) and allows for factorization of singularities in hard amplitudes was provided in Refs. [10,11] and is reviewed in Appendix B. We will use this parametrization to carry out integrations over sectors θ (a) Having introduced partitions and sectors as a tool to identify singularities that may appear in the course of integrating over the angles of the final state gluons, we are now in a position to write the result for the double-real emission cross section as a sum of terms that either can be integrated in four dimensions, or that depend on hard matrix elements of lower multiplicity. The latter contributions still diverge, either explicitly or implicitly, and we will have to combine them with double-virtual and real-virtual contributions to arrive at the finite result.
We can thus rewrite the double-real emission cross section as where the soft-regulated, single-collinear term F s r c s L M reads (1,2,4,5) ,  For example, the contribution of the double-collinear sector i4, j5 follows from an expansion of an identity operator written in the following form: The reason we restrict ourselves to the subtraction of the C 4i and C 5 j collinear projection operators is that in the partition i4, j5 no other collinear singularities appear, thanks to the damping factor w i4, j5 . Similarly, taking e.g. the sector a of the triple-collinear partition w i4,i5 , we write (8.13) because in this case a singularity can only occur either in a triple-collinear limit η 4i ∼ η 5i → 0 or if η 5i → 0 at fixed η 4i . It is worth pointing out a few things in connection with Eq. (8.8).
• The procedure that we used to write Eq. (8.8) is, in principle, process-and phase space parametrizationindependent. We will use a particular parametrization of the phase space to perform the required computation but one should keep in mind that the freedom of changing the parametrization exists and, perhaps, it is worth exploring it in the future.  (1,2,4) and matrix elements of lower multiplicity. This term still contains unregulated singularities that occur when the momentum of the gluon g 4 becomes collinear to the collision axis or to the direction of g 5 . When integrated over gluon energies and angles, this term gives rise to 1/ n poles, with n ≤ 3. • The term F s r c s L M in Eq. (8.9) is the soft-regulated singlecollinear subtraction term. Note that thanks to the damping and θ factors only one kind of collinear singularity per term is present. F s r c s L M involves F L M (1, 2, 4(5)), depending on the partition and the sector and, therefore, contains unregulated collinear singularities related to gluon emissions along the collision axis. It gives rise to 1/ 2 and 1/ poles. • The term F s r c t L M in Eq. (8.10) is the triple-collinear subtraction, where all other singularities are regulated. As we will see, contributions of the double-collinear partitions to F s r c t L M have the "double-convolution" structure. The F s r c t L M term contains 1/ poles in contributions of triple-collinear partitions, and 1/ 2 poles in contributions of double-collinear partitions.
• The term F s r c r L M in Eq. (8.11) is completely regulated, thanks to the nested subtractions. It can be evaluated in four dimensions. It is the only term that involves the full hard matrix element for the process q( p 1 ) +q( p 2 ) → V + g( p 4 ) + g( p 5 ).
Following our general strategy, we need to study the first four terms in Eq. (8.8), which involve matrix elements of reduced multiplicity, and rewrite them in terms of integrable quantities that admit straightforward expansions in the dimensional regularization parameter . We will discuss how to do this in the following subsections.

The double-soft subtraction term
We begin with the discussion of the first term in Eq. (8.8), (1,2,4,5) . It corresponds to the kinematic situation where momenta of both gluons vanish at a comparable rate. The corresponding limit for the amplitude squared is given in Refs. [74,75] and allows us to write where and [75] At this point, we stress again that the hard matrix element F L M (1, 2) corresponds to a tree-level process and that the emitted gluons have no impact on the kinematic properties of the final state V because they decouple from the energymomentum conserving δ-function.
The goal now is to integrate the eikonal factor over the momenta of the two gluons. We note that, at this point, unless put in by hand, the integration over gluon energies becomes unconstrained since the energy-momentum conserving δfunction becomes independent of the gluon momenta after the double-soft limit is taken. It is for this reason that we need to introduce E max as in Eq. (8.3).
To satisfy constraints on gluon energies, we parametrize them as Written in these variables, the eikonal factor becomes The important point is that the dependence on the overall energy scale x 1 factorizes and that the remaining (complicated) function E is independent of energies of the incoming partons. We also use the parametrization of energies Eq.  For the case of a color-singlet final state, this integral is just a constant. 7 The abelian contribution is simple to obtain since it is just the product of NLO structures. In principle it should be possible to compute the non-abelian contribution analytically along the lines of e.g. Refs. [76][77][78]. However, it is also straightforward to obtain it numerically. To do this, we partition the phase space as in residue-improved sector decomposition [10,11]. The corresponding formulas for the angular phase space are given in Appendix B. Performing the required decomposition and integrating Eq. (8.20) numerically, we obtain the -expansion of the double-soft subtraction term, 7 In a more general NNLO problem, this integral is a function of the scalar product of the three-momenta of the two hard partons. (1, 2, 4

SS F L M
Numerical values of the coefficients c SS are shown in Table 1. There we also report numerical results for the abelian contribution, which are in perfect agreement with the analytic calculation. The result for the double-soft subtraction SS F L M (1, 2, 4, 5) does not require any further regularization; we will later combine it with other contributions with tree-level kinematics to cancel the 1/ singularities explicitly.

The single-soft term
We now consider the second term that contributes to Eq. (8.8). It is a double-soft regulated, single-soft singular expression that reads (1,2,4,5) . (8.22) Note that this contribution implicitly depends on F L M (1,2,4) and F L M (1,2), and the hard matrix element that appears in F L M (1, 2, 4) still contains collinear singularities that arise when the momentum of gluon g 4 becomes parallel to the momenta of the incoming partons or to the direction of g 5 . These divergences will have to be extracted and regulated. We start by computing the soft limit for the gluon g 5 . We find (1, 2, 4). (8.23) Since the gluon g 5 decouples from the hard matrix element, we can integrate over its momentum. We find (1, 2, 4) , and We need to simplify Eq. (8.24) because it still contains collinear singularities that appear when the momentum of gluon g 4 becomes parallel to the collision axis. To extract them, we write We reiterate that according to our notational conventions, the collinear projection operators do not act on the phase space element of the gluon g 4 in Eq. (8.27). The first term in Eq. (8.27) is explicitly regulated and can be expanded in powers of ; for this reason, we will only be concerned with the second term. We focus on the projection operator C 41 ; the contribution of the projection operator C 42 is then obtained by analogy. First, we consider how C 41 acts on J 124 . Using η 12 = 1, C 41 ρ 24 = ρ 12 and taking the ρ 41 → 0 limit on K 14 , K 24 we obtain Integrating over the energy and angle of the gluon g 4 we arrive at where z min = 1 − E max /E 4 . The splitting function operator P (−) qq is defined by means of the following equation: Putting everything together, we find the final result for the double-soft regulated single-soft singular contribution to F L M (1,2,4,5) , We need to rewrite Eq. (8.34) in such a way that extraction of the remaining collinear singularities becomes straightforward. We note that F s r c s L M contains contributions from double-and triple-collinear partitions, which we will treat separately. We will start with the double-collinear partitions since they are somewhat simpler.

The double-collinear partitions
In this subsection, we will consider the contribution of the double-collinear partitions to F s r c s L M . We begin with the partition 14, 25. For the first term, we need to compute The function P qq (z) was introduced in Eq. (8.31). We now consider the phase space. According to Eq. (8.34), C 41 acts on the phase space element [dg 4 ]. We introduce In this parametrization, the action of C 41 implies replacing Putting everything together, we obtain Note that, similar to the NLO case, the lower boundary z min is not important when integrating F L M (z · 1, 2, 4) since for z < z min there is no sufficient energy in the incoming partons to produce the required final state. Next we consider the action of the C 52 projection operator. Following the preceding discussion and using z = The sum of Eqs.
We can now combine the contributions of the two doublecollinear partitions. In doing so, it is convenient to always denote the "resolved" (i.e. the non-collinear) gluon by g 4 . Out of the four terms that we need to combine, two correspond to the collinear emission along the direction of the incoming quark p 1 and two along the direction of the incoming antiquark p 2 . We consider terms that belong to the former category first. When combining results, it is important to realize that Note that the lower integration boundary in this formula should be z = z min but we can extend the integration region to z = 0, without making a mistake. This is so because every time F L M (z · 1, . . .) appears in the integrand, the z > z min condition is automatically enforced by the requirement that the initial state should have enough energy to produce the final state. On the other hand, if F L M (z · 1, . . .) does not appear, θ -functions require that z > z 4 > z min . We also note that, thanks to explicit subtractions and constraints due to θ -functions, each term in Eq.  z · 2, 4) .
Note that the second term in the curly bracket only depends on z through the θ -function and so the z-integration of this term can be performed explicitly.

The triple-collinear partition 14, 15
We In parallel to the case of the double-collinear partitions, it is again convenient to always call the resolved gluon g 4 . We combine contributions of sectors (a) and (c), renaming g 5 → g 4 where appropriate, and we obtain (1, 2, 4) .
We now turn to sectors (b) and (d). These sectors are different from the other triple-collinear sectors and from the double-collinear partitions. Indeed, the single-collinear limits that we consider in sectors (b) and (d) correspond to gluons g 4 and g 5 becoming collinear to each other. We consider (8.52) and start with the discussion of how the collinear projection operator C 45 acts on F L M . We find where p 4+5 = p 45 = (E 4 + E 5 )/E 4 · p 4 , i.e. the hard matrix element must be taken in the collinear limit. The splitting functions are and z is the fraction of the total momentum p 45 = p 4 + p 5 carried by gluon g 5 , The vector κ ⊥ is a normalized transverse vector (8.56) defined by the following decomposition: 4 , − p 4 ) and k ⊥ · p 4 = k ⊥ ·p 4 = 0. We now construct these vectors explicitly. For this, we need the parametrization of the angular phase space of the two gluons valid for sectors (b) and (d); it is given in Appendix B.
Here we repeat the relevant formulas and discuss simplifications that occur in the limit where the momenta of g 4 and g 5 become collinear. We write the four-momenta of g 4 and g 5 as where t μ = (1, 0), e μ 3 = (0, 0, 0, 1; 0 . . .), b · t = b · e 3 = 0 and a · t = a · e 3 = a · b = 0. Our goal is to parametrize the phase space in such a way that explicit averaging over directions of k ⊥ can be performed. The phase space parametrization for sectors (b) and (d) can be written as Here, x 4 → 0 corresponds to the 4||5 collinear limit, η 45 = (1 − cos θ 45 )/2, x 3 = ρ 41 /2 and λ is related to the azimuthal angle ϕ 45 . Further details about the parametrization as well as expressions of scalar products in terms of x 3,4 and λ can be found in Appendix B. In this parametrization, the vector κ ⊥ reads 8 Using this expression in Eq. (8.53) together with momenta parametrization Eq. (8.58) and the phase space limit Eq. (8.60), we observe that integrations over λ and the directions of the vector a μ can be performed since neither λ nor a μ appear in the hard matrix element. We define and It follows that averaging over transverse directions introduces an -dependent leftover, as a consequence of the chosen phase space parametrization [10,11].
To write the result of the integration over unresolved phase space variables, it is convenient to define an additional splitting function, Combining the results discussed above, we write an expression for the contribution of the C 45 collinear projection operator in sector (b). We obtain (1,2,4,5), (8.68) where we only need to consider SSr μ r ν F μν L M (1, 2, 4 + 5). We find it using the known soft limits for amplitudes and the explicit form of the vector r given in Eq. (8.62). Indeed, since r · p 4 = 0 and r 2 = −1, r μ is a valid polarization vector of the gluon with momentum 4 + 5, in the collinear 4||5 approximation. For this reason, the soft limit of r μ r ν F μν L M follows from the standard soft limit of the amplitude for qq → V + g, not averaged over gluon polarizations. We obtain (1, 2, 4). (8.71) Collecting all the soft limits, we find (1, 2, 4), (8.72) 123 Fig. 1 Integration region for the (E 4 , E 5 ) → (E 45 , z) change of variables. The colored triangle is the allowed 0 < E 5 < E 4 < E max region. The blue region "A" is the "physical" one, i.e. the one which is not removed by a phase space θ-function inside F L M (1,2,45 where z is defined in Eq. (8.55). This implies We can now use Eqs. (8.72) in (8.68) and integrate over all variables that are not present in the hard matrix elements. This requires different variable transformations in the first and the second terms in Eq. (8.72). To integrate the first term, we change the integration variables from E 4,5 to E 45 and z. We find that the integration region splits into two regions (cf. Fig. 1 (1, 2, 4) .  (1, 2), (8.77) where z 4 = 1 − E 4 /E max and we expressed [dg 4 ] as an integral over energy E 4 and the angular integration variable x 3 = ρ 14 /2. We have also integrated over d − 3 angular variables that do not appear in the hard matrix element and in the splitting function. Finally, we consider sector (d). We need to compute The calculation is similar to what we just described for sector (b), apart from the following modifications of the integration boundaries: . (8.79) Incorporating these changes in Eqs. (8.76) and (8.77) provides us with the result for sector (d).
Significant simplifications occur if the results for the two sectors are added; this happens because the z-integration boundaries in sectors (b) and (d) complement each other. Also, for both I A and I B the z-integration decouples from the rest and can be performed independently of the hard matrix element. In I A , it yields (1, 2, 4) , where we used x 3 = ρ 14 /2 and the constantsγ g ( ),γ g ( , k ⊥ ) are defined in Eq. (A.19). In the integral I B the hard matrix element is that of the leading-order process which implies that integration over all variables related to radiated gluons can be explicitly performed. We find with δ g defined in Eq. (A.20).

Summing double-and triple-collinear partitions
Summing up the above results, we obtain an intermediate representation of F s r c s L M . We write it as a sum of four terms, (1, 2, 4) .  (1, 2, 4) , (8.83) (1, 2, 4) , The four contributions to F s r c s L M described at the end of the previous section require further manipulations because they cannot be expanded in series of as they are. Indeed, all of them exhibit collinear singularities in the limits 4||1 and 4||2, which need to be extracted before expansion in becomes possible. To deal with this issue, we again rewrite the identity operator through collinear projections. For example, we write C 1 (z · 1, 2, 4) = (C 41 + C 42 )C 1 (z · 1, 2, 4) 1, 2, 4) . (8.88) 123 The first two terms can be further simplified by considering respective collinear limits; the last term is regulated and can be Taylor-expanded in . The single-collinear subtraction term can be analyzed in the same way as all the other collinear limits discussed previously. The only new element here is the action of the collinear projection operators on the spin-correlated part. Using the explicit expression for the vector r in Eq. (8.62) we find where z is defined in the usual way z = 1 − E 4 /E 1 . Taking this into account, after tedious but straightforward calculations we arrive at (1, 2) .
In the above equation, we used the following notation: The relevant splitting functions are defined in Appendix A. Also, We will use these results when presenting the final formula for the double-real contribution.
8.6 The double-unresolved collinear limit: double collinear We now turn to the term F s r c t L M and begin by considering the contribution of the double-collinear partitions. It reads (1,2,4,5) . (8.95) Note that, following our notational convention, the collinear projection operators act on the phase space elements [dg 4 ] and [dg 5 ].
We begin with the C 41 C 52 term. Introducing and calculating collinear limits we obtain Since the momenta of gluons g 4 and g 5 decouple from each other, we find where, as usual, the z integrals do not need a lower cut-off whenever z is present in F L M . The term with collinear operators C 51 C 42 in Eq. (8.95) can be simplified in a similar way. Combining the two contributions, we obtain We can rewrite Eq. (8.101) to ensure that all singularities that appear in z andz integrals are regulated with the plus-prescription. This gives the final result for the doublecollinear contribution The relevant splitting functions are given in Appendix A.
8.7 The double-unresolved collinear limit: triple collinear In this section, we consider the contribution of the triplecollinear partitions to F s r c t L M . It reads (1, 2, 4, 5) .
This contribution always contains the triple-collinear projection operator CC i that acts on the hard matrix elements. For i = 1, this gives, schematically, where s 145 = ( p 1 − p 4 − p 5 ) 2 and P ggq (1,4,5) is the known triple-collinear splitting function [75,79,80]. The reduced matrix element in Eq. (8.104) has to be evaluated in the exact collinear limit, i.e.
Other projection operators that appear in Eq. (8.103) provide subtractions that are needed to make the triple-collinear splitting function integrable over the unresolved parts of the (g 4 , g 5 ) phase space. For definiteness, we focus here on the triple-collinear partition where gluons are emitted along the direction of the incoming quark with momentum p 1 ; this corresponds to taking i = 1 in Eq. (8.103).
To proceed, we note that the damping factors in Eq. (8.103) can be removed since the collinear projection operator CC 1 acting on them yields 1. Next, we need to study the triplecollinear limit of the angular phase space. The generic phase space parametrization is described in Appendix B and we use it to compute the triple-collinear limits. We stress that since the phase space parametrization changes from sector to sector, we need to consider all the four sectors separately.
Without going into further detail of the angular integration, it is clear that once this integration is performed, each sector in Eq. (8.103) provides the following contribution to the final integral over energies: where the auxiliary function T We note that the reason that CC 1 is present in Eq. (8.106) is that it still has to act on the phase space; its action on the matrix element has already been accounted for and resulted in the factorized form of Eq. (8.105) and the appearance of the triple-collinear splitting function in Eq. (8.106). We use Eq. (8.105) to write In what follows, we discuss the integration over energies in Eq. (8.107). Our goal is to change variables in such a way that the argument of the hard matrix element becomes z · 1; once this happens, Eq. (8.107) becomes a convolution of a hard matrix element with a splitting function. Although, in principle, changing variables in an integral is straightforward, it turns out that it is beneficial to do so in different ways in the four terms that appear in Eq. (8.107); for this reason, we consider them separately. We emphasize that since the individual contributions to Eq. (8.107) diverge, it is important to keep dimensional regularization in place until the end of the computation.
We begin with the term that contains the identity operator I and change the variables as follows: with r ∈ (0, 1) and z ∈ (0, 1). We note that the lower integration boundary for z can be taken to be z = 0 because F L M (z · 1, 2) always appears in this contribution. As we already discussed several times, this automatically cuts off the integral over z at a proper minimal value. With this in mind, we write 9 Next, we consider the S 5 operator. It describes the limit where E 5 → 0 at fixed E 4 . Calculating this limit with the parametrization in Eq. (8.108) mixes z and r and, therefore, is inconvenient. A better way is to change the parametrization.
We choose where r ∈ (0, 1). In principle, we should use z > z min but, since z enters the hard matrix element, we can extend all the integrals to z ∈ (0, 1). We find Note that the r 2 prefactor ensures that the r → 0 limit of the square bracket exists.
We can also use the change of variables in Eq. (8.110) for terms with operators SS and SS S 5 . The only difference is that since in those terms z does not appear in F L M , we have to keep the lower integration boundary at z = z min . We write . (8.112) Also in this case, the (1 − z) 4 prefactor ensures the existence of the z → 1 limit of the term in the square bracket. The term with an operator SS S 5 in Eq. (8.107) is obtained from Eq. (8.112) by taking the r → 0 limit in the expression in square brackets.
To proceed, it is convenient to define two z-dependent functions and a constant We can now write the result for the integral that we are interested in using A 1,2 (z) and A 3 . We find 10 (1, 2) . (8.116) This integral can be re-written in such a way that all the z → 1 singularities are regulated by plus-prescriptions. We have already discussed how this can be done several times; for this reason, we do not repeat this discussion again and only present the result. It reads 1+4 , The functions A (k) 1,2 (z) and the constant A (k) 3 are given in Eqs. (8.113) and (8.115). We have also used the following notation: (8.119) 10 At this point, we restore the sector label.
We are now in a position to write the contribution of this term in the final form All the terms in Eq. (8.120) can be expanded in power series in the dimensional regularization parameter . The functions R (k) (z) and the constants R (k) + and R (k) δ are calculated numerically.

Pole cancellation and finite remainders
We are now in position to discuss the final result for the NNLO QCD contribution to the cross section. We consider dσ NNLO = dσ RR + dσ RV + dσ VV + dσ ren + dσ CV .
All the different contributions to Eq. (9.1) were considered in the previous sections. It should be clear from these discussions that the result for the NNLO cross section is given by a linear combination of integrated matrix elements with different multiplicities, which may or may not be convoluted with generalized splitting functions. Since, for well-defined observables, the cancellation of soft and collinear divergences occurs point-by-point in the phase space, contributions proportional to F L M (1,2,4,5), F L M (1,2,4), F L M (z · 1, 2, 4) etc. must be separately finite. For this reason, it is convenient to present the result for the NNLO QCD contribution to the cross section as a sum of seven terms (1,2,4,5) (1,2,4,5) This contribution is the only one that involves the matrix element for qq → V + gg. We repeat here the result, already given in Eq. (8.11) dσ NNLO F L M (1,2,4,5) = F s r c r L M = (i j)∈dc (1,2,4,5) .
It follows that dσ NNLO F L M (1,2,4,5) is expressed through a combination of nested soft and collinear subtractions and can be directly computed in four dimensions. (1,2,4) We continue with terms that involve F L M (1,2,4). They are present in the double-real contribution, Eqs. (8.33) and (8.90) and in the real-virtual contribution, Eq. (7.8); they are also found in terms that appear due to ultraviolet Eq. (5.1) and collinear Eq. (5.8) renormalizations of the next-to-leading order cross section. Extracting these terms, we observe that all the 1/ singularities cancel out. The finite remainder reads

Terms involvingÔ N L O F L M
Next, we collect the finite remainders of the one-loop and two-loop virtual contributions to the qq → V process. These contributions appear in the real-virtual, the double-virtual, the collinear subtraction and the ultraviolet renormalization. Upon combining them and expanding the resulting contributions in , we obtain dσ NNLO (9.6) 9.4 Terms of the form P 1 ⊗ dσ ⊗ P 2 Terms of the type P 1 ⊗dσ ⊗P 2 , where P 1,2 are some splitting functions, appear in the double-real contribution as well as in the collinear renormalization. Combining all the relevant terms, we find 123 9.5 Terms of the form P ⊗ dσ These terms appear in the double-real, real-virtual, collinear subtraction and ultraviolet renormalization contributions. We note that starting from O(1/ ), the part of the double-real contribution related to the integral of the triple-collinear splitting function is only known numerically; see the discussion in Sect. 8.7.
Combining all the terms, we observe analytic cancellation of the poles up to 1/ 2 . For the 1/ poles and the finite part, it is useful to split the contribution into a scale-independent and a scale-dependent term We also introduce an expansion of the functions R (k) (z) and the constants R (k) + , which were introduced in Section 8.7, in powers of The scale-independent term reads (9.10) The scale-dependent term reads To arrive at these results, we check the cancellation of 1/ poles in dσ NNLO F L M (z·1,2) and then, assuming that the cancellation is exact, deduce the analytic form of R (0) (z) and R All the different contributions to the NNLO cross section produce terms proportional to F L M (1, 2). These include constants R (k) δ originating from the triple-collinear splitting function, which, as mentioned in the previous subsection, are only known numerically. As before, we introduce an expansion in for these constants.

Numerical results
Having described the subtraction procedure in some detail, we will now study how well it works in practice. We have implemented it in a partonic Monte Carlo program to compute NNLO QCD corrections to the production of a vector boson γ * in proton-proton collisions. 11 The calculation is 11 We remind the reader that we only consider the qq annihilation channel and restrict ourselves to gluonic corrections in this paper.  We begin by comparing the analytic result for the NNLO QCD correction to the pp → γ * → e + e − + X cross section, which we extract from Ref. [87], and the result of the numerical computation based on the formulas reported in the previous section. We wish to emphasize that this comparison is performed using the NNLO contribution to the cross section, and not the full cross section at NNLO, which would have included LO and NLO contributions as well. We take 14 TeV as the center-of-mass collision energy. We include lepton pairs with invariant masses Q in the range 50 GeV < Q < 350 GeV and take μ = 100 GeV for the renormalization and factorization scales. We obtain the NNLO corrections to the cross sections dσ NNLO = 14.471(4) pb, dσ NNLO analytic = 14.470 pb, (10.1) where the first result is ours and the second is extracted from Ref.
[87]. The agreement between the two results is quite impressive; it is significantly better than a permille. To further illustrate the degree of agreement, we repeat the comparison using the kinematic distribution dσ NNLO /dQ, shown in Fig. 2. In the upper pane of Fig. 2, we see a perfect agreement of analytic and numerical results for a range of Q-values where the cross section changes by five orders of magnitude. The ratio of numerical and analytic cross sections is shown in the lower pane of Fig. 2. We see that the agreement is between a fraction of permille and a few percent for all values of Q considered. We reiterate that we plot the NNLO correction to the differential cross section and not the full cross section at NNLO. Given that the NNLO contribution changes the NLO result by about 10%, the permille to percent precision on the NNLO correction leads to almost absolute precision for physical cross sections and simple kinematic distributions. We will further illustrate this point below. Before doing so, we note that we found a similar level of agreement for individual color structures and for individual contributions to the final result. We also note that although we report results for a single scale choice here, using the results in the previous sections and the known amplitudes for qq → e + e − + X , it is easy to check analytically the scale dependence of our result against the one reported in Ref. [87]. Full agreement is found.
As we mentioned in the Introduction, one of the important issues for current NNLO QCD computations is their practicality. For example, with the increasing precision of Drell-Yan measurements, one may require very accurate theoretical predictions for fiducial volume cross sections. It is then important to clarify whether a given implementation of the NNLO QCD corrections can produce results that satisfy advanced stability requirements and, if so, how much CPU time is needed to achieve them.
To illustrate this aspect of our computational scheme, we show the rapidity distribution of the dilepton pair, the rapidity distribution of a lepton, and the lepton transverse momentum distribution in Fig. 3. The plots on the left and on the right provide identical information: the upper panes show nextto-leading and next-to-next-to-leading order predictions for the respective observable, and the lower panes the ratio of the NNLO to NLO distributions. The difference between the plots on the left and the plots on the right is the CPU time required to obtain them; it changes from O(10) CPU hours for the plots on the left, to O(100) CPU hours for the plots on the right. The different run times are reflected in different binto-bin fluctuations seen in both plots. The bin-to-bin fluctuations for the two rapidity distributions are at the percent-level or better in the plots on the left, and they become practically unobservable in the plots on the right. The situation is slightly worse for the transverse momentum of the lepton. However, this observable is rather delicate in the γ * case, as each bin receives contributions from a large range of invariant masses. The introduction of a Z boson propagator will localize the bulk of the cross section in a much smaller invariant mass window, and lead to improved stability in this case. 12 Nevertheless, the results shown in Fig. 3 imply that the numerical implementation of our subtraction scheme allows for high precision computations, while also delivering results that are acceptable for phenomenology even after relatively short run times.

Conclusions
In this paper we described a modification of the residueimproved subtraction scheme [10,11] that allows us to remove one of the five sectors that are traditionally used to fully factorize singularities of the double-real emission matrix elements squared. The redundant sector includes correlated soft-collinear limits where energies of emitted gluons and their angles vanish in a correlated fashion. Once this sector is removed, the physical picture of independent soft and collinear emissions leading to singularities in scattering amplitudes is recovered and the bookkeeping simplifies considerably.
Using these simplifications, we reformulated a NNLO subtraction scheme, based on nested subtractions of soft and collinear singularities that, in a straightforward way, leads to an integrable remainder for the double-real emission cross section. The subtraction terms are related to cross sections of reduced multiplicity; they can be rewritten in a way that allows us to prove the cancellation of 1/ singularities independent of the hard matrix elements. Once singularities cancel, the NNLO QCD corrected cross section is written in terms of quantities that can be computed in four dimensions.
The second is the splitting functionP (1) qq . We emphasize that, for our case, we need the NLO splitting function for a continuous quark line entering the hard matrix element after the emission of two gluons. We then only have to consider the non-singlet NLO splitting function from which the contribution of identical quarks is subtracted. We extract the relevant information from Ref. [87].P (1) qq in Eq. (2.11) must then be identified witĥ P (1) qq, NS (x) = C A C F 3π 2 (1 + x) − 124x − 19 18 We now list the definitions of the various splitting functions used in the text. In these formulas, we use the definition (a) The tree-level splitting function used in the NLO computation is defined as (A.6) (c) The splitting function P qq ⊗ P qq NLO CV (z) reads P qq ⊗ P qq NLO CV (z) = C 2 F 6D 0 (z) + 8D 1 (z) (e) The tree-level splitting function used in the real-virtual contribution is defined as P qq,RV,1 (z) = C F 2 δ(1 − z)L 1 − D 0 (z) (A.10) (f) The one-loop splitting function used in the real-virtual contribution reads (A.11) (g) The splitting function P qq,R R 1 reads (h) The splitting function P qq,R R 2 (z) reads (i) The splitting function P qq,R R 3 (z) reads (A.14) (j) The splitting function P qq,R R 4 (z) reads P qq,R R 4 (z) = C F δ(1 − z)L 2 1 − 2D 1 (z) The splitting function P qq,R R 5 (z) reads (l) The splitting function P qq,R R 6 (z) reads P qq,R R 6 (z) = C 2 F 4D 1 (z) − 12D 2 (z) + (m) The splitting function P qq ⊗ P qq R R (z) is defined as We also require the following constants: and discuss its parametrization.
We take E 4 = E max x 1 , E 5 = E max x 1 x 2 with x 1 ∈ (0, 1) and x 2 ∈ (0, 1) and write the phase space in Eq. We have introduced the following notation: In Eq. (B.13), the term w 14,15 corresponds to the triplecollinear sector where singular radiation occurs along the direction of the incoming quark with momentum p 1 , w 24,25 to the triple-collinear sector where singular radiation occurs along the direction of the incoming antiquark with momentum p 2 , and w 14,25 and w 15,24 to the double-collinear sectors.