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 [Formula: see text] gluonic corrections to the Drell-Yan process of [Formula: see text] annihilation into a lepton pair. We show that this framework leads to fast and numerically stable computation of QCD corrections.

Introduction

One of the most important recent advances in perturbative QCD was the discovery of practical ways to perform fully differential next-to-next-to-leading order (NNLO) QCD computations for hadron collider processes. These methods, which include antenna [1-9], residue-improved sector-decomposition [10-12] (see also [13]) and projection to Born [14] subtraction schemes, as well as [15, 16] and N-jettiness [17, 18] slicing methods, were used to perform an impressive number of NNLO QCD computations relevant for LHC phenomenology [14, 15, 18–62].1

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-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 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 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 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 protonsWe are interested in computing the differential cross section for the process in Eq. (2.1)where is the finite partonic scattering cross section, are parton distribution functions and 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 , the ultraviolet renormalization reduces to the following () redefinition of the strong coupling constant:where is the Euler–Mascheroni constant, andis the leading-order (LO) QCD -function.

Collinear divergences associated with initial state QCD radiation are removed by a redefinition of parton distributions. In the scheme, this amounts to the replacementwhere stands for the convolutionand are the Altarelli–Parisi splitting functions.

As we already mentioned, we focus on gluonic corrections to the annihilation channel,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 annihilation into a vector boson is expanded in series of . We writewhereVarious 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 readwhereand the relevant splitting functions are provided in Appendix A.

The cross section 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 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 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-to-leading 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 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 real-emission contribution. It refers to the processwhere 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) aswhere s is the partonic center-of-mass energy,andIn Eq. (3.4), 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 is an (infra-red safe) observable that depends on kinematic variables of all particles in the process. Also, is an arbitrary auxiliary parameter that has to be large enough to accommodate all possible kinematic configurations for . 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 in Eq. (3.2) is attempted. To this end, we introduce two operators that define soft and collinear projectionswhere and is a unit vector that describes the direction of the momentum of the ith particle in -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 is a short-hand notation for a combination of soft and collinear projection operatorsNote that in Eq. (3.6) soft and collinear projection operators act on , 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 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 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 asIt is easy to see that the first term in Eq. (3.8) vanishes.4 Indeed, in the limit when the gluon becomes soft, we findwhere is the bare QCD coupling, is the QCD color factor, and is closely related to the LO cross sectionThe action of the collinear operators on the givesSince for head-on collision , , we findthis implies .

