Local Analytic Sector Subtraction at NNLO

We present a new method for the local subtraction of infrared divergences at next-to-next-to-leading order (NNLO) in QCD, for generic infrared-safe observables. Our method attempts to conjugate the minimal local counterterm structure arising from a sector partition of the radiation phase space with the simplifications following from analytic integration of the counterterms. In this first implementation, the method applies to final-state massless particles. We show how our method compactly organises infrared subtraction at NLO, we deduce in detail the general structure of the subtraction terms at NNLO, and we provide a proof of principle with a complete application to a simple process at NNLO.


Introduction
The increasing precision of experimental measurements at the Large Hadron Collider (LHC), together with the complexity of the final states currently probed in hadronic collisions, constitute a severe challenge for theoretical calculations. This challenge has driven the development of a number of novel techniques, for precision calculations of scattering amplitudes to high orders, for the study of final-state hadronic jets, and for the accurate determination of parton distribution functions (see, for example, Ref. [1] for a review of recent developments). In particular, a consequence of the current and expected precision of experimental data is the fact that the next-to-next-to-leading perturbative order (NNLO) in QCD is rapidly becoming the required accuracy standard for fixedorder predictions at LHC. A crucial ingredient for the calculation of differential distributions to this accuracy is the treatment of infrared singularities, which arise both in virtual corrections to the relevant scattering amplitudes, and from the phase-space integration of unresolved real radiation. In principle, the problem is well understood. Infrared singularities (soft and collinear) arise in virtual corrections as poles in dimensional regularisation, and all such poles are known to factorise from scattering amplitudes in terms of universal functions, which admit general definitions in terms of gauge-invariant matrix elements [2][3][4][5][6][7][8][9][10][11]. These functions are in turn determined by a small set of anomalous dimensions which, in the massless case, are fully known up to three loops [12,13]. General theorems then ensure that, when considering infrared-safe cross sections, virtual infrared poles must either cancel, when combined with singularities arising from the phase-space integration of final-state unresolved radiation [14][15][16][17], or be factored into the definition of parton distribution functions, in the case of collinear initial-state radiation [18]. Real-radiation matrix elements have also been shown to factorise in soft and collinear limits, and the corresponding splitting kernels are fully known at order α 2 S [19][20][21][22][23][24], with partial information available at α 3 S as well [25][26][27][28][29]. Even with this detailed knowledge of the relevant theoretical ingredients, the practical problem of constructing efficient and general algorithms for handling infrared singularities for generic infrared-safe observables beyond next-to-leading order (NLO) proves to be highly non-trivial. The origin of the difficulty lies in the fact that typical hadron-collider observables have a complicated phase-space structure, nearly always involving jet-reconstruction algorithms as well as complex kinematic cuts; furthermore, real-radiation matrix elements become increasingly intricate, and they cannot be analytically integrated in d dimensions. Integration over unresolved radiation must therefore be performed numerically in d = 4, and all infrared singularities must be cancelled before this stage of the calculation is reached. This cancellation involves a careful use of approximations to the real-radiation matrix elements in the singular regions, and requires a remapping of the real-radiation phase space to match the Born-level configurations.
At NLO, the first fully differential results for jet cross sections were obtained [30,31] by isolating singular phase-space regions and treating them separately, performing the pole cancellation by integrating approximate matrix elements within those regions (a procedure usually described as 'slicing'). Subsequently, two general algorithms were developed, the FKS [32] and CS [33] subtraction methods, based on the idea of introducing local counterterms for all singular regions of phase space, and then integrating them exactly in order to achieve the cancellation of poles without need of slicing parameters (which is usually described as 'subtraction' in a strict sense). These algorithms are currently implemented in full generality in fast and efficient NLO generators [34][35][36][37][38][39][40][41][42], so that the 'subtraction problem' can be considered solved to this accuracy.
At NNLO, numerical and conceptual challenges related to the proliferation of overlapping singular regions become much more significant. This has led to the development of several different methods, which have been successfully applied to a number of simple collider processes. NNLO differential distributions for hadronic final states in electron-positron annihilation were first computed in [43,44], while among the first hadronic processes involving coloured final-state particles to be studied differentially at NNLO there were the production of top-antitop quark pairs, achieved in [45,46] within the Stripper framework [47], and the associated production of a Higgs boson and a jet, achieved with the N-Jettiness slicing technique [48][49][50][51]. A number of hadronic processes with up to two final-state coloured particles at Born level have since been studied at the differential level with various approaches, including q T slicing [52][53][54][55], and Antenna subtraction [56][57][58].
There are several reasons to surmise that existing methods for NNLO subtraction can be generalised and improved: on the one hand, current applications have been computationally very demanding, either in terms of the analytic calculations involved, or because of the large-scale numerical effort required; on the other hand, it is clear that precise NNLO predictions will soon be needed for more complicated processes, such as the production of more than two jets, and it will similarly be useful to compute simple processes at the next order in perturbation theory, N 3 LO.
The need for improved and efficient subtraction algorithms is in fact leading to the development of other methods, or the refinement of existing ones: examples include the CoLoRFulNNLO framework [59][60][61], currently applied to processes with electroweak initial states, the Projection to Born method [62], and the technique of Nested Soft-Collinear subtractions [63,64]. New ideas are also being introduced [65, 66], and the first limited applications to differential N 3 LO processes have appeared [67][68][69].
In this paper, and in a companion paper devoted to the underlying factorisation framework [70], we present a new approach to the subtraction problem beyond NLO, which attempts to re-examine the fundamental building blocks of the subtraction procedure, take advantage of all available information, and build a minimal structure which will hopefully help to streamline and simplify future applications. The ideal subtraction algorithm, in our view, should aim to achieve the following goals: complete generality across infrared-safe observables; exact locality of infrared counterterms in the radiative phase space; independence from 'slicing' parameters identifying singular regions of phase space; maximal usage of analytic information in the construction and integration of the counterterms; and, of course, computational efficiency of the numerical implementation. These are, clearly, overarching goals, and in this paper we present the first basic tools that we hope to use in future more general implementations. In particular, we focus for the moment on the case of massless final-state coloured particles.
In order to achieve the desired simplicity, we attempt to take maximal advantage of the available freedom in the definition of the local infrared counterterms, exploiting and extending ideas that have been successfully implemented at NLO. In particular, a key element of our approach is the partition of phase space in sectors, each of which is constrained to contain a minimal subset of soft and collinear singularities, in the spirit of FKS subtraction [32]. A crucial ingredient is then the choice of 'sector functions' used to build the desired partition: these functions must obey a set of sum rules in order to simplify the analytic integration of counterterms when sectors are appropriately recombined. A second crucial ingredient is the availability of a flexible family of parametrisations of momenta within each sector, allowing for simple mappings to Born configurations in different unresolved regions. Finally, it is necessary to take maximal advantage of the simple structure of factorised kernels in multiple singular limits, which follows in general from the factorised structure of scattering amplitudes: a detailed analysis of this structure will be presented in [70].
With this general strategy in mind, we begin in Section 2 by revisiting the NLO subtraction problem. We define sector functions satisfying our requirements, we introduce local counterterms and appropriate parametrisations, and we integrate the counterterms on the unresolved phase space. Effectively, Section 2 constructs a complete NLO subtraction algorithm for massless final states, which stands out for the simplicity of the required integrations. In Section 3 we attack the NNLO problem, displaying the general structure of subtractions in our approach, defining sector functions, and constructing all local counterterms relevant for massless final states. We then perform the relevant integrations for a specific subset of singularities, and, in Section 4, we use the results to complete a proof-of-concept calculation of NNLO subtraction for the leptonic production of two quark pairs. We conclude in Section 5, outlining the status of our method and the forthcoming steps needed to construct a competitive algorithm. Four appendices contain a number of technical details.
2 Local analytic sector subtraction at NLO

