Analytic results for color-singlet production at NNLO QCD with the nested soft-collinear subtraction scheme

We present analytic formulas that describe the fully-differential production of color-singlet final states in qq¯\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} and gg annihilation, including all the relevant partonic channels, through NNLO QCD. We work within the nested soft-collinear scheme, which allows the fully local subtraction of infrared divergences. We demonstrate analytic cancellation of soft and collinear poles and present formulas for the finite parts of all integrated subtraction terms. These results provide an important building block for calculating NNLO QCD corrections to arbitrary processes at hadron colliders within the nested soft-collinear subtraction scheme.


Introduction
Perturbative computations play an important role in the contemporary exploration of particle physics at the Large Hadron Collider (LHC). In particular, the lack of direct evidence for physics beyond the Standard Model suggests that further progress in particle physics will require a better understanding of hard hadron collisions and a confrontation of precise theoretical predictions with experimental results.
The quality of theoretical predictions for hard processes at the LHC has increased dramatically in recent years thanks to the advent of flexible methods for handling infrared singularities that led to the calculation of next-to-next-to-leading-order (NNLO) QCD corrections for sufficiently complex processes [1][2][3][4][5][6][7][8][9][10]. In spite of these successes, there are ongoing efforts to either simplify and improve existing methods or to devise "better" ones. For example, it is believed that higher efficiency of the slicing formalism of Refs. [6][7][8][9] may come from a deeper understanding of power corrections to soft and collinear limits [11][12][13]. Similarly, improved control of the interplay between soft and collinear dynamics may lead to the formulation of simple "minimal" local subtraction schemes [14,15]. Although there is no guarantee that any of these efforts will result in an absolute breakthrough in fixed-order calculations for LHC processes, it is plausible that these developments will lead to more efficient computational frameworks and enable precise phenomenological descriptions of complex multiparticle final states.
In Ref. [16] we have attempted to simplify the residue-improved subtraction scheme proposed in Ref. [2]. This subtraction scheme is very attractive because it is fully local, completely general and perfectly modular, so that the subtractions for a generic process are built from a relatively small set of basic ingredients. Its main disadvantages include a lack of physical transparency and a certain redundancy, as well as the numerical integration of the subtraction terms that may inadvertently impact its efficiency.
We have argued in Ref. [16] that QCD color coherence removes an interplay between angles and energies of soft and collinear particles in gauge-invariant QCD amplitudes, thus leading to a minimal number of subtraction terms that need to be considered. Perhaps more importantly, since soft and collinear singularities are not intertwined, it is possible to separate them cleanly, removing unnecessary redundancies of the subtraction procedure presented in Ref. [2]. As shown in Ref. [16], it appears to be advantageous to first subtract the double-soft singularities from the full amplitude and then, iteratively, remove the remaining ones. Once this is done, a transparent and physically appealing subtraction scheme is obtained. Moreover, this scheme appears to be very efficient, at least in the color-singlet production case that we have studied up to now [16,17].
Given the improved efficiency and inherent simplicity of the subtraction scheme developed in Ref. [16], it is natural to investigate whether one could obtain analytic results for the integrated counterterms. A successful completion of this task would lead to the formulation of the first subtraction scheme applicable at the LHC that is both fully local and under complete analytic control. Although it is hard to say to what extent these nice features are actually important in practice, we do hope that they will lead to a very efficient subtraction framework for completely generic processes.
It is easy to identify the major obstacles to obtaining fully analytic subtraction schemes. Indeed, any NNLO subtraction scheme involves three "double-unresolved" contributions whose integration is highly non-trivial. They are 1) the double-soft emission of two partons with energies E f 1 ∼ E f 2 √ s, where √ s is the center-of-mass energy of the partonic process; 2) the emission of two partons collinear to one of the incoming legs and 3) the emission of two partons collinear to one of the hard final-state legs. We note that as far as the numerical implementation of NNLO QCD corrections to a generic process is concerned, contributions 1) and 2) are the most problematic. Indeed, for any splitting process, the integrated contribution 3) is a number of the form a f / + b f , so that it can be calculated numerically once and for all. On the contrary, the integrated contributions 1) and 2) are functions of the relative angles between hard partons and the momentum fraction carried into the hard process, respectively. 1 Close to the end-points, these function may develop integrable singularities, which makes their numerical evaluation tedious.
In Ref. [18], some of us presented analytic result for the integrated double-soft subtraction term. In this paper, we will argue that a minor modification of the subtraction procedure described in Ref. [16] greatly simplifies the analytic integration of the triple-collinear subtraction terms. In fact, such an integration of all relevant triple-collinear subtraction terms has recently been performed in Ref. [19]. Thanks to these results, it is now possible to present a subtraction framework for the production of color-singlet particles at hadron colliders that is both fully analytic and fully local.
Although the production of color-singlets at NNLO QCD has been studied many times, including the development of public computer codes, even the simplest versions of these processes such as pp → Z and pp → H are quite useful to us because NNLO QCD corrections to these processes are known analytically [20,21]. This feature allows us to check all the non-trivial ingredients of our computational framework to a very high accuracy. We believe such a validation is important in view of its expected application to more complex cases in the future. Of course, it is also interesting to explore the performance of our subtraction scheme by considering a well-known process, where many benchmarks exist already.
Nevertheless, it should be clear that the goal of this paper is to present analytic formulas relevant for the production of generic color-singlet final states at a hadron collider, written in a way that will allow us to move beyond color-singlet production. For this reason we decided to avoid using simplifications that are particular to the cases of Drell-Yan or Higgs production.
The rest of this paper is organized as follows. In Section 2 we summarize the main features of the nested soft-collinear subtraction scheme of Ref. [16] and explain how we modified it to allow for an analytic integration of the triple-collinear subtraction terms. In Section 3 we provide formulas for the qq-initiated production of the color-singlet final state. In Section 4 we provide formulas for the gluon annihilation into a color-singlet final state. We discuss the validation of our results in Section 5 and conclude in Section 6. A large number of useful formulas are collected in appendices, as well as in an ancillary file attached to this paper.