Hence, the only term that we need to consider in Eq. (3.8) is the collinear subtractionWe will consider the action of the operator on and then infer the result for the operator . First, we find the collinear limitWe note that a new variable is introduced in Eq. (3.14). The notation implies that in the computation of , cf. Eq. (3.10), the momentum is replaced with everywhere, including the energy-momentum conserving -function. We also used to denote the splitting functionTo simplify , we integrate over the emission angle of the gluon , rewrite the integration over its energy as an integral over z and express in terms of the renormalized coupling . After straightforward manipulations we findwhere and we introduced the short-hand notationWe note that in Eq. (3.16) integration over z leads to divergences caused by the soft 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 and split into a piece that is singular at and a regular pieceWe also note that we can extend the integral over z in Eq. (3.16) to since if , will vanish because there is not enough energy to produce the final state. We will use this fact frequently in our NNLO analysis. We writeThe 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:where . It is now straightforward to rewrite Eq. (3.16) in such a way that all soft, , singularities are regulated using the plus-prescription. We use the fact that we are in the center-of-mass frame of the incoming pair, so that . We findThe splitting function in Eq. (3.21) reads5 where is the LO Altarelli–Parisi splitting kernel, see Eq. (A.1), and is defined asThe result for is obtained by a simple replacement in Eq. (3.21). Putting everything together, we find the following result for the real-emission cross section: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 lower-multiplicity 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 obtainwhere 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 obtainwhereand 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 in Eq. (3.26), we find the final formula for the NLO QCD contribution to the scattering cross section for . It readsIt follows that the NLO cross section is computed as a sum of low-multiplicity terms, including those where or are convoluted with particular splitting functions, and the subtracted real-emission term described by . 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 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 (double virtual), from one-loop corrections to the process with an additional gluon in the final state (real-virtual), and from the tree-level process with two additional gluons in the final state (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-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 . 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 .

The NNLO computation: ultraviolet and PDF renormalization

In this section we study the contributions to coming from the ultraviolet and the collinear renormalization, beginning with the former. As mentioned previously, because the process is driven by a conserved current which is independent of , the ultraviolet renormalization contribution is very simple. Combining Eq. (2.10) and the first two lines of Eq. (3.26), it is straightforward to obtainWe proceed to the collinear subtraction. Rewriting Eq. (2.10) to make the convolutions explicit, we obtainTerms that involve convolutions of the various splitting functions with 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.

It is less straightforward to rewrite and in a form convenient for combining them with other contributions to . We focus on , and consider the effect of the convolution on the first two lines of Eq. (3.26).

First, we consider the term proportional to 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 which becomes once the momentum is changed to . On the other hand, the soft remainders of the collinear subtracted terms scale as , with . Hence, in the calculation of , the corresponding contribution scales with z either as or as . Therefore, we have to computewhereFor future convenience, we rewrite Eq. (5.3) aswhere the splitting function is defined in Eq. (A.6).

The other two terms that we need involve convolutions of splitting function with , cf. Eq. (3.26). The first term can be written as a double convolutionwhere the splitting function is defined in Eq. (A.7). The second term is the left–right convolution,Combining all these terms we findEach 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 and write it as an expansion in the renormalized strong coupling constantwhere 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 writewhere in the second line we removed the renormalization contribution that is already accounted for in Eq. (5.1).

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 poles from . Soft and collinear singularities of a generic scattering amplitude are given in Ref. [70]. Using these results, we rewrite asThe sum of the last two terms in Eq. (6.3) is a finite remainder of the contribution to the virtual corrections once its divergent part is written in a form suggested in Ref. [70]. More specifically, is the finite remainder of the one-loop amplitude squared, while is the genuine two-loop finite remainder. Note that, contrary to and , is scale-dependent; the scale-dependent contribution readsAs follows from Eq. (6.3), the singularities of the double-virtual corrections are proportional to the leading-order contribution and to the finite part of the virtual corrections to the NLO cross section , 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 writeWe remind the reader that, according to our notation, soft- and collinear-projection operators in Eq. (7.1) do not act on the phase space of the gluon . 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 findWe need to integrate Eq. (7.2) over the phase space of the gluon . This can easily be done, with the resultNote 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 for this purpose, cf. Eq. (3.3).

The second term that we need to consider is the soft-regulated collinear subtraction termWe will only discuss the collinear projection operator ; the contribution of 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 findwhere the splitting function is given in Eq. (A.8). The soft-collinear limit is easily obtained by taking the collinear approximation in Eq. (7.2). We findIntegrating over emission angles of the gluon and rewriting the result through plus-distributions, following the discussion in Sect. 3, we obtain a convenient representation for the soft-regulated collinear subtraction term. It readswhere the two splitting functions are defined in Eqs. (A.10) and (A.11). We replace with in Eq. (7.7) to obtain the result for .

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 obtainWe 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 singularities when combining Eq. (7.8) with other contributions to . To this end, it is useful to make all the singularities explicit in Eq. (7.8) by writing in terms of and , cf. Eq. (3.25), and asWe used and denoted a finite remainder which does not depend on the scale by .

The NNLO computation: double-real emission

General considerations

In this section, we discuss the double-real emission contributions to . 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 in these regions. In the case of the two-gluon emission in 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 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 writewhere, as in Sect. 3, is the phase space for the final state V, including momentum-conserving delta-function. The gluon phase space elements 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 exampleimplies that the operator acts on and on the phase space element .

We begin by extracting soft singularities from the double-real process, largely repeating what we have done at next-to-leading order.6 We writeIn Eq. (8.5), the last term is soft-regulated, in the second term gluon is soft and soft singularities associated with are regulated, and in the first term both and 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 writewhere the label i runs through the elements of the set . 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 in such a way that when they multiply the matrix element squared, the resulting expression is only singular in a well-defined subset of limits. For example, in the partition 14, 15 collinear singularities in only occur when gluons 4 and/or 5 are emitted along the direction of the incoming quark or when their momenta are parallel to each other. Similarly, in the partition 24, 25, the singularities occur when momenta of and/or are parallel to or to each other. In the partitions 14, 25 and 15, 24, singularities only occur when momenta of and are collinear to and or and , respectively. Apart from these requirements, the specific form of 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 , 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) readswhere, 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 replacement. However, because of the energy ordering , we no longer have the 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 explicitly in what follows.

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 aswhere the soft-regulated, single-collinear term readsthe soft-regulated, triple-collinear terms readsand the fully regulated term readsFor the process under consideration, and . The above results are obtained by combining soft-regulated expression for with multiple partitions of unity for the angular projections, and the understanding of which collinear divergences can appear in each partition and sector. One can easily check, starting from Eq. (8.8), that the collinear projection operators add up to an identity operator for each partition and each sector.

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 and collinear projection operators is that in the partition i4, j5 no other collinear singularities appear, thanks to the damping factor . Similarly, taking e.g. the sector a of the triple-collinear partition , we writebecause in this case a singularity can only occur either in a triple-collinear limit or if at fixed .