Generalities
We restrict our analysis to reactions featuring only massless particles, with n partons appearing in the final state at Born level. We assume a colour-singlet initial state, and we allow for coloured and colourless particles in the final state, the latter not affecting our arguments. Scattering amplitudes involving n final-state partons with momenta k i , i = 1, . . . , n, with k 2 i = 0, are expanded in perturbation theory as n describing the Born process. Correspondingly, differential cross sections with respect to any infrared-safe 1 observable X are schematically written as where, up to NLO, where the virtual correction has been renormalised in the MS scheme. Furthermore, δ i (X) ≡ δ(X − X i ), with X i representing the observable under consideration, computed with i-body kinematics.
In dimensional regularisation, in d = 4 − 2 space-time dimensions, the virtual contribution features up to double IR poles in , while the real contribution, finite in d = 4, is characterised by up to two overlapping singular limits of soft and collinear nature in the radiation phase space. The phase-space integration of such singularities in d dimensions results in explicit poles in , which cancel those of virtual origin if X is infrared safe, ensuring the finiteness of the cross section [15,16].
The NLO-subtraction procedure avoids analytic integration of the full real-radiation amplitudes by adding and subtracting to Eq. (2.4) a counterterm The combination d Φ n+1 K must reproduce all singular limits of the real-radiation contribution dΦ n+1 R, and must be sufficiently simple to be analytically integrated in d dimensions. Note that we allow for the possibility of simplifying the phase-space measure dΦ n+1 to d Φ n+1 in the counterterm, under the assumption that the two coincide in all singular limits. Defining now the (single) radiation phase space as d Φ rad = d Φ n+1 /dΦ n , we may introduce the integrated counterterm and rewrite identically the NLO cross section in Eq. (2.4) in subtracted form as where the first and the second lines are separately finite in d = 4 and do not present any phase-space singularities, allowing an efficient numerical integration.

Sector functions
Our first step in setting up the subtraction formalism at NLO is to introduce a partition of the real-radiation phase space by means of sector functions W ij , inspired by the FKS method [32], and satisfying the following properties where π(ij) = {ij, ji}. S i and C ij are projection operators on the limits in which parton i becomes soft (i.e. all components of its four-momentum approach zero), and partons i and j become collinear (i.e. their relative transverse momentum approaches zero), respectively: the action of these operators on matrix elements and sector functions will be described in detail below. Eq. (2.9) is a normalisation condition that recognises the W ij functions as a unitary partition of phase space. Eq. (2.10) and Eq. (2.11) express the fact that a given sector function W ij selects only one soft and one collinear singular configurations, S i and C ij , respectively, among all those present in the real-radiation matrix element. The sum rules in Eq. (2.12) imply that, upon summing over all combinations of indices associated to sectors that survive in a given soft or collinear limit, the corresponding sector functions reduce to unity. This fact proves crucial for the analytic integration of the subtraction counterterms, as is well known in the FKS method, and as we will further discuss in the following; analytic counterterm integration in turn makes it possible to show in closed form the correctness of the singularity structure of the subtraction terms. There is ample freedom in the choice of sector functions, the only requirement being that they satisfy the relations (2.9) to (2.12). In order to provide an explicit definition of W ij , let us introduce some notation: let s be the squared centre-of-mass energy, q µ = ( √ s, 0 ) the centre-of-mass fourmomentum, and k µ i (i = 1, . . . , n + 1) the n + 1 final-state momenta of the radiative amplitude. We set We now define NLO sector functions as (see also [39]) (2.14) The double sum in Eq. (2.14) runs over all massless final-state partons, including those that are not associated with singular limits. This choice is made in order to ease NNLO extensions, as detailed below. With the definition in Eq. (2.14), it is easy to verify that all properties in Eqs. (2.9) to (2.12) are satisfied, and in particular one finds that from which the desired properties follow.

Definition of local counterterms
As discussed above, properties (2.10) and (2.11) ensure that, in a given sector ij, only the S i and the C ij limits (as well as their product) act non-trivially. A candidate local counterterm K ij for the real matrix element R in this sector can thus be built collecting all terms in the product R W ij that are singular in such soft and collinear limits, and taking care of correcting for the double counting of the soft-collinear region. We define therefore Here and in the following, projection operators are understood to act on all quantities to their right, unless explicitly separated by parentheses: for instance in the expression (S i A) B the soft limit is meant to act only on A, and not on B. In Eq. (2.16), the term featuring the composite operator S i C ij removes the soft-collinear singularity, which is double-counted in the sum S i + C ij ; the order in which the projectors act is arbitrary, as they commute (see Appendix A). As will be detailed in Section 2.4, and can be deduced from the sum rules in Eqs. (2.12), the content of each square bracket in Eq. (2.17) is equal to 1 upon summation over sectors, a crucial property for counterterm integration.
Our candidate counterterm K ij is structurally similar to, and as simple as, the FKS counterterm for sector ij, however it has the advantage of being defined without any explicit parametrisation of the soft and collinear limits. Its constituent building blocks are the universal soft and collinear NLO kernels which factorise from the radiative amplitude in the singular limits. We write where we introduced several notations. Specifically, the prefactor N 1 is defined as where µ is the renormalisation scale and γ E the Euler-Mascheroni constant; {k} is the set of the n + 1 final-state momenta in the radiative amplitude, while {k} / i is the set of n momenta obtained from {k} by removing k i ; when a function takes the argument ({k} / i/ j , k), it depends on the set of n momenta obtained from {k} by removing k i and k j , and inserting their sum k = k i + k j ; finally, B is the Born-level squared matrix element defined in Eq. (2.5), while is the colour-connected Born-level squared matrix element, with T a colour generators, and B µν is the spin-connected Born-level squared matrix element, obtained by stripping the spin polarisation vectors of the particle with momentum k from the Born matrix element and from its complex conjugate.
The NLO soft and collinear kernels are of course well known. In our notation, the eikonal kernel I (i) lm , relevant for soft-gluon emissions, is given by where f i indicates the flavour of parton i, so that δ fig = 1 if parton i is a gluon, and δ fig = 0 otherwise. In order to write the collinear kernels, we begin by introducing a Sudakov parametrisation for the momenta k µ i and k µ j , as they become collinear. We introduce a massless vectork µ , defining the collinear direction, using where k 2 = 2 k i · k j = s ij , and k r is a massless reference vector (for example one of the on-shell momenta of the set {k}, with r = i, j), so thatk 2 = 0. We now write a Sudakov parametrisation of k a (a = i, j), as where we defined the transverse momentak µ a with respect to the collinear directionk, and the longitudinal momentum fractions x a alongk, as The transverse momentak a , for a = i, j, satisfỹ k a ·k =k a · k r = 0 . (2.27) We can now write the spin-averaged Altarelli-Parisi kernels P ij , in a flavour-symmetric notation, as where we defined the flavour delta functions δ f {q,q} = δ f q + δ fq , and δ {fifj }{qq} = δ fiq δ fjq + δ fiq δ fiq . In the following we will use interchangeably the notations P ij , P ij (x i , x j ), or P ij (s ir , s jr ) to denote the collinear kernels of Eq. (2.28), and similarly for the azimuthal kernels Q µν ij and for P µν ij . The Casimir eigenvalues relevant for the SU(N c ) gauge group are C F = (N 2 c − 1)/(2N c ) and C A = N c , consistent with the normalisation T R = 1/2. The azimuthal kernels Q µν ij can be written as We note that the presence of the azimuthal kernels Q µν ij is necessary in order to achieve a local subtraction of phase-space singularities. The collinear kernels satisfy the symmetry properties The final ingredient is the soft-collinear kernel for sector ij, which can be obtained by acting with the soft projector S i on the collinear kernel P ij (indeed, Q µν ij is soft-finite). As detailed in Appendix A, one gets where C fj = C A δ fj g + C F δ fj {qq} . Importantly, the same soft-collinear kernel is obtained also by taking the collinear limit of Eq. (2.23): in other words, the two limits commute, as discussed in detail in Appendix A. Subtracting from the collinear kernels their soft limits, one gets the hard-collinear kernels Although the candidate counterterm K ij defined above contains all phase-space singularities of the product R W ij , with no double counting, the kinematic dependences on the right-hand sides of Eqs. (2.18), (2.19) and (2.20) are not yet suited for a proper subtraction algorithm. Indeed, {k} / i is a set of n momenta that do not satisfy n-body momentum conservation away from the exact S i limit, and, similarly, in the set ({k} / i/ j , k) momentum k = k i + k j is off-shell away from the exact C ij limit. The Born-level squared amplitudes B appearing in the counterterm must instead feature valid (i.e. on-shell and momentum conserving) n-body kinematics for all choices of the n + 1 momenta in the radiative amplitude. A kinematic mapping is thus necessary, in order to factorise the (n + 1)-body phase space into the product of Born (n-body) and radiation phase spaces, thereby allowing one to integrate the counterterms only in the latter. Since the kernels in Eqs. (2.18)-(2.20) are built in terms of Mandelstam invariants, and have not yet been parametrised at this stage, there is still full freedom to choose the most appropriate kinematic mapping in order to maximally simplify the analytic integrations to follow. In particular, at variance with what done in the FKS algorithm, in any given sector one can employ different mappings for different singular limits, or even for different contributions to the same singular limit. In order to take advantage of this freedom, we introduce now a generic Catani-Seymour final-state mapping and parametrisation [33], as follows. Let k a and k b be two final-state on-shell momenta, and let k c be the on-shell momentum of another (massless) parton, with c = a, b. Now one can construct an on-shell, momentum conserving n-tuple of massless momenta {k} (abc) as where s abc = s ab + s ac + s bc , and in particular the condition ensures momentum conservation. Note that the collection of the n light-like momenta {k} (abc) can also be expressed as (2.34) Next, we select different values of a, b, c in different sectors and limits. Consistently with the general structure of factorised virtual amplitudes [70], we treat separately the soft and the hard-collinear limits. For the hard-collinear kernel in sector ij, (C ij −S i C ij ) R W ij , we choose to assign the labels a, b, and c of Eq. (2.32) as a = i, b = j, and c = r: partons i and j specify the collinear sector, while parton r, introduced in Eq. (2.24), is the 'spectator'. For the soft kernel, S i R W ij , we choose to map differently each term in the sum over l, m in Eq. (2.18), with assignments a = i, b = l, and c = m. We then define the local counterterm as where the barred projectors select soft and collinear limits, and assign the appropriate set of on-shell momenta to the kernels. Explicitly where we stress that r = i, j can be chosen differently for different ij pairs, with the constraint that the same r should be chosen for all permutations of ij. The expression in Eq. (2.35) can be rewritten in terms of a sum over sectors of local counterterms K ij , each containing all the singularities of the product R W ij : where it is understood that the action of barred projectors on sector functions is the same as that of un-barred ones, namely S i W ab = S i W ab , and C ij W ab = C ij W ab . To obtain Eq. (2.39) we have used the symmetry under exchange i ↔ j in our definition of C ij R.