Overview of the nested subtraction scheme
In this section we briefly review the method for computing NNLO QCD corrections described in Ref. [16] and explain how to modify it to simplify the analytic integration of the triple-collinear subtraction terms. We consider the production of a color-singlet state V in hadronic collisions. We write the (fiducial) cross section as where n f is the number of light flavors, dσ faf b is the partonic cross section in the f a f b channel, and we employ the following notation for parton distributions functions: f 0 = g and f ±1,±2,±3,±4,±5 = {d/d, u/ū, s/s, c/c, b/b}. Finally, O is a suitable infrared-safe observable that defines the fiducial volume. We consider the perturbative expansion of the partonic cross section and write Here q = 0 for quark-initiated processes and q = 2 for gluon-initiated processes, and we have suppressed the arguments of the functions on the right-hand side. We focus on the NNLO QCD contribution dσ NNLO faf b . It can be written as and dσ ren contains all contributions that originate from the renormalization of input parameters, such as the strong coupling constant α s and the parton distribution functions (PDFs). In Eq. (2.4), N is a normalization factor that takes into account color and spin averages, s is the partonic center-of-mass energy squared, dLips(V ) is the phase space for the final state V , and Here d is the dimensionality of space-time that we use as the regularization parameter, and E max is an arbitrary 2 energy scale that is introduced for future convenience. Each term in Eq. (2.3) is individually divergent. These divergences can either appear explicitly as poles in = (4 − d)/2 or be hidden inside phase-space integrals. The goal of any subtraction scheme is to extract these divergences and to arrive at the following representation of the NNLO contribution to the cross section are finite quantities that involve contributions with V and up to i partons in the final state. We will refer to dσ NNLO V +i , with i = 2, 1, 0, as terms with NNLO, NLO and LO kinematics, respectively.
In Ref. [16], we proposed a framework to achieve the separation of the NNLO contributions to the cross section as in Eq. (2.6). It is based on three ideas: • a multiparticle phase space can be decomposed into a sum of elements (partitions) such that for each partition only a well-definite subset of particles gives rise to collinear singularities upon integration over the phase space of final state partons. An early discussion of this idea can be found in papers on NLO QCD subtractions [22,23]; in the context of NNLO QCD calculations, it was reincarnated in Ref. [2]; • for each of these partitions there exists a phase space parametrization that allows the extraction of both soft and collinear singularities in a fully factorized form [2]; • thanks to gauge invariance and color coherence [24], soft and collinear singularities are not entangled in QCD amplitudes, in contrast to individual diagrams [16].
We argued in Ref. [16] that these three points allow us to follow the so-called FKS subtraction scheme [22,23], developed for NLO QCD computations, and to perform the required soft and collinear subtractions in a nested way. As a consequence, the computational framework becomes very transparent physically and quite efficient numerically. We will illustrate the main ideas of Ref. [16] by considering the double-real emission corrections to the Drell-Yan process qq → V as an example, focusing on the most complicated q(p 1 )q(p 2 ) → V + g(p 4 )g(p 5 ) channel. All other partonic channels can be dealt with along the same lines although the details can be somewhat different. 3 We write the 2 The only requirement on E max is that it should be at least as large as the maximum energy allowed by the momentum-conserving δ-functions in Eq. (2.4). For simplicity, throughout this paper we use s is the partonic center-of-mass energy. 3 Results for all the relevant channels are presented in the next sections.
corresponding cross section as where F LM,qq (1,2,4,5)  . Since the matrix element is symmetric with respect to the permutations of the gluons g 4 and g 5 , we can remove the 1/2! symmetry factor from N . Our goal is to extract singularities from Eq. (2.7). These singularities have different origins. There exist • a double-soft singularity that occurs when energies of the two gluons vanish in such a way that their ratio E 5 /E 4 is fixed; • a single-soft singularity that arises when E 5 vanishes at fixed E 4 . Note that due to the energy ordering in Eq. (2.7) the opposite limit, E 4 → 0 at fixed E 5 , cannot occur; • many different collinear singularities that appear when one or both gluons are emitted along the direction of the incoming quark or the incoming anti-quark, or when the momenta of the two gluons become parallel to each other.
We need to extract all these singularities in an unambiguous way. We explain how to do this in the next two subsections.

Extraction of soft singularities
As we explained in Ref. [16], it is convenient to begin by regularizing the double-soft singularity where S S is an operator that extracts the double-soft λ → 0 singular limit from F LM,qq . To make this statement precise, when the operator S S acts on F LM , it removes the fourmomenta of the gluons from both the energy-momentum conserving δ-function and the observable, and extracts the leading singular behavior from the matrix element squared. The result is S SF LM,qq (1, 2, 4, 5) = g 4 s,b Eik(1, 2, 4, 5) F LM,qq (1, 2), (2.10) where g s,b is the bare strong coupling and Eik (1,2,4,5) is the square of the eikonal factor derived e.g. in Ref. [25]. It is also given in Ref. [16] using notation that is identical to what we use in this paper. Also, F LM,qq (1, 2) is defined analogously to Eq. (2.8); it reads This tree-level matrix element squared integrated over the Born phase-space obviously provides the leading order result for the observable O.
We deal with the two terms on the right-hand side of Eq. (2.9) in very different ways. In the first term, thanks to Eq. (2.10), the hard matrix element decouples and only the eikonal factor needs to be integrated over the two-gluon phase space. In our original paper [16] this integral was calculated numerically but, since then, an analytic computation of this integral has become available [18]. The result reads 4 where we have defined and L = log s µ 2 . (2.14) In Eq. (2.12), the abelian part is known in a closed form (2.16) The equivalent results for gluon-initiated color singlet production can be obtained by simply replacing C F → C A in Eq. (2.12).
We now turn to the second term in Eq. (2.9) where the double-soft divergencies are already regularized but both the E 5 → 0 divergence at fixed E 4 and collinear divergences are still present. To extract them, we repeat the above procedure and subtract the E 5 → 0 singularities at fixed E 4 . We call the corresponding operator S 5 and write (1,2,4,5) .
When the operator S 5 acts on F LM,qq (1,2,4,5), it removes the gluon g 5 from the phase space and the observable, and extracts the leading singularity (2.18) We use the notation ρ ij = 1 − cos θ ij in Eq. (2.18), where θ ij is the relative angle between partons i and j. We have also introduced see Eq. (2.8). From here on, we will omit the subscript on M indicating the partonic process. It is clear that the second term in Eq. (2.17) has a simplified (i.e. independent of g 5 ) matrix element. The integration over the energy and angles of the gluon g 5 can therefore be performed, and the remaining infrared divergences in the matrix element for the process qq → V + g 4 can be dealt with in a way that is similarly to what is usually done in next-to-leading-order computations. On the other hand, the first term in Eq. (2.17) is now free of soft divergences but it still contains collinear singularities. We explain how to extract them in the next subsection.