It is worth pointing out a few things in connection with Eq. (8.8).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 procedure that we used to write Eq. (8.8) is, in principle, process- and phase space parametrization-independent. 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.

The first term in Eq. (8.8), , is the double-soft subtraction term. It contains unregulated soft and collinear singularities and cannot be directly expanded in . However, it only involves the tree-level matrix element which means that emitted gluons decouple from both the hard matrix element and the phase space constraints. When integrated over gluon energies and angles, this term gives rise to poles, .

The second term in Eq. (8.8), , is the double-soft regulated, single-soft subtraction term. It contains and matrix elements of lower multiplicity. This term still contains unregulated singularities that occur when the momentum of the gluon becomes collinear to the collision axis or to the direction of . When integrated over gluon energies and angles, this term gives rise to poles, with .

The term in Eq. (8.9) is the soft-regulated single-collinear subtraction term. Note that thanks to the damping and factors only one kind of collinear singularity per term is present. involves , 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 and poles.

The term 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 have the “double-convolution” structure. The term contains poles in contributions of triple-collinear partitions, and poles in contributions of double-collinear partitions.

The term 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 .

The double-soft subtraction term

We begin with the discussion of the first term in Eq. (8.8), . 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 writewhereand [75]withAt this point, we stress again that the hard matrix element 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 energy-momentum 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 as in Eq. (8.3).

To satisfy constraints on gluon energies, we parametrize them asWritten in these variables, the eikonal factor becomesThe important point is that the dependence on the overall energy scale factorizes and that the remaining (complicated) function is independent of energies of the incoming partons. We also use the parametrization of energies Eq. (8.18) in the phase space to obtainFor 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-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,where . Numerical values of the coefficients 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 does not require any further regularization; we will later combine it with other contributions with tree-level kinematics to cancel the singularities explicitly.

Table 1

Coefficients of the expansion of the double-soft projected real-emission contribution. Full results are given in the first row. Results for individual color factors are given in the second and third rows. Numerical errors are such that their contribution to the final result is below the per mille level

c4SS	c3SS	c2SS	c1SS	c0SS
5.55554(2)	-11.73653(7)	-7.3253(7)	-20.796(5)	-54.65(7)
1.999995(8)CF2	-5.54530(5)CF2	1.1077(3)CF2	1.522(1)CF2	1.961(4)CF2
0.499999(2)C_AC_F	-0.46960(1)C_AC_F	-2.3236(1)C_AC_F	-5.876(1)C_AC_F	-14.52(1)C_AC_F

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 readsNote that this contribution implicitly depends on and , and the hard matrix element that appears in still contains collinear singularities that arise when the momentum of gluon becomes parallel to the momenta of the incoming partons or to the direction of . These divergences will have to be extracted and regulated.

We start by computing the soft limit for the gluon . We findSince the gluon decouples from the hard matrix element, we can integrate over its momentum. We findwhereandWe need to simplify Eq. (8.24) because it still contains collinear singularities that appear when the momentum of gluon becomes parallel to the collision axis. To extract them, we writeWe reiterate that according to our notational conventions, the collinear projection operators do not act on the phase space element of the gluon 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 ; the contribution of the projection operator is then obtained by analogy.

First, we consider how acts on . Using , and taking the limit on , we obtainIntegrating over the energy and angle of the gluon we arrive atwhere . The splitting function operator is defined by means of the following equation:whereand the splitting function is given in Eq. (3.15).

Note that in Eq. (8.29) can be directly expanded in powers of since all the singularities are regulated. The only problem that needs to be addressed is the fact that the integration over z does not start at , as is the case for the convolutions. The lower integration boundary must be kept in Eq. (8.29) because of the subtraction term . Indeed, if is replaced with zero, the integration over the gluon energy for this term extends to the region , in contradiction with the original phase space parametrization. The extension of the integration region in Eq. (8.29) is accomplished following steps discussed in the context of the NLO QCD computation in Sect. 3. Effectively, this leads to a redefinition of the splitting functionwhere is given in Eq. (A.12). We note that the contribution of the collinear operator to the second term in Eq. (8.27) is computed in a similar way; the computation leads to the same result as in Eq. (8.29) up to an obvious replacement .