Counterterm integration
The counterterm defined in Eq. (2.39) is a sum of terms, each factorised into a matrix element with Born-level kinematics, multiplying a kernel with real-radiation kinematics. The analytic integration of the latter in the radiation phase space proceeds by first summing over all sectors, as done in FKS. This operation matches the fact that the integrated counterterm must eventually cancel the singularities of the virtual contribution, which obviously is not split into sectors. Upon summation over sectors, the integrand becomes independent of sector functions. In fact In the soft term we have considered that the kinematic mapping is j-independent, and performed the sum over j, exploiting the soft sum rule in Eq. (2.12); in the hard-collinear contribution we have used the symmetry of the kinematic mapping and of the collinear operator C ij under the interchange i ↔ j, exploited the collinear sum rule in Eq. (2.12), and the fact that where dΦ (abc) n is the n-body phase space for partons with momenta {k} (abc) , φ is the azimuthal angle between k a and an arbitrary three-momentum (other than k b , k c ), taken as reference direction, and we have set We first consider the integral I hc of the hard-collinear counterterm Each term in the double sum in K hc is parametrised assigning labels a = i, b = j, and c = r, as detailed below Eq. (2.33). We have where ς k indicates the symmetry factor associated to the k-body final state. We note that the integral does not receive any contribution from the azimuthal kernels Q µν ij , as the latter integrate to zero in the radiation phase space. In our chosen parametrisation, the variable z coincides with the collinear fraction jr . The analytic integration of the counterterm is therefore straightforward, and can be carried out exactly to all orders in . By defining one finds where in the last step we replaced the sum over i, j with a sum over 'parent' partons p (which has absorbed the ς n+1 /ς n symmetry factor), carrying momentumk (2.32)), we included a 1/2 Bose-symmetry factor in the C A term, accounting for gluon indistinguishability, and we considered N f light qq pairs. The invariantη pr is defined asη pr =s Notice that the result contains only a single 1/ pole, consistently with the fact that soft singularities are excluded.
Next we turn to the integral I s of the soft counterterm We parametrise it by assigning different labels to each term in the eikonal sum, with a = i, b = l and c = m, as detailed below Eq. (2.33), obtaining In our chosen parametrisation s lm /s im = (1 − z)/z, and s il = ys (ilm) lm : the soft counterterm can then be analytically integrated, once again to all orders in . By defining, for each term of the eikonal sum, we get the simple result We can finally combine soft and hard-collinear integrated counterterms, obtaining, up to O( ) corrections, where we introduced the spin-dependent one-loop collinear anomalous dimension The integrated counterterm in Eq. (2.55) successfully reproduces the pole structure of the virtual NLO contribution (see for example [4]), which provides a check of validity of the subtraction method. Moreover, we note the simplicity of the integrated counterterms to all orders in , which is a direct consequence of having optimally adapted term by term the kinematic mapping and parametrisation.
We conclude this Section with three considerations on the structure of the counterterm. First, the strong coupling α S has been treated as a constant throughout the computation. A dynamical scale for the coupling can simply be reinstated in the counterterm by evaluating it with the Bornlevel kinematics {k}. Second, in the counterterm definition in Eq. (2.35) we have chosen to apply projectors S i and C ij only on the product R W ij , while treating exactly the phase-space measure dΦ rad . In other words, the counterterm phase space is exact, and coincides with that of the realradiation matrix element. We stress that this feature is not crucial to our method: one could as well consider approximate phase-space measures d Φ rad , provided they correctly reproduce the exact dΦ rad in the singular limits. In the massless final-state case detailed in this article, as evident from the above calculation, no computational advantage would result from such an approximation, however the latter may become relevant in more complicated cases. Analogously, restrictions on the counterterm phase space could be applied in order to improve the convergence of a numerical implementation. We leave these possibilities open for future studies.
Third, Eq. (2.39) and Eq. (2.40) are analytically equivalent, but they underpin different philosophies in the implementation of the subtraction scheme. In Eq. (2.39), which is our preferred choice, subtraction is seen as the incoherent sum of terms, each of which features a minimal singularity structure and is separately optimisable, in the same spirit of the FKS method but, we believe, featuring enhanced flexibility. Eq. (2.40), which in what we have presented is employed only for analytic integration, represents a single local subtraction term containing all singularities of the real matrix element, hence it has the same essence as CS subtraction, but with much simpler counterterms. Our method at NLO represents thus a bridge between these two long-known subtraction methods, aiming at retaining the virtues of both, and not being limited by the mutual suboptimal features.