Extraction of collinear singularities
In the previous subsection we extracted soft singularities from the double-real emission contribution by writing it as The procedure continues with the extraction of collinear singularities. This requires an additional step, similar to the energy ordering in Eq. (2.7). Indeed, we need to split the phase space into regions such that in each region only a limited subset of momentum configurations can lead to collinear singularities.
Doing that involves the first two points in the itemized list that we presented after Eq. (2.6). The first point is the phase space partitioning; our goal is to split the phase space so that collinear singularities are localized in a clean and physical way. For example, we may want to focus on the collinear emissions off the incoming quark or the collinear emissions off the incoming anti-quark, or on the collinear emission of the gluon g 4 off the quark and the gluon g 5 off the anti-quark etc.
We can do that by introducing a partition of unity and using it to split the phase space. We write 1 = ω 41,51 + ω 42,52 + ω 41,52 + ω 42,51 . (2.21) For the double-collinear partitions {4i, 5j}, i = j, the damping factor ω 4i,5j is engineered in such a way that collinear singularities in ω 4i,5j F LM,qq (1,2,4,5) arise only if momentum p 4 is parallel to p i and/or the momentum p 5 is parallel to p j . Conversely, in the triple-collinear partitions {4i, 5i}, i = 1, 2, the damping factor w 4i,5i is designed in such a way that only the p 4 ||p i , p 5 ||p i and p 4 ||p 5 momentum configurations lead to a singularity. Apart from these conditions, there is significant freedom in choosing the partition functions; we will present a possible choice in the forthcoming sections. 5 Contributions from the double-collinear partitions ω 41,52 , ω 42,51 can be computed right away since the singular limits are easy to establish and no overlapping singularities are present. For example, in case of ω 41,52 , it is sufficient to use the angle between the threemomenta p 4 and p 1 and the angle between the three-momenta p 5 and p 2 as independent variables to describe the collinear singularities in this partition.
The situation is more complex for the triple-collinear partitions, where overlapping singularities are present. The complexity stems from the fact that different hierarchies between ρ 4i , ρ 5i and ρ 45 lead to inequivalent limits in this case. To identify these limits and extract them in a unique way, we further partition the phase space into four sectors. Taking as an example the w 41,51 partition, we introduce four sectors as follows (2.22) The four sectors in the partition w 42,52 are constructed analogously. It is clear that Eq. (2.22) acts in such a way that in each of the four sectors only a small number of singular collinear limits occurs. We then expect that by choosing an appropriate parametrization for each of the four sectors, these singularities can be isolated and extracted. A convenient phase space parametrization for each of the four sectors can be found in Ref. [2]. In each of the four sectors shown in Eq. (2.22), the nested subtraction of these collinear limits can then be performed similar to what we discussed in connection with the soft limits. We sketch how to do this by considering sector where C 51 is an operator that extracts the most singular contribution in the collinear 5||1 limit from the quantity on the left-hand side of Eq. (2.23) and enforces this collinear limit on the damping factor w 41,51 , the reduced matrix element, the momentum-conserving δfunction and the observable O. This amounts to the replacements ρ 51 → 0 and p 1 → p 1 = p 1 (E 1 − E 5 )/E 1 in these quantities. The result reads where P qq is the Altarelli-Parisi splitting function P qq (z) 25) and the "1 " notation in F LM,qq refers to the p 1 → p 1 substitution that we just described. Compared to soft limits, there is an additional subtlety. Indeed, in our construction the angular part of the phase space is non-trivial. To unambiguously define the C 51 operator, we must specify its action on the gluons' phase space [df 4 ][df 5 ]. A convenient choice, adopted already in Ref. [16], is to let C 51 act on it, i.e. to take the ρ 51 → 0 limit of the The right hand side of Eq. (2.23) includes a term with reduced kinematics, which can be dealt with using methods similar to the ones used in NLO computations, and another term that only contains a triple-collinear singularity. The latter occurs whenever 4||5||1, without further hierarchy between ρ 51 , ρ 41 and ρ 45 . To regulate this last singularity, we introduce a triple-collinear operator C C 1 and write where we used C C 1 w 41,51 = lim ρ 51 →0,ρ 41 →0,ρ 45 →0 = 1, which immediately follows from the definition of w 41,51 . Similar to the single-collinear case, the operator C C 1 extracts the most singular behavior from the matrix element in the limit ρ 41 ∼ ρ 51 ∼ ρ 45 → 0 and sets where s 145 = s 45 − s 14 − s 15 and P ggq is a triple-collinear splitting function [25] that depends on the invariants s ij = 2E i E i ρ ij and the momentum fractions Note that in the triple-collinear limit the only effect of the gluon emission on the reduced matrix element and the kinematics of the initial state is the boost where P (z) is the integral of the (soft-regulated) splitting function over the phase space of the unresolved gluons, with the constraint E 4 + E 5 = (1 − z)E 1 and E 4 < E 5 . We note that the second term in Eq. (2.28) is free of all singularities, and can be integrated in four dimensions using standard Monte-Carlo techniques.
Although this discussion is valid for any triple-collinear operator C C 1 that extracts the corresponding triple-collinear singularity from the matrix element squared, we must specify the action of C C 1 on [df 4 ][df 5 ] and on the P ggq function itself to unambiguously define the subtraction framework. In Ref. [16], we let C C 1 act on both [df 4 ][df 5 ] and on the splitting function, i.e. we evaluated all the s ij invariants in Eq. (2.27) and the angular factors in the [df 4 ][df 5 ] phase-space in the triple-collinear limit. While this is a valid option, it is not the only one. In fact, this choice makes the analytic integration over angles of the unresolved partons rather complicated, since it constrains the internal rotational symmetry of the unresolved phase space and does not allow for simple reparametrizations.
To overcome these issues, we now define the operator C C 1 in such a way that it does not act on either [df 4 ][df 5 ] or on P ggq . Rather, C C 1 acts on the momentum-conserving δfunction and on the observable, and extracts the leading triple-collinear singularity from the matrix element according to Eq. (2.27), but it leaves the angular factors in the [df 4 ][df 5 ] phase space and all the s ij invariants in Eq. (2.27) untouched. This modification of the subtraction scheme leads to a simpler integration of the triple-collinear splitting function over the unresolved phase space. Indeed, such a calculation has recently been performed for all relevant triple-collinear splitting functions in Ref. [19].
The results of Ref. [19], combined with integrated double-soft subtraction terms presented earlier in Ref. [18], allow us to promote the fully local subtraction framework or Ref. [16] to a fully analytic scheme. This implies that we can now check the cancellation of all infrared poles analytically and achieve faster and more stable physical predictions by using analytic formulas for all the integrated subtraction terms.
We will present analytic formulas required for the computation of NNLO QCD corrections to the production of color-singlet final states in the remaining parts of this paper. However, before we do that, a general comment is in order. Indeed, as should be clear from the discussion in this section, our framework is highly modular; we believe that this modularity ensures that its generalization beyond color-singlet production will proceed seamlessly. Indeed, the only differences between the color-singlet production and the general case with colored partons in the final state are: 1. compared to color-single production, a generic process has a more complicated color structure and requires double-soft integrals that are functions of relative angles of pairs of hard emittors, rather than pure numbers as in Eq. (2.16). The results relevant for this case have been presented in Ref. [18]; 2. a generic process also involves triple-collinear final state splitting. While the analytic integration of the relevant splitting functions over the unresolved triple-collinear phase space has not been performed for all possible splittings, in Ref. [19] it was shown that techniques used to deal with initial state splittings can be successfully applied to final splittings as well.
It follows that the most general ingredients required for computing NNLO QCD corrections to generic partonic processes at the LHC can be obtained. From that perspective, the analytic formulas presented in this paper provide important building blocks for such a generic computation and give an excellent starting point for its generalization that will be addressed in the future. For now, we will proceed with presenting analytic formulas for all partonic channels that may contribute to the production of color-singlet final states at a hadron colliders.