Putting everything together, we find the final result for the double-soft regulated single-soft singular contribution to ,

The single-collinear term

Next, we consider the soft-regulated, single-collinear contribution to We need to rewrite Eq. (8.34) in such a way that extraction of the remaining collinear singularities becomes straightforward. We note that 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 . We begin with the partition 14, 25. For the first term, we need to computewhere does not depend on the momentum of gluon anymore. To further simplify Eq. (8.35), note that collinear and soft projection operators commute with each other and thatThis implies that we can drop the term in Eq. (8.35). We use the collinear limit for obtained by a straightforward generalization of Eq. (3.14). We define and obtainThe function was introduced in Eq. (8.31).

We now consider the phase space. According to Eq. (8.34), acts on the phase space element . We introduce to getwithIn this parametrization, the action of implies replacing with . Putting everything together, we obtainNote that, similar to the NLO case, the lower boundary is not important when integrating since for there is no sufficient energy in the incoming partons to produce the required final state.

Next we consider the action of the projection operator. Following the preceding discussion and using , we obtainThe operator was introduced in Eq. (8.30) andThe sum of Eqs. (8.40) and (8.41) gives the required result for the collinear sector 14, 25.

The partition 15, 24 is obtained from the results for 14, 25 after a few obvious replacements. We findandWe can now combine the contributions of the two double-collinear partitions. In doing so, it is convenient to always denote the “resolved” (i.e. the non-collinear) gluon by . Out of the four terms that we need to combine, two correspond to the collinear emission along the direction of the incoming quark and two along the direction of the incoming antiquark . We consider terms that belong to the former category first.

When combining results, it is important to realize that . We will denote it by . After straightforward manipulations we findNote that the lower integration boundary in this formula should be but we can extend the integration region to , without making a mistake. This is so because every time appears in the integrand, the 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 does not appear, -functions require that . We also note that, thanks to explicit subtractions and constraints due to -functions, each term in Eq. (8.45) vanishes if or . Finally, we stress that and commute since they act on different variables.

We can write a similar equation for the sum of the two terms where the collinear gluon is emitted along the direction of the antiquark . Finally, putting everything together, we obtain the contribution of the double-collinear partitions to . We findNote 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 consider the triple-collinear partition 14, 15 and study the contribution of single-collinear limits in Eq. (8.34). We begin with sector (a). The relevant expression readsThe calculation is identical to the case of the double-collinear partition except that we need to account for the constraint that defines sector (a) when integrating over the angle of gluon . Writing and , and taking , we findUsing this result, we obtainA similar calculation for the sector (c) givesIn parallel to the case of the double-collinear partitions, it is again convenient to always call the resolved gluon . We combine contributions of sectors (a) and (c), renaming where appropriate, and we obtainWe 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 and becoming collinear to each other. We considerand start with the discussion of how the collinear projection operator acts on . We findwhere , i.e. the hard matrix element must be taken in the collinear limit. The splitting functions areand z is the fraction of the total momentum carried by gluon ,The vector is a normalized transverse vectordefined by the following decomposition:where and .

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 and become collinear. We write the four-momenta of and aswhere , , and . Our goal is to parametrize the phase space in such a way that explicit averaging over directions of can be performed. The phase space parametrization for sectors (b) and (d) can be written aswithwhereHere, corresponds to the 4||5 collinear limit, , and is related to the azimuthal angle . Further details about the parametrization as well as expressions of scalar products in terms of and can be found in Appendix B. In this parametrization, the vector reads8 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 can be performed since neither nor appear in the hard matrix element. We defineUsingandwe findIt 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 collinear projection operator in sector (b). We obtainwhereIt follows from Eq. (8.68) that we need to know how an operator acts on . We recall that the action of on energy variables implies that at fixed . Computing the soft limits is simple and standard except for the spin-correlated part that we address below. In principle, we need to know three soft limits and . However, sincewe only need to consider . 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 and , is a valid polarization vector of the gluon with momentum , in the collinear 4||5 approximation. For this reason, the soft limit of follows from the standard soft limit of the amplitude for , not averaged over gluon polarizations. We obtainCollecting all the soft limits, we findwhere z is defined in Eq. (8.55). This impliesWe 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 to and z.