Generalities
The NNLO contribution to the differential cross section with respect to a generic IR-safe observable X can be schematically written as where RR, V V , and RV are the double-real, the UV-renormalised double-virtual, and the UVrenormalised real-virtual corrections, respectively, with In dimensional regularisation, V V features up to a quadruple IR pole in ; RR is finite in d = 4, but it is affected by up to four singularities in the double-radiation phase space, stemming from configurations that feature up to two soft and/or collinear emissions; RV has up to a double IR pole in , originating from its one-loop nature, on top of a double singularity in the single-radiation phase space. The sum of these three contributions is finite due to the IR safety of X and to the KLN theorem. It is however clear that the difficulty of evaluating and integrating complete radiative matrix elements in arbitrary dimension at NNLO is significantly more severe than at the NLO, hence the necessity of a subtraction procedure.
Subtraction at NNLO amounts to modifying Eq. (3.1) by adding and subtracting three sets of counterterms: single-unresolved, double-unresolved, and real-virtual, which we write as and which can be characterised as follows. The single-unresolved counterterm d Φ n+2 K (1) features the subset of phase-space singularities of dΦ n+2 RR which correspond to configurations where only one parton becomes unresolved, analogously to what happens at NLO. The sum d Φ n+2 K (2) + K (12) contains all singularities stemming from kinematic configurations where exactly two partons become unresolved. At NNLO, this exhausts all possible phase-space singularities. We note that the Dirac delta functions associated with these two counterterms mirror their physical meaning, with δ n+1 (X) associated with K (1) , and δ n (X) with (K (2) + K (12) ). The distinction between K (2) and K (12) will be described in detail in Section 3.3. The third subtraction term, d Φ n+1 K (RV) cancels the phase-space singularities of dΦ n+1 RV .
Denoting the corresponding phase-space-integrated counterterms with NNLO cross section can be identically rewritten as In the third line of Eq. (3.5), all terms are separately finite in d = 4, and their sum is finite in the double-radiation phase space, making this contribution fully regular and integrable numerically. In the second line, I (1) features the same poles in as RV , up to a sign, so that their sum is finite in d = 4. The counterterm K (RV) locally subtracts the phase-space singularities of RV ; it contains however explicit poles in , and the local counterterm K (12) is such that the integral I (12) cancels those poles; furthermore, the finite sum RV +I (1) features phase space singularities, and these must be cancelled by the finite sum K (RV) − I (12) . In total, the sum of the four terms in the second line of Eq.

Sector functions
As in the NLO case, we start by partitioning the phase space in sectors, each of which selects the singularities stemming from an identified subset of partons. We thus introduce sector functions W abcd , with as many indices as the maximum number of partons that can simultaneously be involved in an NNLO-singular configuration. We reserve the first two indices for singularities of singleunresolved type, implying that b, c, and d differ from a. As far as double-unresolved configurations are concerned, in particular those of collinear nature, they can involve three or four different partons, hence either indices b, c, and d are all different, or two of them are equal. Without loss of generality we choose the third and the fourth indices to be always different, so that the allowed combinations of indices, that we refer to as topologies, are Since the sector functions must add up to 1, in order to represent a unitary partition of phase space, they can be defined as ratios of the type There is a certain freedom in the definition of σ abcd . Analogously to the NLO case, we design them in such a way as to minimise the number of IR limits that contribute to a given sector. In addition, at NNLO there is another property to be required, new with respect to NLO, and related to the fact that the integrated single-unresolved counterterm I (1) must be combined with the real-virtual contribution, to cancel its explicit poles in , as detailed in Section 3.1. Since RV , as any term with (n + 1)-body kinematics, is split into NLO-type sectors, the same must be true for I (1) . This implies that, roughly speaking, sector functions with four indices must factorise sector functions with two indices in the single-unresolved limits, in order for the cancellation of poles to take place NLO-sector by NLO-sector.
A possible expression for σ abcd with the required properties is With the sector functions defined in Eq. (3.7) and Eq. (3.8), the list of singular limits acting non-trivially in each NNLO sector includes the single-unresolved projectors S a and C ab , already considered at NLO, as well as the following double-unresolved limits: C abc : w ab , w ac , w bc → 0 , w ab /w ac , w ab /w bc , w ac /w bc → constant (uniform double-collinear configuration of partons (a, b, c)) , (ordered soft (first) and collinear configuration of partons a and (b, c)) , (ordered collinear (first) and soft configuration of partons (a, b) and c) .
Notice that only the first two limits of the list (3.9) are genuinely double-unresolved 2 , namely they cannot be reduced to compositions of single-unresolved limits when acting on the double-real matrix elements; the remaining three configurations are compositions of single-unresolved limits when acting on matrix elements, but not when they are applied to the sector functions in Eq. (3.7), therefore they have to be introduced as independent limits. In Appendix B we show that, among the single-and double-unresolved limits that we are considering, only a subset give a non-zero contribution in the various topologies. They are In Appendix B we also show that all the limits reported in Eq. (3.10) commute when acting on the sector functions, and that the combinations of these limits exhaust all possible single-and doubleunresolved configurations in each sector. We stress that this structure depends on our choice of sector functions; with other functions the surviving limits would in general be different.
It is now necessary to study the properties of the sector functions defined in Eq. (3.7) and Eq. (3.8) under the action of single-unresolved limits. As noted above, in these configurations the NNLO sector functions must factorise into products of NLO-type sector functions. To this end, let us define so that the NLO sector functions in Eq. (2.14) are given by W ij = W (11) ij , and similarly σ ab = σ (11) ab . One easily verifies that the functions W (αβ) ij satisfy all the requirements that must apply to NLO 2 In the literature the configuration C abc is sometimes referred to as triple-collinear. We call it double-collinear, following [71], in order to consistently specify the type of configuration as being double-unresolved, rather than indicating the number of partons that become collinear. sector functions. It is now straightforward to verify that the NNLO sector functions defined in Eq. (3.7) and Eq. (3.8) satisfy where W [ab]c is the NLO sector function defined in the (n + 1)-particle phase space with respect to the parent parton [ab] of the collinear pair (a, b). Finally, the NNLO sector functions satisfy sum rules analogous to the NLO ones in Eq. (2.12), and which stem from their definition in Eq. (3.7). One may verify that where by π(ijk) we denote the set {ijk, ikj, jik, jki, kij, kji}. Sum rules for composite doubleunresolved limits, that follow from those reported in Eqs. (3.13)-(3.15), will be further detailed in Section 3.5, where we describe the structure of the double-unresolved counterterm. We stress that the properties in Eqs. (3.13)-(3.15), in full analogy with the NLO case, allow one to perform sums over all the sectors that share a given set of double-unresolved singular limits, eliminating the corresponding sector functions prior to countertem integration. This feature, distinctive of our method at NNLO, is crucial for the feasibility of the analytic integration of counterterms.