Quark-initiated color-singlet production
In this section, we consider the production of a color-singlet final state to NNLO QCD accuracy for reactions that are quark-initiated at leading order. We refer to these processes as "Drell-Yan processes", however, we emphasize that the results presented in this section are applicable to any color-singlet production process which is quarkinitiated at LO. Typical examples include pp → Z, W + , γ * , ZZ, W + W − , W Z, W H, ZH and so on. Starting from Eq. (2.1), we find it convenient to group the different partonic channels in three categories We omit the dependence of dσ DY f on the renormalization and factorization scales µ R , µ F and the observable O to shorten the notation. The first term in Eq. (3.2), dσ DY faf b , receives contributions from quark channels and is present at LO. The terms on the second line, dσ DY fag and dσ DY gfa , start contributing at NLO, and the last term dσ DY gg appears for the first time at NNLO. In what follows, we will consider the LO, NLO and NNLO contributions in turn. To simplify the notation, we will omit the superscript DY for the rest of this section.
For the NNLO contribution, we will consider the different channels defined in Eq. (3.2) separately. We also find it convenient to split each of these channels further, according to the highest final state multiplicity that they involve, (cf. Eq. (2.6)) Finally, we separate the above terms into those involving only tree-level matrix elements and those terms involving loop corrections, by writing 6 The term dσ NNLO 1245,faf b receives contributions from processes with NNLO-like kinematics (i.e. with two additional resolved partons in the final state), and corresponds to the fully subtracted real-real contribution. The remaining terms arise from integrated subtraction terms, α s and parton distribution function renormalizations, and real-virtual and purely virtual corrections. The terms dσ NNLO 124,faf b and dσ NNLO 12,faf b only involve tree-level matrix elements squared, while dσ NNLO virt 124 ,faf b and dσ NNLO virt 12 ,faf b also involve finite remainders of virtual amplitudes. It is important to emphasize that all of the different terms in Eq. (3.4) are separately finite, so that we can discuss them separately. In what follows, we will present results for each of these terms.

LO and NLO
We start by discussing the quark channel dσ faf b , with a, b = 0, which is the only channel contributing at leading order. The LO cross section reads NLO QCD corrections to the quark channel can be written as (3.6) 6 We note that certain partonic channels only contain a subset of these terms.
with a, b = 0. TheÔ NLO operator that appears in this formula renders the contribution of the single-gluon emission process finite. It readŝ The other quantity in Eq. (3.6), F fin LV,faf b (1, 2) , refers to the finite remainder of the (UVrenormalized) one-loop virtual correction. Its definition is given in Appendix A. Finally, P (0) qq,R and P qq are related to the LO Altarelli-Parisi splitting function, and γ q = 3C F /2, see Appendix C for explicit formulas.
The qg and gq channels start contributing at NLO. They read and analogously

NNLO: quark channels
In this section we consider the NNLO corrections to dσ faf b , with a, b = 0. This includes the partonic processes q iqj → V + gg, q iqj → V + q kql and q i q j → V + q k q l . Of these partonic processes, the q iqj → V + gg has the most complicated singularity structure; it was discussed in detail in Ref. [16] and reviewed in Section 2. Recall that we introduced an energy ordering E 4 > E 5 (cf. Eq. (2.7)), which is natural since the amplitude is symmetric under the exchange of the two final state gluons. The singularity structure is much simpler for final state quarks, where one could use only two sectors to separate the collinear singularities. Nevertheless, we find it convenient to treat the gluon and quark final states on an equal footing. We therefore need to symmetrize the amplitudes involving the final state quarks explicitly, since they are not symmetric in general; we do this by writing If one wishes to consider the final state gluons and quarks separately, one could do away with the energy ordering and the symmetrization of the quark amplitudes. We emphasize that in this case, the formulas in the forthcoming sections would require modifications.
As mentioned in Section 2, an important part of the subtraction scheme is the separation of the phase space into partitions such that in each partition, only a limited number of kinematic configurations leads to collinear divergences, cf. Eq. (2.21). Throughout this paper, we choose the partition functions to be , w 42,51 = η 41 η 52 η 45 η 45 + η 42 + η 51 , where we have used η ij = ρ ij /2. It is straightforward to check that these functions restrict the collinear singularities as discussed in Section 2, and also that they sum up to one, cf. Eq. (2.21). We now present results for the different terms in Eq. (3.4) that arise in the quark channel.