We find that the integration region splits into two regions (cf. Fig. 1)Following this separation, we split the integral into two partswhere the integral corresponds to the region and the integral to the region ; see Fig. 1. We obtain an integral representation for starting from Eq. (8.72), changing variables in the first term, and in the second term in Eq. (8.72) and, finally, renaming . We obtainTo compute , we notice that only the soft term of Eq. (8.72) contributes. Again renaming and using Eq. (8.71), we obtainwhere and we expressed as an integral over energy and the angular integration variable . We have also integrated over angular variables that do not appear in the hard matrix element and in the splitting function.

Fig. 1

Integration region for the change of variables. The colored triangle is the allowed region. The blue region “A” is the “physical” one, i.e. the one which is not removed by a phase space -function inside . The orange region “B” only contributes to the soft limit, since there no -function from is preventing it. Lines of fixed are shown in solid red (for ) and dashed orange (for ). In dot-dashed blue, lines of constant z are shown. In the “physical” region, only the condition is relevant. In the “B” region, we also have to impose , to prevent the integration to go outside the triangle (see the intersection of blue and orange lines)

Finally, we consider sector (d). We need to computeThe calculation is similar to what we just described for sector (b), apart from the following modifications of the integration boundaries: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 and the z-integration decouples from the rest and can be performed independently of the hard matrix element. In , it yieldswhere we used and the constants , are defined in Eq. (A.19). In the integral 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 findwith defined in Eq. (A.20).

Summing double- and triple-collinear partitions

Summing up the above results, we obtain an intermediate representation of . We write it as a sum of four terms,These terms are defined as with

The single-collinear term: extracting the last singularities

The four contributions to 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 writeThe 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 findwhere z is defined in the usual way . Taking this into account, after tedious but straightforward calculations we arrive at In the above equation, we used the following notation:The relevant splitting functions are defined in Appendix A. Also,andWe note that in Eq. (8.90) the first curly bracket is fully regulated, while the second contains subtraction terms. Note also that since , if we are only interested in the poles, we can substitute in Eq. (8.90). We also note that, for the process of interest, terms that contain can easily be integrated over the relative angles of the gluon with respect to the collision axis. We findWe will use these results when presenting the final formula for the double-real contribution.

The double-unresolved collinear limit: double collinear

We now turn to the term and begin by considering the contribution of the double-collinear partitions. It readsNote that, following our notational convention, the collinear projection operators act on the phase space elements and .

We begin with the term. Introducingand calculating collinear limits we obtainSince the momenta of gluons and decouple from each other, we findAs the result, the original expression simplifiesPerforming the angular integrations and accounting for the hierarchy of energies , we obtainwhere, as usual, the z integrals do not need a lower cut-off whenever z is present in .

The term with collinear operators in Eq. (8.95) can be simplified in a similar way. Combining the two contributions, we obtainWe can rewrite Eq. (8.101) to ensure that all singularities that appear in z and integrals are regulated with the plus-prescription. This gives the final result for the double-collinear contributionThe relevant splitting functions are given in Appendix A.

The double-unresolved collinear limit: triple collinear

In this section, we consider the contribution of the triple-collinear partitions to . It readsThis contribution always contains the triple-collinear projection operator that acts on the hard matrix elements. For , this gives, schematically,where and 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 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 ; this corresponds to taking in Eq. (8.103).

To proceed, we note that the damping factors in Eq. (8.103) can be removed since the collinear projection operator acting on them yields 1. Next, we need to study the triple-collinear 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 in Eq. (8.105) is defined asWe note that the reason that 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 writewhere in the collinear approximation.

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 ; 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 and change the variables as follows:with and . We note that the lower integration boundary for z can be taken to be because 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 write9 Next, we consider the operator. It describes the limit where at fixed . 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 choosewhere . In principle, we should use but, since z enters the hard matrix element, we can extend all the integrals to . We findNote that the prefactor ensures that the limit of the square bracket exists.

We can also use the change of variables in Eq. (8.110) for terms with operators and . The only difference is that since in those terms z does not appear in , we have to keep the lower integration boundary at . We writeAlso in this case, the prefactor ensures the existence of the limit of the term in the square bracket. The term with an operator in Eq. (8.107) is obtained from Eq. (8.112) by taking the limit in the expression in square brackets.