Definition of local counterterms
As reported in Eq. (3.10), a limited number of products of IR projectors is sufficient to collect all singular configurations of the double-real matrix elements in each sector. By subtracting these products from the matrix element, one gets, for the different topologies, the finite expressions where we separated the action of the single-unresolved limits L The order with which the various operators are applied to matrix elements is irrelevant, as all limits commute. In Appendix B we show that this property is also respected by the sector functions defined in Eq. (3.7). Candidate double-real local counterterms for the various topologies T can thus be defined, in analogy with Eq. (2.16), as The different contributions are naturally split according to their kinematics. All terms containing only single-unresolved limits are assigned to K (1) , the single-unresolved counterterm; terms containing only double-unresolved limits are assigned to K (2) , which we refer to as pure doubleunresolved counterterm; all remaining terms, containing overlaps of single-and double-unresolved limits, while still featuring double-unresolved kinematics, are assigned to K (12) , which we refer to as mixed double-unresolved counterterm. A direct characterisation of mixed double-unresolved counterterms in terms of factorisation kernels will be discussed in Ref. [70]. We write therefore, for each topology T , . This correspondence is strict: indeed, if one imagines removing from a given process all n-body contributions, for instance by means of phasespace cuts, the original NNLO computation reduces to the NLO computation for the process with n + 1 particles at Born level, with RR playing the role of single-real correction, and RV that of virtual contribution; in this scenario, K (1) becomes exactly the candidate NLO local counterterm. As for the double-unresolved contributions, K (2) is to be integrated in d Φ rad,2 , giving rise to up to four poles in , multiplied by Born-like matrix elements, analogously to V V ; the single-unresolved structure in K (12) , on the other hand, makes it suitable for integration in d Φ rad,1 ; once this is achieved, its double-unresolved projectors naturally become single-unresolved projectors for the parent parton which originated the first splitting, thus reproducing the structure of K (RV) . This is necessary, since the integral of K (12) must compensate the explicit poles in of K (RV) . This cancellation also relies on the factorisation properties of sector functions, presented in Eq. (3.12), as will be further detailed below. The double-unresolved kernels appearing in the counterterm definitions of Eqs. (3.19)-(3.21) can be derived from soft and collinear limits of scattering amplitudes, which are universal, and for the massless case relevant to this article they were computed in Refs. [20,22]. General expressions for the kernels can also be derived starting from the factorisation of soft and collinear poles in virtual corrections to fixed-angle scattering amplitudes, as will be discussed in detail in Ref. [70].
Here we just write symbolically In the double-soft limit, B cdef is the doubly-colour-connected Born matrix element, defined for instance in Eq. (113) of [22]; the eikonal kernels I  96) and (110) of [22] 3 . In the non-factorisable double-collinear limit C ijk , the set of momenta ({k} / i/ j / k , k) refers to a set of n partons obtained from {k} by removing k i , k j , and k k , and inserting their sum k = k i + k j + k k . The expressions for the double-collinear spin-averaged kernels P ijk and for the azimuthal kernels Q µν ijk , all symmetric under permutations 4 of i, j, and k, can be easily extracted from [20,22], but the expressions are long and therefore will not be reproduced here. We note however that Q µν ijk can always be cast in the form where, in analogy with Eq. (2.26), 27) and k µ r is a light-like vector which specifies how the collinear limit is approached. The Lorentz structure in Eq. (3.26), identical to the NLO one in Eq. (2.29), is such that the radiation-phasespace integral of the double-collinear azimuthal terms vanishes identically. Hence, once more, the analytic integration of the counterterms involves only spin-averaged kernels. The factorisable double-collinear limit C ijkl features the doubly-spin-correlated Born matrix element B µνρσ , with a kinematics obtained from {k} removing k i , k j , k k , and k l , and inserting the sums k ij = k i + k j , and k kl = k k + k l ; the corresponding kernel is defined as Finally, the soft-collinear limit SC ijk features a colour-and spin-correlated Born contribution B cd µν , obtained from the colour-correlated Born matrix element B cd by stripping external spin polarisation vectors.
We now note that, while Eqs. (3.19)-(3.21) are quite natural, they contain a certain degree of redundancy. In fact, the double-real matrix element RR can feature at most four phase-space singularities, hence not all of the projectors relevant to a given topology, listed in Eq. (3.10), carry independent information on its singularity structure. These redundancies can be eliminated by exploiting the idempotency of projection operators: for instance, once SC icd has been applied to the double-real matrix element, further action on the latter by S i does not produce any effect, and analogously if the limit C ij is applied after the action of CS ijk . This ultimately stems from the factorisable nature of SC icd RR, and of CS ijk RR, namely Even if this factorisation property does not hold when the SC icd and CS ijk limits are applied to the sector functions of Eq. (3.7), the commutation relations discussed in Appendix B are sufficient to prove that 5 As a consequence of Eq. (3.30), the candidate mixed double-unresolved counterterm K (12) simplifies to 31) and the sum of K free of any contribution from the limit L valid both on matrix elements and on sector functions, to rewrite After the simplifications just discussed, we are finally in a position to write down the definition of 5 Also the limit C ijkl has a factorisable nature when applied on the double-real matrix element, however, in this case, the relevant commutation relations are not sufficient to obtain the analogue of Eq. (3.30).
the candidate local counterterms for all contributing topologies T = ijjk, ijkj, ijkl: The final step for the construction of the NNLO counterterms, analogously to what happens in the NLO case discussed in Section 2.3, is to apply kinematic mappings to Eq. (3.35). There is ample freedom in the choice of these mappings, and in principle different mappings can be employed for different kernels, or even for different contributions to the same kernel. The detailed definition of the kinematic mappings we employ for each counterterm is given in Sections 3.4 and 3.5 where, as usual, all remapped quantities will be denoted with a bar. Finally, the real-virtual counterterm has formally the same structure as the NLO counterterm of Eq. (2.39), with the replacement R → RV , and will be sketched in Section 3.6.

Single-unresolved counterterm
We start by separating the hard-collinear and the soft contributions to the candidate singleunresolved counterterm: Using the factorisation properties (3.12) we can proceed as done at NLO. We define the appropriate counterterms with remapped kinematics, where in this case barred projectors apply not only to matrix elements, but also to sector functions: The kinematic mapping of sector functions, once the integrated counterterm is considered, allows to factorise the structure of NLO sectors out of the radiation phase space, and integrate analytically only single-unresolved kernels. Explicitly where R ab and R µν are the colour-and spin-correlated real matrix elements and if both k = j and l = j, the same r should be chosen for all permutations in π(π(ij) π(kl)).

Integration of the single-unresolved counterterm
As done at NLO, we now integrate the single-unresolved counterterm in its radiation phase space. We first get rid of the NLO sector functions W (αβ) ij using their NLO sum rule, obtaining fully analogous to its NLO counterpart in Eq. (2.50). The integral of K (1, s) similarly yields where, in the last step, all identical soft-gluon contributions have been remapped on the same real kinematics {k}, and the sum i δ fig has absorbed the symmetry factor ς n+2 /ς n+1 . The combination of hard-collinear and soft contributions is straightforward, as in the NLO case, yielding where indices h and q run over the NLO multiplicity, barred momenta and invariants refer to NLO kinematics, and r = a. Eq. of the sum RV + I (1) . It is important to note, however, that in Eq. (3.49), as well as in RV , the full structure of NLO sector functions W hq is factorised in front of the integrated singularities, which means that the cancellation of 1/ poles between RV and I (1) occurs sector by sector in the (n + 1)-body phase space.

Double-unresolved counterterm
The double-unresolved counterterm with n-body kinematics consists of two parts: the pure doubleunresolved counterterm K (2) , which must be integrated in the double-radiation phase space, and the mixed double-unresolved counterterm K (12) which must be integrated in a single-radiation phase space. From Section 3.1 we see that, while their integration has to be performed independently, the non-integrated counterterms K (2) and K (12) appear only combined in the last line of Eq. (3.5).
Owing to the simplifications discussed at the end of Section 3.3, the sum K (2) + K (12) is much simpler than the two terms taken separately, and it reads Exploiting Eqs. (3.12), together with one is able to recast the above expression in a form that explicitly features only sums of sector functions that add up to 1, according to the sum rules of Eqs. (2.12), (3.13)-(3.15), and to as well as Introducing remapped kinematics for the double-real matrix element and for the sector functions W ab , analogously to what done for the single-unresolved counterterm, the double-unresolved counterterm finally reads We stress that in each contribution the kinematics of the double-real matrix element undergoes a different mapping onto the Born one, so as to maximally adapt the parametrisation of the integrands to the kinematic invariants that naturally appear in the respective kernels. The explicit definition of the barred limits appearing in the first three lines of Eq. (3.54) is where we have introduced the remapping where the same r = i, j and r = k, l should be chosen for all permutations in π(π(ij) π(kl)). We have introduced the remapping The remaining composite limits of RR appearing in Eq. (3.54) are listed in Appendix C.