Terms with NNLO kinematics
This (hard) regularized contribution is the only one that involves the full matrix element In this equation, dc = {(1, 2), (2, 1)} and tc = {1, 2} refer to the double-and triple-collinear partitions, respectively, while the sectors (a)-(d ) are defined by the angular ordering conditions in Eq. (2.22). The operators S 5 , S S, C ij and C C i have been discussed in great detail in Ref. [16], and in Sec. 2. We note that dσ NNLO 1245 faf b is computed numerically in four dimensions. In order to do so, we must provide the explicit parametrization of the phase space for the complete final state which includes two radiated partons and a vector boson (or its decay products). It is clear that there are many different ways to do so. We find it useful to describe the phase space using tree-level variables, i.e. the invariant mass M 2 V and the rapidity Y of the vector boson, but other choices are possible. In addition, we have to choose the energies of the two final state partons and the relative angles between them and the hard emittor 7 as independent variables, in order to extract singularities in the same way as in the computation of the integrated subtraction terms, which are presented in the forthcoming subsections. For this reason, there is less freedom in choosing how to parametrize the momenta of the radiated partons. We have discussed this point in some detail in Appendix B of Ref. [16]. We will not repeat this discussion, instead, our goal here is to provide a guide for a numerical implementation of Eq. (3.13).
We work in the center-of-mass frame of the colliding partons. As the first step, we determine the center-of-mass collision energy squared s. To do so, we parametrize the energies of the radiated partons as 8 where x 1 , x 2 ∈ [0 : 1], and use momentum conservation There is an obvious constraint s < S h , where S h is the center-of-mass energy squared of the colliding hadrons, that we have to impose while generating the events. The choice of angular variables depends on the partition and the sector; for the sake of definiteness, we will discuss the sector "a" of the partition w 41,51 . In this case, the scalar products η ij = (1 − cos θ ij )/2 may be parametrized by the variables x 3 , x 4 , λ ∈ [0 : 1] according to The identity of the "hard emittor" depends on the partition. 8 We remind the reader that we have chosen see [2]. We note that, since 0 ≤ x 4 ≤ 1, the angular ordering η 51 < η 41 /2 is assured.
In addition to the invariant mass, we also fix the rapidity of the vector boson in the laboratory frame, Y . This allows us to determine fractions of hadron energies carried by the colliding partons, ξ 1,2 . We find where We require that 0 < ξ 1,2 < 1 and that both the numerator and the denominator in the argument of the logarithm in Eq. (3.19) are positive definite.
We are now in position to write down the four-momenta of the QCD partons in an event. We do so in the partonic center-of-mass frame. The knowledge of ξ 1,2 then allows us to boost momenta to the laboratory frame where all kinematic constraints are defined. The momenta read where cos θ ij = 1 − ρ ij and [2] sin The four-momentum of the vector boson is obtained by momentum conservation p V = p 1 + p 2 − p 4 − p 5 . If needed, further details of the colorless final state can be described.
For example, in case of V → l + l − , the phase space for leptonic decays is generated in the V -rest frame and the lepton momenta are boosted back into the partonic center-of-mass frame using the known p V and M 2 V . For the chosen partition and sector, the phase space weight reads where w LO is the weight for the Born f a f b → V process, which depends in general on a set of variables {y i } that describe the V final state. The contribution of the generated hard event to the phase-space integral is then The matrix element squared can be calculated either in the center-of-mass frame or in the laboratory frame since all required boosts are defined at this point. The contribution that we just described corresponds to the product of identity operators in Eq. (3.13); below we discuss how the subtraction terms in Eq. (3.13) are to be calculated. To this end, we consider first the class of terms in Eq. (3.13) where the double-soft operator S S appears. We will start with the term S SF LM and describe the weight of the counter-event produced by this term. To compute the weight, we take the limit x 1 → 0 everywhere; this corresponds to E 4,5 → 0 at E 5 /E 4 held fixed. We obtain  Suppose we consider terms in Eq. (3.13) where, in addition to the double-soft operator S S, some other operator acts on F LM,qq . In this case, we should just set the relevant variables(s) to zero in Eq. (3.26) and, if necessary, change the way the four-momenta are generated. For example, consider a term C 51 S SF LM,qq (1,2,4,5). For this sector, the operator C 51 implies that x 4 should be set to zero everywhere after the leading 1/x 4 asymptotic is extracted. The four-momenta are then unchanged, except for p 5 which becomes Computing the 5||1 limit of the double-soft eikonal function, we arrive at the contribution from the C 51 S SF LM,qq (1, 2, 4, 5) kinematic configuration where We note that, according to these equations, s 15 = 2(p 1 · p 5 ). This is so because the 1/s 15 term describes the leading x 4 → 0 singularity, so x 4 = 0 is kept there and set to zero everywhere else. As the last example, consider the triple-collinear limit, which corresponds to x 3 → 0, cf. Eqs. (3.16, 3.20). In this case, the partonic center of mass collision energy squared is Similar to the case of hard event, we require 0 < s C C < S h . The four-momenta of partons to be used in the matrix element and the momentum-conserving δ-function are 31) and the vector boson four-momentum is As we emphasized in Sec. 2, we define C C i in such a way that it does not act on the phase-space. The weight of the counter-event is then identical to the one for the hard process and the parameter y relevant for calculating fractions of hadron energies carried by the incoming partons reads The momentum fractions themselves are then given by Combining the different ingredients, we derive the weight of the triple-collinear counterevent where p 1 = (1 − x 1 (1 + x 2 ))p 1 . The arguments of the triple-collinear splitting function P ggq are then computed as and We stress that the above scalar products in the splitting function are evaluated with x 3 = 0, i.e. not in the triple-collinear limit.
The remaining contributions to the fully-subtracted cross section dσ NNLO 1245,faf b are computed along the same lines. The important thing is that we always take the leading singularity in the relevant variables and employ the limiting behavior of amplitudes squared to calculate weights of the subtraction terms. We also make sure that the subtraction counter-terms that make the hard matrix element finite are identical to the subtraction terms that have been analytically integrated.

Tree-level terms with NLO kinematics
In this section, we consider the term with NLO kinematics and tree-level matrix elements, dσ NNLO 124,faf b . The general structure of this contribution is where the various splitting functions are defined in Appendix C. Although we only need ∆ q in Eq.(3.38), it appears to be convenient to introduce a more general object ∆ i∈{q,g} . It is defined as In the case i = q, we find C q = C F , γ q = 3C F /2 and X q = 3C A /2, cf. Appendix. B. In addition Note that to obtain Eq. (3.39), we used η 41 + η 42 = 1. In Eq. (3.39), F µν LM,faf b is analogous to F LM,faf b but with the polarization vector for gluon 4 removed (3.41) F µν LM is contracted with r µ r ν where r µ is a unit vector that spans the two-dimensional space orthogonal to p 4 , see Ref. [16] for further details. If p 4 is parametrized as in Eq. (3.44)

Tree-level terms with LO kinematics
We now turn to the contribution involving terms with LO kinematics and tree-level matrix elements, dσ NNLO 12,faf b . Accounting for the boost of the initial state along the collision axis, it can naturally be split into

45)
a, b = 0. We now consider each of these terms separately.

Terms involving
(3.46) Once again, to illustrate this equation we consider the case of Z production. Here, this contribution is only relevant for the (f a , f b ) = (q,q) channel, where it reads (3.47) 2. Terms involving F LM (z · 1, 2) and F LM (1, z · 2): This term has a non-trivial flavor structure. To simplify it, we employ the notation used to describe the NNLO QCD contributions to the Altarelli-Parisi splitting functions, and write the functions T in terms of nonsinglet, singlet, and vector functions (3.49) Similar to the Altarelli-Parisi splitting functions, we have T S qq = T S qq through NNLO, and we will always use the latter in what follows. Again, we consider the example of Z production. For the (f a , f b ) = (qq) channel, Eq. (3.48) becomes (3.51) The transition function T NS qq is explicitly shown in Appendix D. All other T ij functions are presented in an ancillary file.