To proceed, it is convenient to define two z-dependent functions and a constantWe can further simplify the function if we realize that the limit of the triple-collinear splitting function is homogeneous in . This implies that the limit of the two functions in the formula for in Eq. (8.113) are related and can be combined. Changing variables in the first term on the r.h.s. of the integral for , we obtainSince the term in the square bracket no longer depends on r after the limit is taken, we can perform the r integrations in the first line to getWe can now write the result for the integral that we are interested in using and . We find10 This integral can be re-written in such a way that all the 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 readswhereThe functions and the constant are given in Eqs. (8.113) and (8.115). We have also used the following notation:We are now in a position to write the contribution of this term in the final formAll the terms in Eq. (8.120) can be expanded in power series in the dimensional regularization parameter . The functions and the constants and 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 considerAll 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 , , 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 termswhich are individually finite. Each of the individual terms in Eq. (9.2) has a subscript that indicates the highest multiplicity matrix element that it contains. Below we collect all the different contributions to and present finite remainders for terms with different multiplicities. For simplicity, we fix the arbitrary parameter .

Terms involving

This contribution is the only one that involves the matrix element for . We repeat here the result, already given in Eq. (8.11)It follows that is expressed through a combination of nested soft and collinear subtractions and can be directly computed in four dimensions.

We continue with terms that involve . 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 singularities cancel out. The finite remainder readswhere

Terms involving and

Next, we collect the finite remainders of the one-loop and two-loop virtual contributions to the 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

Terms of the form

Terms of the type , where are some splitting functions, appear in the double-real contribution as well as in the collinear renormalization. Combining all the relevant terms, we find

These terms appear in the double-real, real-virtual, collinear subtraction and ultraviolet renormalization contributions. We note that starting from , 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 . For the poles and the finite part, it is useful to split the contribution into a scale-independent and a scale-dependent termWe also introduce an expansion of the functions and the constants , which were introduced in Section 8.7, in powers of The scale-independent term readsThe scale-dependent term readsTo arrive at these results, we check the cancellation of poles in and then, assuming that the cancellation is exact, deduce the analytic form of and . These analytic results are then used in the scale-dependent term. Thus the only numerical contributions needed for the finite part are and .

All the different contributions to the NNLO cross section produce terms proportional to . These include constants 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.Furthermore, for the double-soft contribution we know the abelian constants of Eq. (8.21) analytically, but only have numerical results for the non-abelian constants, which are reported in Table 1. Thus, for each order in , we check the cancellation of terms numerically and then, assuming that the cancellation is actually exact, we deduce an analytic form of the triple-collinear splitting and double-soft constants at this order. This form is then used in determining the cancellation at lower orders in . Thus, the only numerical constants appearing in our final formula are and . The final formula reads

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 fully differential; we consider decays of the virtual photons to massless leptons and study NNLO QCD corrections to lepton observables. We extracted the relevant matrix elements from Refs. [81, 82] as implemented in [83], and from Refs. [84, 85]. For all computations reported below we employ the NNLO parton distribution functions from the NNPDF3.0 set [86].

We begin by comparing the analytic result for the NNLO QCD correction to the 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 as the center-of-mass collision energy. We include lepton pairs with invariant masses Q in the range and take for the renormalization and factorization scales. We obtain the NNLO corrections to the cross sectionswhere 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 , 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.

Fig. 2

Comparison of the NNLO QCD contribution computed in this paper with the analytic results in Ref. [87]

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 , 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 next-to-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 CPU hours for the plots on the left, to CPU hours for the plots on the right. The different run times are reflected in different bin-to-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.

Fig. 3

Upper panes Rapidity distribution of the vector boson, rapidity distribution of a lepton and distribution of a lepton at different orders of perturbation theory. Lower panes the ratio of NNLO/NLO prediction for a given observable. Plots on the left the runtime of CPU hours; plots on the right the runtime of CPU hours. Note that the dip in the ratio of NNLO/NLO lepton distribution at is a physical feature and not a fluctuation

Conclusions

In this paper we described a modification of the residue-improved 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 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.

Although we believe that this framework is applicable for generic NNLO QCD computations, in this paper, for the sake of simplicity, we studied dilepton pair production in quark–antiquark annihilation and computed gluonic contributions to NNLO corrections. We implemented our formulas in a numerical program and used it to calculate NNLO QCD corrections to the production cross section of a vector boson in hadron collisions with a sub-percent precision. We also showed that kinematic distributions, including the lepton rapidity and transverse momentum distributions, can be computed precisely and efficiently. We look forward to the application of the computational framework discussed in this paper to more complex processes, relevant for the LHC phenomenology.