Integration of the mixed double-unresolved counterterm
The mixed double-unresolved counterterm features n-body kinematics but, peculiarly, it needs to be integrated analytically only in the phase space of a single radiation. This operation is necessary to show that such an integral features the same explicit 1/ singularities as the K (RV) counterterm, and, at the same time, it features the same phase-space singularities is I (1) .
We start by considering the hard-collinear contribution to K (12) . Following Eqs. (3.35) we have We stress that in the last expression we have kept the CS ijk terms: these cancel out in the sum K (2) + K (12) , but do contribute to the integrals I (2) and I (12) , which have to be evaluated separately.
The explicit computations reported in Appendix D show that the phase-space integral I (12, hc) of the hard-collinear contribution can be recast in the simple form where the integral J hc ij is defined in Eq. (2.49), and the barred limits on R are given by jr , ), do not hamper numerical integrability. We now consider the K (12, s) counterterm, which is obtained combining the soft contributions of the last three equations of (3.35). The result is The explicit computations reported in Appendix D show that the phase-space integral I (12, s) of the soft contribution can be recast as I (12, s) where the integral J s is defined in Eq. (2.53), and the limits in this case are defined by where r = k, l, and the same r should be chosen for kl and for lk. The same considerations that were applied below Eq. (3.65) hold in this case as well, referring now to the comparison between Eq. (3.69) and Eq. (3.48). Combining soft and hard-collinear contributions, the final expression for the integrated counterterm I (12) is where the soft and collinear limits are meant to be applied on matrix elements and on sector functions, but not on the logarithms lnη ij , while barred momenta and invariants refer to NLO kinematics, and finally one must choose r = a.
Since I (12) collects the same phase-space singularities as I (1) , and I (1) in turn features the same explicit 1/ poles as RV , it follows by construction that I (12) also contains the same 1/ poles as K (RV) , as necessary in order to compute the second line of Eq. (3.5) in d = 4. We stress that these considerations hold separately in each NLO sector W hq .

Integration of the pure double-unresolved counterterm
The candidate pure double-unresolved counterterm, summed over NNLO sectors, follows from Eq. (3.35) and reads We work on this expression by symmetrising indices, and exploiting the sum rules in Eqs. (3.13)-(3.15), as well as Eq. (3.53), together with Introducing remapped kinematics for the double-real matrix element, the pure double-unresolved counterterm can be finally cast in the form which is manifestly free of NNLO sector functions. The counterterm in Eq. (3.76) is thus suitable for analytic integration over the double-unresolved phase space, upon definition of the barred limits. First, we note that the barred limits appearing in the first line of Eq. (3.76) have already been defined in Eqs. (3.55)-(3.57). Next, we consider all terms in Eq. (3.76) containing the four-particle double-collinear barred limits C abcd . Their contribution can be rewritten as Defining the barred limits in terms of soft and collinear kernels, Eq. (3.77) becomes (3.80) Note that K (2) only involves simple combinations of soft and collinear kernels, all remapped in an optimal manner so as to make their analytic integration as straightforward as possible. The complete integration of the pure double-unresolved counterterm, along with other details of the implementation, will be presented in a forthcoming publication. Here we will limit ourselves, for the sake of illustration, to the computation, in Section 4, of the subset of the terms that enter the T R C F contribution to e + e − → qq at NNLO.

Real-virtual counterterm
The real-virtual NNLO contribution RV features a structure of explicit poles dictated by its nature of virtual one-loop matrix element, namely where the indices k and l run over real-radiation multiplicities, and G( ) denotes the collection of terms that are non-singular in the → 0 limit, encoding process-specific information.
The corresponding real-virtual counterterm K (RV) contains all phase-space singularities appearing in Eq. (3.81). Analogously to what done at NLO in Eq. (2.39), it is defined as In this paper we do not aim at giving a final expression for the integrated real-virtual counterterm I (RV) , which will instead be detailed in a subsequent publication, together with the completion of the integrals contributing to I (2) ; we limit ourselves to stressing that such an analytic integration is of a comparable or lower complexity with respect to that of the pure double-unresolved counterterms, hence it does not pose any new significant computational challenges. Indeed, as far as the -singular contributions in Eq. (3.81) are concerned, they are proportional to real or colour-connected real matrix elements, hence their IR limits in Eq. (3.82) involve single-soft and single-collinear kernels of NLO-level complexity. The structure of the -finite remainder G( ) is slightly subtler: it can be further split into the sum of a process-specific regular contribution, plus a universal phase-spacesingular term. The IR limits of the latter, in particular, involve kernels which represent integrands of a higher complexity than the NLO ones, but still can be handled analytically in full generality. We leave the completion of these contributions to future work.

Proof-of-concept calculation
In order to demonstrate the validity of our local subtraction method, in this Section we apply it to di-jet production in electron-positron annihilation, as a test case. We consider radiative corrections up to NNLO, restricting our analysis to the contributions proportional to T R C F . The production channels available in this case are B, V, V V : e + e − → qq , R, RV : e + e − → qq g , RR : e + e − → qq q q .
where y and z are the Catani-Seymour variables relative to the secondary-radiation phase space, and x parametrises the azimuth between subsequent emissions. In the chosen parametrisation, four out of the six involved binary invariants have simple expressions, while the remaining two involve square roots related to azimuthal dependence. The explicit expressions are s ab = y y s abcd , where, for the process at hand, the invariant s abcd =s (abcd) cd coincides with the squared centreof-mass energy s. In this parametrisation, all integrations for the process we are considering are straightforward. For the case of double-soft radiation the relevant integral is [22] dΦ rad,2 S ij RR = N 2 , (4.10) where {ij} = {34}, according to Eq. (4.6). Different terms in the eikonal sum can be remapped to the same Born kinematics, and, performing the relevant colour algebra, the result is (4.11) The double-collinear contribution (before the subtraction of the soft-collinear region) can be similarly computed, and it yields where, following [20,22], we have set . The subtraction of the double-counted soft-collinear limit is very simple in this case, since one has dΦ rad,2 S ij C ijk RR = dΦ rad,2 S ij RR , (4.14) as can be deduced from Eq. (3.55) and Eq. (3.57) in the case of two soft quarks, in a process featuring only two partons at Born level, identified here with k and r. Adding up all contributions to the pure double-unresolved integrated counterterm, we get Next, we consider the integration of the single-unresolved counterterm, applying the general formula, Eq. (3.49), and restricting our analysis to the case in which only the single-collinear limit is non-zero. We find where the real-radiation matrix element R involves n + 1 = 3 particles, the indices h and q take values in the set {1, 2, 3 ≡ [34]}, and we can choose r = 1 or r = 2 when h = 1, q = 2, while r = 2 − h in the other cases. The result in Eq. (4.16) must be combined with the RV contribution, and we can explicitly check that their sum is finite in d = 4, sector by sector in the NLO phase space. Indeed The next ingredient is the mixed double-unresolved contribution, which can be read off the general formula, Eq. (3.73). In sector hq it reads The combination of Eq. (4.18) with the real-virtual local counterterm in the same NLO sector must be finite in d = 4. Indeed we find that The final ingredient for subtraction is the integral of the real-virtual counterterm. In the present case, it is given by where I C F , n=2 denotes the NLO counterterm given in Eq. (2.55), considered in the particular case of two non-gluon final-state partons at Born level. All required ingredients for NNLO subtraction for the process at hand are now assembled, and we can proceed to a numerical consistency check.