Terms involving virtual corrections
Finally, we consider the two terms in Eq. and dσ NNLO virt 12 ,faf b . The former corresponds to the real-virtual corrections, which have NLO kinematics. As such, they have singularities that appear when the radiated parton becomes unresolved. These singularities can be subtracted as at NLO, so that this term reads where F fin LV,faf b (1, 2, 4) is a finite remainder of the one-loop amplitude, see Appendix B. The other term corresponds to virtual contributions with LO kinematics. It reads where γ q is defined after Eq. (3.39) and F fin LVV , F fin LV 2 and F fin LV are defined in Appendix A.

NNLO: quark-gluon channels
In this section, we describe the NNLO contributions to the qg channel, see Eq. (3.2). Similar results hold for the gq channel. In principle, this channel could be treated in the same fashion as the quark channels discussed in the previous section. However, its singularity structure is much simpler, and so we need to consider a smaller number of limits. Indeed, no double-soft singularities are present in this case. Because of this, we find it convenient not to order the energies of partons 4 and 5. We write and parametrize E 4,5 = x 1,2 E max . However, the structure of the collinear singularities is similar to that discussed in Sec. 3.2, so we use the same angular parametrization and partitioning as defined there. There is another important difference compared to the qq channel discussed in Sec. 3.2, namely that in the qg channel spin correlations appear in the collinear emissions off the incoming gluon. We postpone their discussion to Sec. 4.1, where we consider the most general case of spin correlations. Apart from this, the structure of the result is very similar to the one discussed previously, so we limit ourselves to reporting the relevant formulas.
We write (3.59) We consider the case with a = 0, and discuss each term separately. Note that in this channel there are no terms proportional to F LM (1, 2) and to F LM (z · 1, 2) since the process qg → V at leading order is impossible. For the other terms, we obtain the following results.
1. Tree-level terms with NNLO kinematics: which is analogous to Eq. (3.38) for the quark channels. The splitting functions in Eq. (3.61) are defined in Appendix C, and ∆ (qg) is given by  3. Tree-level terms with LO kinematics involving F LM (z · 1,z · 2): analogous to Eq. (3.56). The finite remainder F fin LV is defined in Appendix A.
Results for the gq channel can be obtained from the above formulas in a straightforward manner, by replacing labels 1 ↔ 2.

NNLO: gluon-gluon channel
In this section, we describe the gg channel, see Eq. (3.2). This case is particularly simple, since no soft or triple-collinear singularities are present. As the result, only double-collinear configurations need to be considered. As a consequence, it is not necessary to partition the phase space in any way and the singularity structure can be dealt with as in NLO computations described earlier. We parametrize the energies and angles of the emitted partons as (3.67) Singularities only appear when x 3,4 are equal to either zero or one. Because of the simplicity of the singularity structure, we treat the two cases at once. To deal with them, we use a similar but simpler strategy to the one discussed in the preceding sections. We write  At variance to the cases discussed above, virtual corrections do not contribute in this channel. The individual terms read as follows.

Gluon-initiated color-singlet production
In this section, we consider the production of a color-singlet final state H in gluon fusion through NNLO QCD. We refer to this process as "Higgs production", although, similarly to the "Drell-Yan process" discussed in the previous section, these results are applicable to the production of any color-singlet final state which proceeds through gluon fusion at LO. The procedure of extracting the infrared divergences is identical to what has already been discussed in the case of the qq annihilation and we do not repeat it. However, in gluon-initiated processes, initial state radiation leads to spin correlations that we did not discuss up to now. In the next section we show how to deal with this complication.

Spin correlations
We have discussed spin correlations relevant to the computation of NNLO QCD corrections to the qq → V process in Ref. [16]. In that case, the spin correlations appeared because of the splitting of a virtual gluon to two final state partons, g * → f 4,5 . For the gg → H process, the situation is different in that spin correlations also appear in the initial state radiation, including its triple-collinear limit. In this section, we discuss this case. We will begin with the discussion of NLO QCD corrections to gg → H. The computation proceeds exactly as for the qq → V process, cf. Sec. 3.1. However, when the collinear operator acts on the matrix element squared, we find where the splitting function reads The transverse momentum k ⊥ is defined using the Sudakov decomposition where k ⊥ · p 1 = k ⊥ · p 2 = 0. The metric tensor of the transverse space g ⊥,d−2 satisfies g µν ⊥,d−2 p 1,ν = g µν ⊥,d−2 p 2,ν = 0 and g µν ⊥,d−2 k ⊥,ν = k µ ⊥ . We write the four-momenta of the QCD partons as for the qq → V process. We introduce d−dimensional vectors t µ = (1, 0) and e µ 3 = (0, 0, 0, 1, 0), and an additional vector a µ that is orthogonal to both t and e 3 and is normalized a 2 = −1. We write the four-momenta in terms of t, e 3 and a and obtain p 1,2 = E 1 (t ± e 3 ), p 4 = E 4 (t + cos θ 41 e + sin θ 41 a). (4.4) By comparing the two parametrizations of the vector p 4 , we find Since the transverse components of the gluon four-momentum decouple from the hard matrix element in the collinear limit, the vector a µ only appears in the splitting function. We can then integrate the splitting function over the directions of the vector a µ using We find dΩ is the spin-averaged splitting function. It follows that averaging over the directions of the transverse components of the gluon momentum naturally appears in our construction at NLO; as the consequence, the rest of the NLO QCD calculation is identical to the qq case. Before discussing spin correlations in the computation of NNLO QCD corrections, we note that, in the particular case of Higgs boson production, spin correlations are actually not needed. Indeed, the spin-correlated matrix element for gg → H is proportional to g µν ⊥,d−2 . As the result, the spin-averaged splitting function P gg,µν g µν ⊥,d−2 naturally appears in the calculation. We emphasize, however, that this feature is particular to the process gg → H, so that understanding spin correlations is necessary in a more general context.
We then consider the generic NNLO case. Here, the situation is more complex since we have to consider the momenta of the two gluons becoming collinear to the direction of the incoming partons. In the double-collinear partitions the situation is identical to the NLO case since the averaging over the transverse spaces of the two gluons is performed independently. The triple-collinear partitions require some discussion. We consider the case when collinear singularities arise because of the emissions of gluons g 4,5 off the gluon g 1 . We parametrize the four-momenta of the final-state gluons as [16] p 4 = E 4 (t + cos θ 4 e 3 + sin θ 4 a) , p 5 = E 5 (t + cos θ 5 e 3 + sin θ 5 (cos ϕ 45 a + sin ϕ 45 b)) , (4.9) where the vectors t, a, e 3 have already been defined in the discussion after Eq.(4.3) and the vector b satisfies t · b = e 3 · b = a · e 3 = 0, as well as b 2 = −1.
We begin with the double-collinear limits that develop spin correlations. There are three possibilities: g 4 is collinear to g 5 , g 5 is collinear to g 1 and g 4 is collinear to g 1 . The first case is identical to the qq → V process; it was discussed in Ref. [16] and we do not repeat this discussion here. The second case, p 5 ||p 1 , relevant for sector (a), is discussed below. After that we comment on the third case, relevant for sector (c).
Starting from the angular part of the phase space for sector (a) and considering the limit x 4 → 0, corresponding to θ 51 → 0, we find lim see Ref. [16]. Also, in this limit This follows immediately from the definition of sin ϕ and by inverting the definition of λ in terms of cos ϕ. We then rewrite We also have (4.14) These identities imply that in the x 4 → 0 limit Hence, in case of double-collinear limits with respect to the incoming partons, integration over the transverse directions of the collinear gluons always leads to spin-averaged splitting functions. This implies that subsequent computational steps are conceptually identical to those of the qq → V process. Finally, we note that the above discussion can be repeated verbatim also in case p 4 ||p 1 , relevant for sector (c), if instead of Eq. (4.9) we use p 4 = E 4 (t + cos θ 4 e 3 + sin θ 4 (cos ϕ 45 a + sin ϕ 45 b)) , p 5 = E 5 (t + cos θ 5 e 3 + sin θ 5 a) . It remains to discuss the triple-collinear limit that corresponds to the splitting g 1 → g 4 + g 5 + g * . This splitting is described by the P µν ggg splitting function that contains spin correlations, see e.g. Ref. [25]. This splitting function is a symmetric rank-two tensor constructed from g µν ⊥,d−2 , and the vectors k µ 4(5),⊥ . These vectors read, c.f. Eq. (4.9), k µ 4,⊥ = E 4 sin θ 4 a µ , k µ 5,⊥ = E 5 sin θ 5 (a µ cos ϕ 45 + b µ sin ϕ 45 ) . In the triple-collinear limit, we need to integrate over the directions of the vectors a and b. Note that the integration over the angle ϕ 45 is non-trivial since 2(p 4 · p 5 ) = s 45 depends on it. To describe the integration of different tensor structures over the directions of a and b, we introduce the notation We find the following results for the four tensor structures that contribute to P µν (4.19) We write the spin-correlated splitting function as 20) and observe that Eq. (4.19) leads to The same result is obtained upon replacing the spin-correlated splitting function with its spin-averaged version P µν ggg → P ggg,αβ g αβ Once this is done, the triple-collinear splittings in the gg → H process can be treated in exactly the same way as in the Drell-Yan qq → V case.