Collection of results
The heart of the subtraction procedure is the combination of analytic results with numerical integration of the finite remainder of the real-radiation squared matrix element, to get physical distributions and cross sections. For this proof of concept, we will simply reconstruct numerically the total cross section for the production of two quark pairs of different flavours. We emphasise however that the formalism we constructed is completely general and local: a detailed numerical implementation for all processes involving only final state massless partons is being developed and will be presented in forthcoming work. The cross section is constructed in general, as shown in Eq. (3.5), as a sum of three finite and integrable contributions, given by The subtracted double-virtual contribution is computed analytically, and is finite in d = 4. In this case, it is given by where, for definiteness, in the second line we have randomly chosen µ 2 /s = 0.35. For real radiation, we have written a Monte Carlo code to integrate numerically the remaining two terms in Eq. (4.21), obtaining The rescaled NNLO correction, evaluated numerically by means of the subtraction method, is then (4.24) to be compared with the analytical result For completeness, we also show in Fig. 1 that also the logarithmic renormalisation-scale dependence is correctly reproduced with the same accuracy.

Conclusions
In this work we have presented a new scheme to perform local analytic subtraction of infrared divergences up to NNLO in QCD. The method has for now been developed and applied to processes featuring only massless partons, and not involving coloured partons in the initial state, as a first significant step towards a general formulation. Our subtraction procedure is conceived with the aim of minimising complexity in the definition of the local IR counterterms, aiming for their complete analytic integration in the unresolved phase space, and working towards an optimal organisation of the numerical integration of the observable cross section.
Our local IR counterterms are defined through a unitary partition of the phase space into sectors, in such a way as to isolate in each sector a minimal number of phase-space singularities, associated with soft and collinear configurations of an identified set of partons (up to two at NLO, and up to four at NNLO). In each sector, the counterterms are built out of a collection of universal kernels, written in terms of kinematic invariants, which can be defined in terms of gauge-invariant operator matrix elements, as detailed in [70], or can be obtained as limits of radiative matrix elements in the dominant soft and collinear configurations. Overlapping singularities are fully taken into account by suitable compositions of such singular limits, with no need to resort to sector-decomposition techniques.
The sector functions that realise the phase-space partition are engineered in such a way as to satisfy fundamental relations that allow to achieve the main goals of the method. A number of sum rules, stemming from the definition of the sector functions, allow one to recombine various subsets of sectors, prior to performing counterterm integration, eventually yielding integrands that in all cases are solely made up by sums of elementary infrared and collinear kernels. Moreover, through factorisation relations, NNLO sector functions reproduce the complete structure of NLO sectors in all relevant single-unresolved limits, allowing to subtract, sector by sector in the NLO phase space, the singularities of the NNLO contributions featuring NLO kinematics.
The kinematic mappings necessary for phase-space factorisation, as well as the parametrisations of the radiation phase space over which the counterterms are integrated, are devised by maximally exploiting the freedom one has in their definition. They are not only chosen differently for different sectors, but also, importantly, for different counterterm contributions in the same sector. This allows us to employ parametrisations that are naturally adapted to the kinematic invariants that appear in each singular contribution, yielding simple integrands to be evaluated analytically.
In this article we have integrated all needed counterterms over the exact phase-space measures, without exploring the possibility of approximating the latter in the relevant soft and collinear limits. While this possibility would not have resulted in any analytic simplification in the cases considered here, this might instead be the case for general hadronic reactions (for example when including initial-state partons, or for a generalisation to the massive case). This possibility will be investigated in dedicated future studies, which are beyond the scope of the present paper.
At NLO, we have shown that the proposed subtraction method works in the general case of massless QCD final states, with the integrated counterterms reproducing analytically the full structure of virtual one-loop singularities. Moreover, as a test of the power of the method, we have shown that the NLO counterterm integration can be performed exactly to all orders in the dimensional regulator , which bears witness to the extreme simplicity of the integrands involved.
At NNLO, we have deduced the structure of the subtraction scheme in full generality for massless QCD final states. All single-unresolved and mixed double-unresolved counterterms of doublereal origin have been integrated analytically to all orders in , as simply as in the NLO case, and the properties of sector functions have allowed us to show that these integrals correctly reproduce, sector by sector, the explicit poles and phase-space singularities of real-virtual contributions. We stress that this is a highly non-trivial test of the consistency of the scheme, and of the delicate organisation of different contributions to the cross section. As for double-unresolved counterterms, we have deduced their structure in general, and performed the relevant integrations in a proof-ofconcept case, the T R C F contribution to e + e − → qq at NNLO, which has been detailed explicitly. Integrals relevant for the full set of counterterms for generic final-state massless configurations are not expected to present special difficulties, nor to pose an obstacle to a more general application of the method.
This article represents a first step towards the formulation of a general, local, analytic, and minimal subtraction scheme, relevant for generic multi-particle hadronic processes at NNLO in QCD. To reach this goal, a number of important steps still need to be taken, including the analytic integration of the remaining double-unresolved counterterms for final-state processes, the generalisation to include initial-state massless partons, and the extension to the massive case, as well as the completion of an efficient computer code implementing the subtraction method in a fully differential framework. We believe however that the present work lays a solid foundation for these future developments.

A Commutation of soft and collinear limits at NLO
In this Appendix, as an example, we explicitly show the commutation of the soft and collinear limits S i and C ij , and, in the process, deduce the form of the soft-collinear kernel S i P ij appearing in Eq. (2.20). The action of operators S i and C ij on ratios of elementary massless invariants s ij is given by We start by verifying that the sequential action of the singular projectors on sector functions does not depend on their ordering. To this end note that where in Eq. (A.3) we used the fact that only l = j gives rise to a singular contribution 1/w il in the collinear limit, while in Eq. (A.4) we have noted that e i → 0 in the soft limit. Next, we consider the action of the composite projector S i C ij on the physical real-radiation amplitude squared, where, without loss of generality, we drop all kinematic dependences in the real and Born-like matrix elements. Starting from Eq. (2.19) we find We now note that Q µν ij , defined in Eq. (2.29), is not singular in the soft limit for parton i, hence S i Q µν ij = 0. The same happens for all terms in P ij which do not contain a denominator 1/x i . We now rewrite the remaining contributions in terms of Mandelstam invariants, using the definition of x i and x j in Eq. (2.26), with the result where the ellipses denote terms that remain regular as parton i becomes soft. Taking now the S i limit according to Eq. (A.1), we get which is Eq. (2.30). In particular, we note that the soft limit S i does not correspond to taking x i → 0, rather to taking s ir → 0 (the two definitions differ by subleading soft terms). The softcollinear limit is thus We can now verify commutation by considering the two singular limits in reversed order. We find Among all the terms in the double sum, only those with k = j or l = j are singular in the collinear limit, hence According to property (A.2), in the collinear limit C ij the ratio s jl /s il is independent of l: we can therefore set l = r and get where in the last two steps we have used colour algebra, and the definition of the Casimir operator

B Soft and collinear limits of sector functions
In this Appendix we explore the properties of the NNLO sector functions defined in Eqs. (3.7) and (3.8), including their relation to the NLO-like functions defined in Eq. (3.11). We begin by establishing which limits, among S a , C ab , S ab , C abc , C abcd , SC abc , and CS abc , are non-vanishing in the three sector topologies W ijjk , W ijkj and W ijkl . To this end, we start by analysing the behaviour of the sector-function denominator σ (see Eq. (3.7)), in these limits. We find where [ij] denotes the parent parton of i and j and we have used the definition of the NLO-like sector functions in Eq. (3.11). Now we note that a singular limit L gives a non-zero result, when applied to the sector functions W abcd , only if the numerator of the latter, σ abcd , appears as one of the addends of L σ. Inspection of Eq. (B.1) then proves that the limits reported in Eq. (3.10) exhaust the surviving ones in each sector.
Next, we show that all of the limits in Eq. (3.10) commute when acting on σ. This is a crucial step for our method, since commutation of limits drastically reduces the number of independent configurations one needs to explore. Furthermore, one must note that, while commutation can be understood from physical considerations when limits are taken on squared matrix elements, sector functions are a crucial but artificial ingredient of our method, and commutation of limits is nontrivial in this case. We list below all relevant ordered limits, acting on the denominator function σ, beginning with those involving the single-soft limit S i .
Next, we list ordered limits involving the single-collinear limit C ij , and not considered above.
Moving on to ordered limits involving the double-soft limit S ab , and not considered above, we find Coming to double-collinear limits of type C ijk and C ijkl , we get Finally, the mixed soft-collinear limits SC ijk and CS ijk satisfy The relations in Eqs. (B.2)-(B.6), where the limits are applied to the sector-function denominator σ, are sufficient to prove that all non-vanishing limits in the different topologies commute when acting on the sector functions. The same commutation relations hold when applied to the physical double-real matrix elements, as can be proved analogously to what was done in Appendix A. The next step in our analysis is to prove that the compositions of the limits given in Eq. (3.10) exhaust all single-and double-unresolved configurations in each sector. In other words, there are no leftover singular phase-space regions after all combinations of limits in Eq. (3.10) have been applied. We start by denoting with L i a generic set of soft and collinear limits, corresponding to configurations where some physical quantities λ i , which could be collections of energies, or angles, or similar, approach zero. Compositions of two (or more) such limits can be either 'uniform' or 'ordered', with the two cases being defined as ⇐⇒ uniform composition of L i and L j ; All single-and double-unresolved configurations in each sector can then be systematically generated by combining in all possible ways the single-soft and single-collinear limits selected by the sector functions, namely S a , S c , C ab , and C cd 6 in sector W abcd . Owing to the prescription α > β > 1 in Eq. (3.8), the action on σ of a uniform composition involving soft and collinear limits is equivalent to the corresponding ordered composition where the soft limits act first: where L s (L c ) are collections of soft (collinear) limits, while L, and L are generic combinations of limits. The remaining uniform compositions involve either a pair of collinear or a pair of soft limits 7 , which can be directly identified with the limits given in Eq. (3.10): We conclude that all possible single-and double-unresolved singular configurations can be obtained as ordered compositions without repetition 7 of the limits • S i , S j , C ij , C jk , S ij , and C ijk for topology W ijjk ; • S i , S k , C ij , C jk , S ik , and C ijk for topology W ijkj ; • S i , S k , C ij , C kl , S ik , and C ijkl for topology W ijkl .
To conclude, we reduce this list of limits, topology by topology, to that given in Eq. (3.10).
• Topology W ijjk According to Eqs. (B.2)-(B.6), the S j limit commutes with all other limits in the list except S i . Therefore, when appearing in a generic composition of limits, it can be moved to the right until it encounters S i . At this point one can use valid for generic limits L and L , to remove S j . If S i is not present at the right of S j , the latter can be moved to the rightmost position, where it vanishes: Since the action of S j either gives zero or can be replaced by that of S ij , S j can be simply removed from the list. 6 Note that compositions of limits involving both C ij and C jk automatically also involve the limit C ik . Indeed 7 Repeated limits can in all cases be readily simplified. Given a generic limit L, one has for example Considering now C jk , we note that it commutes with C ijk and with S ij , and it satisfies so that C jk can either be moved to the rightmost position, where it gives zero, or replaced by C ijk or SC ijk . Consequently, one can remove C jk from the list of limits, and add SC ijk in its stead. The list of singular limits is thus reduced to the first line of Eq. (3.10), Besides commuting with C jk , S ik , and C ijk , the single-soft limit S k satisfies Since S k can be either moved to the rightmost position, where it gives zero, or replaced by S ik or CS ijk , one can remove it from the list of contributing limits. A similar statement holds for C jk , which commutes with S ik , and C ijk , and satisfies As a consequence, C jk can either be moved to the rightmost position, where it gives zero, or replaced by C ijk or SC ijk . The list of singular limits in sector W ijkj can thus be reduced to the second line of Eq. (3.10), The discussion of the S k and C kl limits holds unchanged with respect to the one relevant for S k and C kj in topology W ijkj . These limits can either be moved to the rightmost position, where they yield zero, or be replaced by limits that are already present in the list, (S ik or CS ijk in the case of S k , C ijkl or SC ikl in the case of C kl ). The final list of contributing limits thus coincides with the third line of Eq. (3.10),

C Composite IR limits of the double-real matrix element
In this Appendix we list the composite soft and collinear limits of the double real-radiation squared matrix element needed for the evaluation of the double-unresolved counterterm K (2) + K (12) in Eq. (3.54), including the detailed dependence on the remapped phase-space variables described in Section 3.5. The remappings described in the following apply also to the corresponding sector functions W ab in Eq. (3.54). First, we consider composing a double-collinear limit and a collinear limit. We find where the same r = i, j, and r = j, k should be chosen for all permutations of ijk.
Composition of a double-soft limit and a collinear limit (on the same pair of particles) yields , and the same r = i, j should be chosen for all permutations of ijk (recall that index k appears in the sector function associated with these contributions in Eq. (3.54)). Next, we consider the composition of a four-particle double-collinear limit with a soft limit on one of the collinear particles. We get S i C ijkl RR = 2 N 1 C fj I where the same r = i, j and r = k, l should be chosen for all permutations in π(π(ij) π(kl)).
Taking successively a double-soft limit and a single-soft limit, we get Finally, composing a double-collinear limit and a single-soft limit, we get where r = j, k, and the same r should be chosen for jk and for kj.

D Results for the mixed double-unresolved counterterm
In this Appendix we show explicitly how the terms in the integrated mixed double-unresolved counterterm organise themselves in the form of single-unresolved limits in the NLO phase space. Starting from Eq. (3.63), using Eqs. (3.12), and (3.51), and introducing remapped kinematics for the double-real matrix element and for the sector functions W ab , the hard-collinear contribution to the mixed double-unresolved counterterm can be cast in the form Using the NLO sector-function sum rules, and appropriate symmetrisations, the latter becomes K (12, hc) = − i, j>i k =i,j C jk W jk + W kj C ij C ijk + l =i,j,k The singular limits in Eq. (D.2), as well as their phase-space integrals, are explicitly computed in the following. For brevity, in all contributions to the hard-collinear counterterm we do not display kinematic dependences, writing P hc ij for P hc ij (s ir , s jr ), and similarly for Q µν ij , while the real matrix element is written as R ≡ R {k} (ijr) , and similarly for R µν . We note that all limits are accompanied by single-and double-soft subtractions, guaranteeing the hard-collinear character of the counterterm. Terms containing C ij S ij give where the hard-collinear integral J hc ij is defined in Eq. (2.49), and we exploited the fact that azimuthal terms integrate to zero in the radiation phase space. The soft-collinear limit CS ijk contributes The nested collinear limit C ijk C ij contributes jr , C jk R . (D.5) The nested collinear limit C ij C ijkl contributes Combining Eq. (D.5) with a double-soft limit we get Acting on Eq. (D.4) with a three-particle double-collinear limit one finds Finally, replacing the three-particle double-collinear limit in Eq. (D.8) with the four-particle one we get the result which can be straightforwardly rewritten as Eq. (3.64). We next turn to the soft term in Eq. (3.68). Using Eq. (3.12), together with SC ikl C ijkl RR = S i C ijkl RR , (D.11) and introducing, as usual, remapped kinematics for the sector functions and for the limits of the matrix element, we obtain the expression By means of Eq. (3.52), and renaming indices, we finally get The singular limits in Eq. (D.13), as well as their phase-space integrals, are explicitly computed below. For brevity, in the following we set R ab ≡ R ab {k} (iab) unless explicitly stated otherwise. Let us begin by considering the iteration of a soft limit and a double-soft limit. We find cd , S k R cd , (D.14) where the soft integral J s is defined in Eq. (2.53). Next, we consider the non-collinear projection of the soft-collinear limit SC ijk , obtaining