Definition of partonic channels
After discussinng spin correlations, we proceed with setting up the NNLO QCD calculation for color-singlet production in gluon fusion. Starting from Eq. (2.1), we find it convenient to write the cross section for the generic process as (4.24) The first term is the gg channel which is the only partonic channel contributing at LO. The terms on the second and third line correspond to the quark-gluon channels and quarkantiquark channel respectively and enter at NLO. The q i q j channel, where q i and q j can be both identical or different (anti)quarks first appears at NNLO. We will discuss each of these channels separately in the following subsections. For simplicity, we will omit the "H" superscript from now on. We express our results in terms of fully renormalized amplitudes for the pp → H + X process, where H is a generic color-singlet state. In the case of Higgs production in the m t → ∞ approximation, this implies that our results include both the divergent and the finite renormalization of the Hgg Wilson coefficient, see App. A for more details.
with γ g = β 0 = 11C A /6 − n f /3, and F fin LV,gg and the various splitting functions defined in Appendix A and Appendix C, respectively. NLO corrections to the gq and qg channel read and respectively. Finally, the qq channel starts contributing at NLO but it is finite finite at this order and can simply be written as qq (1, 2, 4) .

NNLO: gluon channel
This channel has the same singularity structure as the Drell-Yan quark channels, cf. Sec. 3.2, and we use the same phase-space parametrization and partitioning described there. This means that the structure of the result is identical to what was discussed in Sec. 3 We now list all these terms separately.
2. Tree-level terms with NLO kinematics: which is analogous to Eq. (3.38). The splitting functions in Eq. (4.33) are defined in Appendix C, and ∆ g is given in Eq. (3.39) with C g = C A and γ g = X g = β 0 .

NNLO: quark-gluon channels
The structure of this channel is analogous to the qg channel for the Drell-Yan process, discussed in Sec. 3.3. We don't repeat the discussion here, and limit ourselves to presenting final results. To closely follow Sec. 3 We display the individual contributions below.

NNLO: qq channel
The singularity structure of this channel can be organized as follows. There are purely collinear singularities, coming from configurations where the q and theq emit two t-channel gluons that produce a Higgs. These have the same structure as those appearing in the gg channel for the Drell-Yan process. Apart from these, there are other singular contributions, that don't have a Drell-Yan equivalent. They stem from extra gluon emission from the s-channel qq → H + g process. These are of NLO origin, so they don't pose any particular challenge.
Because of its simple singularity structure, we don't discuss this channel in detail. For completeness, we present final formulas that are obtained using the same setup that we employed for the gg channel. We write   qq (1, 2, 4) .
3. Tree-level terms with LO kinematics involving F LM (z · 1,z · 2): The splitting functions used in these equations are defined in App. C, ∆ q is given in Eq. (3.39) with C i = C F , γ q = 3C F /2, X q = 3C A /2, and the F fin LV finite remainder is defined in App. B. Note that, contrary to all the cases discussed so far, F fin LV,qq does not require any additional regularization.

NNLO: quark channels
The singularity structure of this channel is the same as the gg channel for the Drell-Yan process. Because of this, we use the same parametrization described in Sec. 3.4. We write dσ NNLO (4.54) Repeating the same steps discussed in Sec. 3.4 we obtain the following results.

Validation of results
In this section, we describe the numerical checks that have been used to validate results described in the preceding sections. We use the processes pp → Z and pp → H as test cases, since for both of these processes NNLO QCD corrections to the inclusive cross sections are known analytically [20,21]. This allows us to perform a high-precision check of our formulas. We begin by describing our setup. We consider proton-proton collisions with 13 TeV center-of-mass energy. We use m Z = 91.1876 GeV, m H = 125 GeV and m t = 173.2 GeV for the Z, the Higgs and the top quark masses, respectively. We derive the weak coupling constant from g 2 W = 4 √ 2m 2 W G F , with m W = 80.398 GeV and G F = 1.16639×10 −5 GeV −2 . The weak mixing angle is computed from sin 2 θ W = 1 − m 2 W /m 2 Z . We will consider both on-shell pp → Z production and pp → e + e − production. In the latter case, we include both the Z and the γ * contributions, and we use Γ Z = 2.4952 GeV.
For the Higgs case, we consider pp → H production in the m t → ∞ approximation and describe Higgs coupling to gluons using the effective Lagrangian L I = −λ Hgg HG (a) µν G µν,(a) , see App. A for details. The coupling λ Hgg depends on the Higgs vacuum expectation value v. For our results, we use v 2 = (G F √ 2) −1 . All computations are done using the NNPDF3.0 parton distribution set [26], with 5 active flavors. We employ LO/NLO/NNLO sets for LO/NLO/NNLO predictions, respectively. We use the value of the strong coupling and its evolution provided by the PDFs sets, with α s (m Z ) = 0.118 at (N)NLO, and α s (m Z ) = 0.130 at LO.
We first consider fully inclusive on-shell Z production. We compare results obtained within our framework to the analytic results of Ref. [20] that we implemented in HOPPET [27]. We study each partonic channel individually, and additionally split some of them into contributions with different color factors in order to validate all the different singularity structures independently. We show results for a single fixed factorization and renormalization Channel Color structures Numerical result (nb) Analytic result (nb) scale µ = 2m Z , although we have performed the same check for different scales as well. We note that these inputs are not chosen for their phenomenological relevance, but rather to provide stringent checks on our results. We summarize our findings in Table 1. It shows that our framework allows for extremely high precision results, with numerical errors at the per mille level or better 9 . These results are always fully compatible with the analytic ones within the numerical uncertainties. We remind the reader that these numbers refer to the NNLO coefficients, which implies absolute precision on the physical cross section. Analogous results for the case of Higgs production for equal renormalization and factorization scales µ R = µ F = m H /2 are shown in Tab. 2. Again, we find it convenient to perform numerical checks for different partonic channels independently. Also in this case, our numerical results have tiny uncertainties and are in perfect agreement with the analytic values obtained from Ref. [21].
Having fully validated our results, we now briefly investigate the performance of the framework when applied to the computation of physically relevant predictions. Specifically, we explore the computational effort required to obtain predictions for physical quantities at the per mille level. We start by considering inclusive Higgs production, at the 13 TeV LHC. For this study, we set µ R = µ F = m H . Running for less than an hour on a single core of a standard laptop, we obtain  We now move to fiducial cross sections. We consider pp → Z/γ * → e − e + production in the fiducial volume defined by symmetric lepton cuts studied in Ref. [28]. Specifically, we require that the transverse momentum and rapidity of each lepton satisfy p T, > 25 GeV |η | < 2.47, (5.2) and that the invariant mass of the lepton pair is in a window 66 GeV < m e − e + < 116 GeV.
In this case, we use µ R = µ F = m Z . Running on a single core of a standard laptop for about an hour, we obtain We note that in this case the error is at the few per mille level. We compared the NNLO K-factor against the benchmark result presented in Ref. [28], and found agreement within the numerical precision. As a final comment, we note that although the processes studied here are very simple, which makes it difficult to predict how the framework will perform for more complicated ones, these results are very encouraging.

Conclusions
In this paper, we presented compact analytic formulas that describe fully-differential production of color-singlet final states in hadron collisions. We studied final states that, at leading order, can be produced either in qq or in gg annihilation.
Our calculation employs the nested soft-collinear subtraction scheme that we developed earlier in Ref. [16]. However, compared to its original formulation, we found it useful to modify it to allow for a simpler analytic integration of the triple-collinear limits. We explained the required changes in Section 2.
We validated our results by using them to numerically compute NNLO QCD contributions to total cross sections of Z and H production in proton collisions and comparing them with results based on the convolution of known analytic results for partonic cross sections relevant for these two processes with parton distribution functions. In both cases, we found an agreement at a better-than-per-mille level for NNLO coefficient functions. We also showed that our computation can deal with fiducial cuts quite efficiently.
As far as we know, the results presented in this paper provide the first implementation of a fully local and fully analytic NNLO QCD subtraction scheme. Although -as we already emphasized in the introduction -it is difficult to say to what extent these nice features of the subtraction scheme will help with its efficiency, we hope that they will be helpful in that respect. However, we also believe that, regardless of the performance issues, physical clarity and overall transparency of the obtained formulas gives us hope that analytic, local NNLO QCD subtractions for arbitrary complex hadron collider processes are within reach.

A Purely virtual contribution: definitions
We consider the UV-renormalized amplitude for the process p 1 + p 2 → V A(p 1 , p 2 ; {p V }) = A 0 (1, 2) + α s (µ) 2π A 1 (1, 2) + α s (µ) 2π T R n f . To express the virtual contribution to the cross-section, it is also useful to define The NLO virtual correction can then be written as with F fin LVV and F fin LV 2 finite and proportional to 2Re A 0 A * 2,fin and |A 1,fin | 2 , respectively.

A.1 Finite remainder: Drell-Yan
In this section, we report the finite remainders for the Drell-Yan process, see e.g. [30]. We obtain with Q 2 = p 2 V and α s = α s (Q). Results for generic µ can easily be obtained from renormalization group evolution (RGE) arguments.

A.2 Finite remainder: Higgs
In this section, we report the finite remainders for the Higgs process. More precisely, we consider a theory where the Higgs is coupled directly to gluons through the effective interaction Lagrangian L I = −λ Hgg HG (a) µν G µν,(a) , (A.11) where the (bare) Hgg coupling is given by In this formula, α s = α s (µ) is the renormalized coupling in a theory with 5 light flavors, v is the Higgs v.e.v. and the divergent (Z eff (α s )) and finite (C(α s )) parts of the Wilson coefficient renormalization are given in the MS scheme by (A.13) see e.g. [31]. In Eq. (A.13), m t is the top-quark mass, β 0 has been defined in Eq. (A.4) and Combining the result for the Hgg form factor in e.g. [30] with the finite part of the Wilson coefficient renormalization, we obtain for the Higgs finite remainders with Q 2 = p 2 H and α s = α s (Q). Results for generic µ can easily be obtained from RGE arguments.

B Real-virtual contribution: definitions
We consider the one-loop amplitude for the process with F fin LV,ij (1, 2, 4) finite. The explicit form of I 124 depends on the color structure of the process. Since the process Eq. (B.1) must involve either a gluon and a qq pair or 3 gluons, we can classify the most general case according to the position of a gluon and write Here s ij = 2E i E j ρ ij , with E i,j > 0.