Geometric IR subtraction for final state real radiation

A scheme is proposed for the subtraction of soft and collinear divergences present in massless final state real emission phase space integrals. The scheme is based on a local slicing procedure which utilises the soft and collinear factorisation properties of amplitudes to produce universal counter-terms whose analytic integration is relatively simple. As a first application the scheme is applied to establish a general pole formula for final state real radiation at NLO and NNLO in Yang Mills theory for arbitrary multiplicities. All required counter-terms are evaluated to all orders in the dimensional regulator in terms of Γ — and pFq hypergeometric — functions. As a proof of principle the poles in the dimensional regulator of the H → gggg double real emission contribution to the H → gg decay rate are reproduced.


Introduction
As the LHC is entering its high precision phase the need for precise calculations of higher order perturbative corrections in QCD is becoming ever more important. A major bottleneck in such calculations is due to infrared (IR) divergences, which are related to non-integrable singularities in scattering amplitudes which arise in soft and/or collinear momentum configurations. While the KLN theorem [1,2] guarantees that these singularities cancel (in JHEP08(2018)006 the sum over real and virtual emissions for final state radiation) collinear divergences associated to partons in the initial states are handled instead by a renormalisation of the parton densities [3][4][5].
Since the IR divergences can be regulated dimensionally their cancellation is less problematic for analytic calculations of total inclusive quantities where one integrates over the entire phase space volume. The situation is different for differential quantities, where the domain of integration is taken over an arbitrary IR safe region of phase space. Then it may be unfeasible to carry out the integral analytically and a subtraction procedure is required to render the integrand finite. To accomplish this task has inspired the construction of a large number of different subtraction procedures.
Existing methods fall into either one of two conceptually quite different approaches. These are referred to as subtraction and slicing methods. The subtraction method is based on making the divergent integrand finite by subtracting from it a suitable counterterm whose singular behaviour matches that of the original integrand. This counter-term is subsequently added back in, analytically or numerically, integrated form. Methods designed for dealing with real emissions at the next-to-leading order (NLO), which fall into this category, are the Catani-Seymour (CS) dipole method [6,7], the Frixione-Kunszt-Signer (FKS) subtraction method [8,9] as well as the Nagy-Soper subtraction method [10]. While all of these methods rely on the universal soft and collinear factorised limits of amplitudes to construct suitable counter-terms, they are implemented in quite different ways.
In the FKS method sets of non-overlapping divergences are separated by a partition of unity in the form of partial fractions. In each partition a suitable energy and angle parameterisation is chosen to factorise its soft and collinear divergences and subtract their singular parts via residue subtraction. The FKS method therefore defines its counter-terms not just at the level of the squared amplitude, but at the level of phase space measure times squared amplitude.
The CS method instead constructs its counter-terms purely at the level of the squared amplitude. The counter-terms are constructed by combining together soft and collinear limits which are promoted into the full phase-space by the so-called Catani-Seymour momentum mapping. This mapping also allows to analytically integrate the counter-terms over the singular emission phase space.
Due to more complicated overlapping divergences at the next-to-next-to-leading order (NNLO) neither the FKS nor CS methods can be naively generalised. The problem for FKS-like methods based on residue subtraction is that parameterisations which completely factorise all divergences present in a given partition do not appear to exist. The most successful approaches based on residue subtraction therefore make use of sector decomposition [11][12][13][14][15] to factorise the divergences. Such approaches, especially those based on sector decomposition, have been used successfully in calculations at the current state of the art; see, e.g., [16,17]. While these methods can be efficiently implemented numerically they come with their own set of disadvantages: they are parameterisation dependent and the integration of the counter-terms remains, for now, numerical.
Other approaches based on residue subtraction have been limited either to simpler applications, in which the divergences were factorisable [18][19][20], or were based on topol-JHEP08(2018)006 ogy dependent parameterisations [21], which is difficult to apply to more complicated final states.
Another class of subtraction methods used at NNLO are closer to the CS idea. A prominent such method is antenna subtraction [22][23][24]. Its counter-terms are based on physical matrix elements. Instead the colourful subtraction method relies on combining the universal soft and collinear limits [25] into suitable counter-terms. The advantage of the antenna method is that it leads to a comparably small number of counter-terms, whose analytic integration has been achieved. A disadvantage of the antenna method are that it is not fully local in the phase space, as certain spin correlations are ignored. This makes the phase space generation quite cumbersome. Despite this the method has been applied successfully in many state of the art NNLO calculations; see , e.g., [26,27]. In comparison the colourful subtraction method, see, e.g., [28][29][30], has so far been applied only to final state radiation in, e.g., [31]. While this method is fully local it comes with the disadvantage that the analytic integration of its counter-terms is highly non-trivial due to the appearance of Jacobians introduced by the mapping. It has thus been relying in part on numerical integration techniques for its counter-terms.
Slicing methods instead render divergent integrals finite by slicing out the singular regions. As in the subtraction method the integration over the singular region -which takes the role of the counter-term -is subsequently added back in (analytically or numerically) integrated form. The original slicing method, implemented at NLO, was based on imposing cuts on the Mandelstam variables [32,33]. These methods were never fully generalised to NNLO, although an extension was applied in mixed QCD-QED corrections in [34]. More recent developments at NNLO include kt-subtraction [35] and N-jettiness subtraction [36,37]. Both of these methods have been implemented in an impressive number of fully differential NNLO calculations; see, e.g., [38,39]. A clear advantage of these methods is the comparable ease of their implementation, since one simply implements a measurement function to cut out the singular parts of the phase space. While the ktsubtraction method has been applied primarily only to colorless final states, as well as final state radiation and massive partons for small multiplicities [40,41], it has the advantage that its counter-terms are relatively simple to integrate analytically. N-jettiness subtraction can be applied to more general processes; the integration of the required soft function is however challenging and requires a numerical implementation. An advantage of both the kt-subtraction and the N-jettiness methods is that the singular limits can be obtained from general factorisation theorems. A disadvantage of these methods is that the slicing parameters must be chosen small enough for the factorisation formula to be valid; and this may lead to numerical instabilities; see e.g. [42,43]. To address this problem one may need to add higher orders in the expansion around the slicing parameters. The challenge with this is that the structure of the sub-leading singular limits is more complicated; for instance derivatives of amplitudes will be required, which could be difficult to obtain for complicated processes.
Subtraction methods have also already been employed at next-to-NNLO (N 3 LO) in two incidences: an application of the projection to Born method [44] in DIS jet production [45] and a novel application [46] of the reverse unitarity approach [47] in Higgs production JHEP08(2018)006 (relying, for now, on the first two terms in an expansion around the threshold). While these are impressive calculations of so far unrivalled complexity, the methods they relied upon can not be (at least naively) applied for more general final states.
Despite the large number of existing methods, we propose -in the hope that it may overcome the shortcomings of existing methods -a new approach in this work. To accomplish maximally simple counter-terms, from the point of view of analytic integration, an FKS like residue subtraction procedure will be employed based on a Feynman diagram dependent slicing observable. By employing a slicing approach we can overcome the limitations of parameterisations which require singularities to be factorised.
By explicitely introducing a slicing parameter for each possible soft or collinear singular configuration it should be possible to promote this slicing scheme to a subtraction method by employing suitable phase space mappings not unlike in the CS, antenna or colourful subtraction methods. Indeed this idea has already been proposed in [48] in the context of the conventional NLO slicing method. Here we will only demonstrate this idea, for a simple example in section 3 in the context of the new slicing observable. The purpose of this paper is instead to establish the general principles of the method. In particular we work out the combinatorics of a generalised soft-collinear subtraction formula in section 4 -which uniquely defines a simple prescription for the soft and/or collinear final state counter-terms. How to extend the scheme to initial states is briefly sketched. We apply this framework for final state radiation at NLO and NNLO in pure Yang-Mills theory for arbitrary multiplicities in section 5 and perform the integration of all counter-terms required. We complete the section by reproducing the poles of the gluonic double real emission correction to the gluonic Higgs decay at NNLO. Possible extensions and future developments of the method are discussed in section 6.

Notation
In this section we introduce some of the notation used throughout the later sections. It is convenient to define the normalisation factor Of particular importance will be the D = 4 − 2 dimensional n-particle differential Lorentz invariant phase space measure: Here Q is taken to be an off-shell time-like vector, i.e., Q 2 > 0. Final state particles are constrained to be on-shell with positive energy by the distribution We mostly deal with massless vectors and masses occur in phase spaces only through the Mandelstam variables formed by squaring sums of momenta. For this purpose it is JHEP08(2018)006 convenient to introduce the shorthand Mandelstam variables are defined as follows: We thus mostly suppress the dependence on masses and employ the shorthand for all p i massless. The total phase space volume is defined as A massive sum of momenta, e.g. p 12 , with mass s 12 , is indicated by bracketed notation (12), e.g., dΦ (12)34...n (Q; s 12 , 0, . . . , 0) = dΦ (12) upon which much of this paper rests. Although this notation is intuitive and compact care has to be taken with identities such as p 1...n = p 1 + . . . + p n which, when substituted in eq. (2.7), could lead to appearances of δ (D) (0). Rather identites such as eq. (2.3) should be interpreted to arise in eq. (2.7) as a consequence of momentum conserving δ-functions. Similar considerations apply to the Mandelstam variables s ij .

Motivation
Before describing the general method in section 4 we will illustrate it here in the context of a simple example of a divergent phase space integral: This integral could appear in the real emission process γ * →q(1)q(3)g(2) at NLO in massless QCD. The integral is problematic in D = 4 due to soft and collinear singularities respectively located at 2 → 0 and 1||2, 2||3. It is instructive to see where the singularities are located in the space of Mandelstam variables. Let us express the phase space measure in terms of the variables s 12 , s 13 and s 23 :  where We thus see that the 3-particle phase space can be represented as the area constrained on the surface together with s ij ≥ 0. This physical region with the locations of its singularities, embedded in the three dimensional s ij -space, is shown in figure 1. We now wish to construct a subtraction scheme with which to isolate the finite part of the integral in eq. (3.1). Since there exists quite some freedom how to define such a finite part it is natural to ask: how can we define the divergent part in the simplest possible way? In other words, how can the evaluation of the divergent parts be maximally simplified. A scheme that accomplishes this for the UV divergences is the minimal subtraction (MS) scheme within dimensional regularisation. This is also what is commonly used to cancel the soft and collinear divergences between real, virtual and counter-term contributions in higher order calculations. The problem with the MS-scheme is however that it does not provide us with a prescription for defining a manifestly finite integrand which can be expanded around D = 4 before integration. This, for the practitioner of perturbative QFT, is indeed the dilemma of dimensional regularisation. Although it allows to define a unique finite part -it is not obvious how to obtain this finite part without performing the integral first in D-dimensions and then expanding the analytic functions around D = 4.
For the simple example the integration is easily carried out in terms of Γ-functions, whose expansion around D = 4 is trivial, to obtain:

JHEP08(2018)006
But when increasing the perturbative order very complicated integrals appear which evaluate to complicated hypergeometric series. For such functions a Laurent expansion around D = 4 may be difficult to obtain. Key is that one may not always be interested in integrating the integral over the entire phase-space. In fact for various reasons, such as detector efficiencies or large signal-swamping backgrounds, only part of the phase space may be experimentally accessible, or interesting. The phase space integration could then be constrained to an in principle arbitrary (infrared safe) region. It is therefore highly desirable to have a procedure to extract the finite part of an integral which does not rely on carrying out the integration. While, as summarised in the introduction, such procedures have been developed in the past, based on subtraction, sector decomposition and phase space slicing, we believe that none of these approaches fully matches the simplicity of the divergent parts defined via the MS-scheme. Inspired by algebro-geometric schemes based on blow-ups, which have been applied in the subtraction of UV-and IR-divergences of Euclidean loop integrals [49,50], we intend to provide a prescription which meets this criterion as closely as possible in the following.

Normal coordinates and phase space factorisation
We start by identifying a set of suitable variables -we call them normal coordinates or slicing parameters -with which to separate the phase space into singular and finite regions. Here we shall use the word normal in the sense in which it was introduced by Sterman, see, e.g., [51][52][53]; and we can identify these normal coordinates with the variables which we introduce to take soft and collinear limits. As a guide to find these variables we will use the phase space factorisation property of eq. (2.7) in terms of Mandelstams. It turns out that this property alone allows one to obtain suitable Lorentz invariant factorised soft and collinear limits of the phase space measure.
Let us first see how this works for the case of collinear divergences. A choice of a variable to parameterise the collinear limit C 12 is given by s 12 . The collinear limit is approached linearly as s 12 → 0. The factorisation property in eq. (2.7) then allows us to take this limit as follows: Since the 2-particle phase space measure dΦ 12 (s 12 ) can not be simplified further, the limit operation acts only on the remaining phase space measure dΦ (12)3 , which has support in the limit s 12 → 0. Here it is useful to introduce the Sudakov parameterisation such that we can parameterise p 1 = z 1 p 12 + s 12 z 2 2p 12 .n n + √ s 12 z 1 z 2 e ⊥ , p 2 = z 2 p 12 + s 12 z 1 2p 12 .n n − √ s 12 z 1 z 2 e ⊥ . with z 1 + z 2 = 1 and e ⊥ being a space-like unit length (|e ⊥ | = 1) vector transverse to both p 12 and n. This allows us, at the expense of the massless reference vector n, to control how the off-shell vector p 12 present in dΦ (12)3 approaches the massless vector p 12 lim s 12 →0 We thus obtain the following factorisation of the phase space measure in the collinear limit: with the collinear phase space defined as What is characteristic to this factorised limit is that the remaining Q dependence is present only in the reduced measure dΦ 123 (Q). In contrast the s 12 -dependence is present only in dΦ 12 (s 12 ). This in effect means that the variable s 12 is no longer bounded from above. Since the factorisation is only valid in the limit of small s 12 it is sensible to introduce "by hand" a small upper cutoff for s 12 to make the integral over this measure well defined. We will describe below how to do this in a consistent manner. We now come to the parameterisation of the soft limit. A Lorentz invariant variable suitable to parameterise the soft limit p 2 → 0 is (3.12) This variable, which linearly approaches zero as p 2 → 0, is also directly proportional to the energy of p 2 in the rest frame of p 13 . Due to the relation eq. (3.4) we can also identify the limit s 2(13) → 0 with the limit s 13 → Q 2 . This allows us to derive the soft phase space factorisation from the factorisation property eq. (2.7), as follows: Only the term dΦ 13 (s 13 ) has further support in this limit, and we find: with the soft phase space measure defined as: where we have used ds 13 = ds 2 (13) . Note that the soft measure depends on the hard momenta p 1 and p 3 . Thus even after integration over the soft momentum the soft limit retains a dependence on the variable s 13 . In contrast the collinear limit factorises entirely and retains no dependence on any hard momenta.

Geometry of regions
So far we have not discussed how to construct a subtraction method which consistently combines the soft and collinear limits we introduced in the previous subsection. A way to attack this problem is to use a phase space slicing approach. The idea here is that the phase space can be separated into a finite region (F ) and singular, that is soft (S 2 ) and/or collinear (C 12 , C 13 ), regions. To accomplish this decomposition we will associate a set of small dimensionless parameters a i for S i and b ij for C ij to bound the slicing parameters of each singular region from above. This procedure will naturally lead to a classification of the overlap of soft and collinear regions. We will associate: S 2 : {s 2(13) < a 2 s 13 }, Here we have used s 13 as the scale entering the upper bound of the soft variable s 2 (13) . First let us remark that for small a 2 and therefore also s 2(13) we have s 13 ∼ Q 2 , and so we could have equally well used Q 2 here. However s 13 is the natural hard scale appearing in the soft phase space as already commented on above, and as we will see later it will turn out important in more complicated situations to use this s 13 instead of Q 2 . For this reason we introduce this notion already here. The decomposition of the phase space into its different finite and singular regions is best understood geometrically and is illustrated in figure 2. In the figure it is shown how the regions can be chosen such that the overlap of the two collinear regions C 12 and C 23 is contained inside the soft region S 2 . Using this feature and the additivity of areas we arrive JHEP08(2018)006 at the following partition of unity: (3.16) with the different regions defined by:

Counter-term integration
Let us now come to the evaluation of the integrals corresponding to the singular regions. We will start with the soft. A convenient parameterisation of the soft phase space measure is given by which allows us to obtain We continue with the evaluation of the collinear region. A convenient parameterisation for the collinear region is given by: The collinear limit of eq. (3.1) then leads us to the following integral: We continue with the evaluation of the soft-collinear overlap contribution. We can simplify the calculation of the overlap contribution by demanding that b 12 a 2 . In the limit of small b 12 , which in turn forces a small s 12 , the soft region then simplifies to lim s 12 →0 Θ(s 12 + s 23 < a 2 s 13 ) = Θ(z 2 s 123 < a 2 s 123 ) = Θ(z 2 < a 2 ) , (3.22) where we also used that in the soft region z 1 ∼ 1. In other words the soft-collinear limit, in this setting, is conveniently computed by taking the limit z 2 → 0 in the collinear phase space. The soft-collinear phase space measure is thus given by (3.23)

Slicing method
We will now test the finite part of the subtraction terms numerically using the slicing method. The advantage of the slicing approach is the simplicity with which the finite part, defined by can be implemented in a numerical simulation in D = 4. To implement the hierarchy between soft and collinear limits we introduce a single slicing parameter λ, such that We can then define the quantity measures the percent difference by which this quantity differs from zero for a small finite value of λ. While we know both I Singular (Q) and I(Q) analytically we can compute the finite integral I F (Q) numerically; this is plotted over a range of values of λ in figure 3. The figure on the left clearly shows the formation of a plateau in the range 10 −3 < λ < 10 −6 . As is expected for a slicing scheme the numerical accuracy deteriorates as λ is decreased and improves for larger values of λ, where however the counter-terms can not be used as a reliable approximation, and power corrections in λ would be required. It is thus evident, given the simplicity of the example, that this approach gives rise to a rather poor numerical accuracy. Even after sampling 10 8 points using the Cuba implementation of Vegas [54] (albeit in a non-optimal phase space parameterisation) we are not able to arrive at an accuracy much better than 1% for the most optimal values of λ.

Subtraction method
An alternative to the slicing approach, presented in section 3.4, is to rewrite the region approximants as local counter-terms. In this subsection we will illustrate -for the simple example -how a slicing method can be promoted to a fully local subtraction method. This idea is not new and was already discussed in, e.g. [37,48]. The promotion can be accomplished in several different ways. One method would be to employ a momentum map, of the kind which has been employed by Catani  in the dipole subtraction method [6]. Alternatively one can use, what we shall call, the re-weighting approach, this mimics in some sense how the limit subtraction is embedded in the full phase space in the FKS subtraction method. In our approach we can use a similar idea, based on the particular phase space factorisation property used to derive a certain limit.
In the context of the simple example the subtraction method can be implemented very easily by noting that we can match soft and collinear measures with the full phase space measure by multiplying with the corresponding factors which were simplified to unity in the limit taking procedure. This leads us to the following relations: where we have made explicit choices for the reference vectors n in the different collinear limits. Such that the momentum fractions z ij become: Using these relations for the measures we can derive the following alternative representation for the finite part defined in eq. (3.27) in D = 4:

JHEP08(2018)006
A numerical evaluation of this finite part for a range of suitable values of is presented in figure 4 using again the Cuba Vegas implementation. The clear advantage of the subtraction method is that it allows to use arbitrary values for λ ∈ (0, 1], since the integrated counter-terms evaluate by construction to those used in the slicing method for any value of a i and b ij .
In contrast to the slicing method better convergence is observed for large values of λ. In fact it appears that λ = 1, which corresponds to the counter-terms ranging over the entire phase space, is the optimal choice for this example. Notably the accuracy reached for this value of λ is 100 times as good as that reached by the slicing method using the same number of numerical evaluations. The subtraction method -to which we have promoted the slicing method -therefore appears far superior in terms of numerical accuracy and stability when compared to its parent slicing method.

Normal coordinates and phase space factorisation
In the previous section we applied the geometric subtraction method in a simple example. Let us briefly summarise the idea behind the procedure. We introduced the variables s ij = 2p i .p j to define a collinear region and the variable s i(jk) = 2p i .p jk to define a soft region. With these variables we then derived suitably factorised limits of the phase space measure and furthermore partitioned the phase space volume into finite and singular (soft and/or collinear) regions.
This procedure can be generalised to arbitrary perturbative order. For instance at NNLO we can use the variable s ijk to subtract the triple collinear limit i||j||k and we can use the variable to subtract the double soft limit ij → 0 sensitive to the hard momenta k and l. The phase space factorisation property can be used to determine the corresponding factorised phase space volume limits, with which to integrate the singular limits of amplitudes. This leads us to the following phase space factorisation in the triple collinear limit: with the triple collinear phase space measure defined by

JHEP08(2018)006
To parameterise the momenta in this collinear limit we use, as before, the Sudakov parameterisation (expressed in terms of Mandelstams, rather than transverse momenta): As before e ⊥ i are space-like unit length (|e ⊥ i | = 1) vectors transverse to both p ijk and the reference vector n. The off-shell vector p ijk approaches the massless vector p ijk in the limit of vanishing s ijk : lim In the double soft limit we obtain the phase space factorisation: with the double soft phase space measure given by: The pattern of these measures follows those defined at NLO. However there is subtle difference between the soft and double soft measures. The double soft measure is not simply dΦ ij(kl) , since there exist further support in the limit ij → 0. Instead an explicit form for it is given by: The phase space measures at yet higher order, e.g. at N 3 LO, can be defined similarly. For the m-collinear limit we would use the variable with the following phase space factorisation and m-collinear phase space measure

JHEP08(2018)006
Similarly we may define the m-soft variable 14) and the m-soft phase space measure dΦ (k,l)

Soft and collinear forests
From the conceptual side it remains to determine the overlap contributions. Since the phase space volume of the example given in section 3 was simple enough, we were able to construct the partition based on the geometry of the phase space in Mandelstam variables. For higher dimensional phase spaces it will be advantageous to have at hand a formalism which allows to derive this partition in a more algebraic manner. We will employ the properties of Heaviside step functions for this purpose. Employing the normal coordinates defined above we associate Θ-functions for each region as follows: Here we assume that all soft regions are defined in some rest frame of p kl . This rest frame could be chosen differently for different soft divergences. We will focus our discussion more on this point in section 5, where we show how different eikonal factors can have their regions bounded in different rest frames. Our starting point for now is again the equation: which follows given that finite (F ) and singular regions are to cover the entire phase space volume. Let us now define the set R as the set of all possible singular regions, excluding their overlaps, such that for the example of eq. (3.1) we would have R = {C 12 , C 23 , S 2 }. We can then write which combined with eq. (4.18) leads to Here we sum over all nonempty subsets U , each of size |U |, of the set R. This expression still lacks knowledge of the geometric structure of the soft and collinear regions as well as the perturbative order.

JHEP08(2018)006
To get a first feel for this equation let us study its consequences using R = {C 12 , C 23 , S 2 } as input. We shall use the notation Θ(A)Θ(B) = Θ(A ∩ B) if regions A and B depend on common momenta. Using this notation we then find: which agrees with eq. (3.16), if we apply the relation Indeed this relation follows from the geometric construction introduced in figure 2, since the soft region contains the intersection C 12 ∩C 23 . By demanding the hierarchy a i b ij we can thus guarantee its validity. But eq. (4.22) would not hold for other parameter choices such as a i b ij . This shows how, in a simple example, the geometry of regions plays an important role.
Let us continue with our exploration of eq. (4.20). To be able to include more complicated final states into our discussion let us introduce the measurement function J (l) 1...n+l , with l denoting the maximum number of unresolved partons, which are permitted by this measurement function. In particular the measurement function obeys the relations: (4.23) Here the notation J ... i... indicates that J does no longer depend on the soft momentum i, while J ... ij... indicates that the collinear momenta i and j have been merged into the massless momentum ij. The (purely) real emission contribution to an arbitrary N l LO observable may then be defined as with M 1...n+l an n + l parton tree-level amplitude.
In the following we wish to define a set U (l) , whose elements are themselves sets of possible singular regions which may pass the criteria of the measurement function J (l) , and where we omit all those overlap regions which cancel by vitue of eq. (4.22) and its generalisations at higher order. We can thus write where now cancellations among different regions have taken place and only those regions are included which can pass the criteria of the measurement function. Let us start by studying the case l = 1. The only divergences allowed by J 1...n+1 are collinear i||j and soft i → 0 or their overlap with common partons. Let us recall at this point the slight mismatch between our region definitions C ij and S i and the locations of singularities i||j and i → 0 respectively. Since C ij is the region defined by s ij < b ij Q 2 JHEP08(2018)006 even a soft singularity can fake the region C ij ; an apparent paradox which is resolved by subtracting the overlap contribution C ij ∩ S i . It follows that care must be taken when considering the possible regions which may pass the criteria of the measurement function. While the NLO measurement function J (1) may not allow for a singularity i||j and j||k which would correspond to a triple collinear singularity, the region C ij ∩ C jk can still be mimicked by a soft singularity j → 0. As before we escape this apparent new region by relying on the cancellation in eq. (4.22), as in the simple example. At NLO we thus define: (4.26) where the notation {C ij } is not meant to indicate the set of all collinear divergences, but rather one set for each collinear region. Given for example the regions R = {C 12 , C 23 , S 2 } we would then have Let us now come to the definition of U (2) , the set of all possible singular regions which pass the criteria of the NNLO measurement function J (2) . To define U (2) more precisely we must establish what the geometric cancellation identities are. Since these identities depend on the geometric properties of the soft and collinear regions, they depend on the parameters a i... and b ij... . To make progress we choose the following order as it produces simple iterated phase space volumes for the counter-terms. Further this ordering gives rise to the following region cancellation identities: and and is therefore consistent with eq. (4.27). Here A, B and D are sums of momenta, while i and j are single (massless) momenta. We do not claim that these are all the identities which are fulfilled, but they suffice to establish the desired region cancellations at NNLO. Pictorially these identities are illustrated in figures 5 and 6. A derivation is sketched in appendix A. The region cancellation identities are useful since they considerably reduce the number of different singular regions, or equivalently counter-terms, which are required in the subtraction. At NLO they allow us to obtain  from eq. (4.30), since terms containing C ijk are rejected by the NLO measurement function. At NNLO we obtain instead the following set: But furthermore these regions allow us to identify the set U (2) as the Cartesian product of a set of soft forests U S times a set of collinear forests U C modulo measurement function constraints. We write this statement as with U

JHEP08(2018)006
and U (2) The sets U (2) S and U (2) C can be thought of as the soft and collinear analogues of Zimmermann's set of forests of UV-divergent subgraphs [55]. This should not come as a complete surprise, since it is well known that the subtraction of both soft and collinear divergences follows similar patterns as those found in the UV; see for instance [56][57][58][59] and references therein. Nevertheless it is not completely trivial that the regions defined with the particular ordering of eq. (4.27) follow this pattern as well. One may associate it to the observation that also in the context of subtracting UV-divergences an order is chosen with which to subtract singular subgraphs; such that the smaller subgraphs are subtracted before the larger. Similarly we may have chosen the relative ordering of soft and collinear divergences according to such a pattern. Regardless of this similarity eqs. (4.28), (4.29) and (4.30) are highly dependent on the relative ordering of soft and collinear regions, which has no analogue in the subtraction of UV divergences alone.
Multiplying out the product in eq. (4.34) and discarding those regions which can not pass the criteria of the measurement function we arrive at: (4.37) The size of this list shows the enormous increase of complexity which is encountered at NNLO, when compared to NLO. Before moving on to the definitions of the asymptotic phase space measures in the various regions let us conclude this subsection with the conjecture that eq. (4.34) is valid also for the case of l potentially unresolved emissions at N l LO: as long as the order is chosen, with U

Asymptotic phase space measures
Having fixed the order of limits we are now in a position to determine the asymptotic measures associated to the singular regions. To compute the asymptotic measure associated to a particular region U one should take the following sequence of limits of the expression: (4.42) there are two different cases, of which the more complicated one is: There exists a subtlety when taking the limit a ij → 0 in the third line, which forces the single soft limit j → 0, since the argument of the factor Θ(a i s jk > s i(jk) ) does JHEP08(2018)006 not simplify in this limit. To understand this better let us study the scaling with a ij which the different factors portray: The two terms on the right hand side are thus (at least naively) not of the same size; their exist however no support to further simplify this expression. The only plausible interpretation is that the upper bound on the left forces both Mandelstams on the right, i.e. s ij and s ik to scale as ∼ a i a ij . Given that both i and j are sufficiently soft it thus appears that i is forced to become collinear to k for s ik to be of similarly size as s ij .
A different way to come to a similar conclusion comes from the result of the integral: where S (j,k) i corresponds to the singular limit of an amplitude. It is thus clear that this integral is homogeneous in p j , although the constraint eq. (4.46) does not appear to be so.
We will leave it as an exercise for the reader to work out the asymptotic forms of the measures corresponding to other regions. Compact expressions for the resulting measures for all required regions can be found in appendix C. In particular we find that at NLO and NNLO all regions can be written using the following primitive phase space measures: It is left to future work to show that an analogous statement will hold also at higher orders, i.e. that one requires a certain set of new primitive measures corresponding to, e.g., {C ijkl }, {S ijk } and {C ijkl , S ijk } at N 3 LO, while all other regions will factorise into the measures already present at lower orders.

Soft and collinear master integrals at NNLO
The integration of singular limits of amplitudes over the primitive NNLO phase spaces ({C ijk }, {S ij } and {C ijk , S ij }) is more complicated than the integration over NLO primitive phase spaces or their iterated limits which occur at NNLO. While the latter evaluate to Γ-functions, the primitive NNLO limits can lead to hypergeometric functions of type p F q .

JHEP08(2018)006
In order to minimise the total number of integrals it is advantageous to express them in terms of a basis of master integrals. This basis can be constructed by solving the relevant IBP identities using the Laporta algorithm. Using reverse unitarity [47,60] to treat Dirac δ-functions which appear in the definition of their associated phase space measures, it is a well defined task to use public codes to set up the IBP reduction for these limits. How to use this technique in the presence of Heaviside step functions is briefly reviewed in appendix B. To accomplish this task we have employed the public softwares FIRE [61,62] and AIR [63].
For the double soft this procedure leads to exactly two master integrals, they can be easily extracted from the two independent double soft Master integrals which were presented and evaluated in [64], where they appeared in the threshold limit of the Higgs boson cross section in gluon fusion. In our conventions this leads to: The double soft-triple collinear master integrals are in one-to-one correspondence to those in the double soft limit, and results for these slightly simpler integrals can be read off from their double soft counter parts using the phase space parameterisation which was presented in section 5.2 of [65]. In the triple collinear limit we will require four master integrals. These, among two other simpler ones, were presented and evaluated in [66], where they appeared in the context of jet fragmentation. In our conventions these results can be expressed as follows The expansions around = 0 of the hypergeometric functions, can be obtained with the HypExp package [67].

Example at NNLO
Let us now consider an example where we can apply the ideas we developed in the last section in a simple setting. The example we shall consider is where we have set Q 2 = 1 and We have analytically evaluated this integral applying again reverse unitarity, IBP reduction as well as the results for the master integrals given in [68]. This integral contains the following set of singular regions: • {C 34 , S 34 }: Summing the counter-terms, weighted with the appropriate signs, we obtain: where we use the notation L n z = log n z . (4.76) A numerical evaluation of I 1F using the slicing method for a range of values of a parameter λ, which fixes the cut-off parameters via is shown in figure 7. The strict hierarchy which the slicing cut offs must satisfy unfortunately limits the range of possible values of λ which can be chosen without risking loss of numerical stability. Nevertheless good numerical convergence is observed in the range λ ∈ [0.1, 0.001] for this particular example. In general we may not expect such good convergence in the range λ ∼ 0.1 which is likely due to the trivial numerator.

Generalisation to initial states
Let us briefly comment here also on the generalisation of the geometric subtraction procedure to divergences in the inital states. Since soft divergences are always associated to momenta in the final state, no new complications are encountered here in the initial state. Indeed the double-soft counter-terms we propose here have already been used in initial state calculations before, see [69].
The new feature is that there can be collinear divergences, where final state partons become collinear to partons in the initial state. One can define the slicing parameters for these regions analogously to those in the final states, as the Mandelstams s ij = −2p i .p j , where now i (j) is in the initial (final) state. Phase space factorisation in these spacelike Mandelstams similarly provides one with a way to define suitable singular phase space volumes. Let us briefly sketch how this would work for a simple NLO process:

JHEP08(2018)006
with H a higgs boson, or a collection of "hard momenta", p 1 , p 2 the initial state partons and p 3 a final state parton. A collinear divergence can arise in the limit p 3 ||p 1 . To define a counter-term for this region we define the collinear region s 13 < −b 13 Q 2 . Phase space factorisation for the counter-term yields: Let us now introduce a parameterisation for the collinear phase space dΦ 1→(13)3 . We parameterise the momentum as: such that Then we can derive the following factorised limit: The massless collinear momentum p 1 = z 3 p 1 , which enters the Born amplitude in the collinear limit, is thus modified only by the longituidinal momentum fraction z 3 . This feature indeeds leads to the typical convolution integral, which one expects from initial state collinear divergences. Further more our prescriptions for dealing with the soft collinear overlap will lead to subtraction in the variable 1−z 3 with an upper bound a 3 . For the value 1 a 3 = 1 this is indeed identical to the usual plus prescription with which these divergences are commonly subtracted, see e.g. [5]. The geometric subtraction procedure can thus be extended in a straight forward manner to the subtraction of collinear divergences in the initial states. In fact the integrated triple collinear counter-terms at NNLO have already been presented in [66]. We will further explore this extension in the future.

Counter-terms for final state real emissions in Yang Mills theory
In the following section we will employ the geometric subtraction formalism to construct suitable counter-terms for tree-level Yang Mills amplitudes; that is QCD amplitudes without quarks, which we ignore here for simplicity. While amplitudes factorise completely in collinear limits, soft limits are color correlated. To make the counter-terms as simple as possible we will employ tailor-made soft volumes for the different color correlated eikonal factors, which make up a particular soft limit.

JHEP08(2018)006
To accomplish this task we correlate the sum over singular regions with individual interference terms contributing to the squared amplitude: where the sum over k, m labels different color projected (sets of) Feynman diagrams contributing to the matrix element M, such that M = k M k . We can then use eq. (4.25) individually for each component of Θ(Singular). In order to ensure gauge invariance in the counter-terms it is sufficient that sets of Feynman diagrams multiplying a particular entry of Θ(Singular) conspire to singular limits which are gauge invariant, such as the Altarelli-Parisi splitting function or the eikonal factor. In the following we shall present how to define Θ(Singular) explicitly at NLO and at NNLO. Further more we will provide all integrated counter-terms for final state emissions at these first two orders.

Counter-terms at NLO
At NLO we define the sum over singular regions as with the set of regions defined by To define a suitable Θ-matrix it is sufficient to define how it behaves in the singular limit.
In the soft limit we require it to behave as follows: where the eikonal factor is given by Here T i denote the (by now standard) color charge operator; see, e.g., [25]. In the collinear limit we require the Θ-matrix to factorise completely:

JHEP08(2018)006
Here (P ij ) µ 1 µ 2 is the standard spin correlated gluonic Altarelli-Parisi splitting function defined in, e.g., eq. (12) of [25] and |M µ 1 µ 2 ... ij... | 2 is the spin correlated squared matrix element, denoted T µν in [25]. Let us remark also that the different soft volumes Θ(a n+1 s ij − s n+1(ij) ) all collapse to the same limit in the soft-collinear limit; since To write down the integrated counter-term it is convenient to define the functions: as well as the following linear combination: It is straight forward to show that eq. (5.14) agrees with the corresponding one-loop pole operator given by Catani in, e.g., [70].

Counter-terms at NNLO
At NNLO we define the sum over singular regions similarly as with U (2) defined similarly, although not identically due to the more elaborate soft structure, as in eq. (4.37). To define the limits Θ(C ij ) and Θ(S i ) we simply use the NLO definitions. The triple collinear limit is defined similarly to the double collinear:

JHEP08(2018)006
with (P ijk ) µ 1 µ 2 the triple collinear Altarelli-Parisi splitting function, defined, e.g., in eq. (66) of [25]. The double soft limit receives two independent color-correlated contributions: into the color off-diagonal terms. The double soft limit thus contains two terms. The first term factorises over the soft momenta and contains color-kinematic correlations with up to four hard partons (Wilson lines). Instead the second term contains kinematic correlations of the two soft momenta and color-kinematic correlations with up to two hard partons. It is interesting to note that while both terms are singular in triple collinear regions i||l||k and j||l||k only the second term contributes to the limit containing k||l and only the first term contributes to limits containing i||k, i||l, j||k and j||l.
A natural measure for the second term is dΦ S kl (s ij , a kl ) which we introduced earlier. Instead we shall treat each of the eikonal factors in the first term with a single soft phase space measure, i.e. with dΦ S l (s rt , a l ). The "true" double soft measure will thus be associated only to the second term, while the first term is naturally associated to the product of two single soft limits. These considerations lead us to define: and lim a kl →0 lim (a k ,a l )→0 With this distribution of the theta functions it follows that the double soft limit k, l → 0 is not entirely controlled by the limit a kl → 0, instead also a k , a l → 0 is required for both terms JHEP08(2018)006 in eq. (5.18) to diverge. The region cancellation between the regions Θ(S ij )Θ(S i )Θ(S j ) and Θ(S i )Θ(S j ), which was given in eq. (4.28), therefore only holds for the second term in eq. (5.18), since the first does not contribute to the region Θ(S ij ). Let us now consider what happens in the strongly ordered double soft limits corresponding to {S ij , S i }. One can show, by taking the successive soft limits, that this limit becomes: Thus the different single eikonal factors which contribute to the strongly ordered limit of the non-Abelian double soft limit come with their distinct single soft phase spaces. A caveat of the method is that in the strongly ordered soft limit certain collinear limits such as {S kl , S k , C il }, which would usually not survive in the non-Abelian double soft factor, e.g. {S kl , C ik } is not singular, now survive: The re-scaling invariance of the eikonal factor, ensures that the last two terms in eq. (5.24) would cancel, if it was not for the differing Θ-functions which break the re-scaling invariance upon which the cancellation mechanism relies. The chosen distribution of single and double soft Θ-functions similarly splits the various overlapping soft-collinear limits. For instance the triple collinear double soft limit splits into a non-Abelian part: and an Abelian part:

JHEP08(2018)006
where n denotes the collinear reference vector such that, e.g., The single soft triple collinear limits instead are split into three different eikonal factors: It is interesting to note how the different eikonals contributing to this limit come with their individual phase space volume constraints. The first term can be interpreted as a degeneration of the limit {C ijk , S k } into something like an ordered {C ij , S k } with the soft scale smaller then the collinear one. The second term on the other hand takes the form of a soft collinear limit where the soft and collinear scales are of the same order. Let us briefly analyse the kinematics of this second limit. Taking the limit a k → 0 in Θ(a k z i − z k ) forces the momentum fraction z k to vanish; in turn this means that and so s ijk → s ij +s k ij . It appears almost as something of a miracle that the triple collinear splitting function produces a factor (s ij + s k ij ) 2 in the numerator which precisely cancels the overall denominator 1/s 2 ijk . A feature which leads to a welcome simplification for the integrated counter-term associated to this limit. The strongly ordered double soft limit receives contributions only from the non-Abelian double soft limit of the triple collinear region and yields Other limits can be worked out similarly starting from these expressions and using the soft and collinear limits of amplitudes. All in all, with this choice of the single and double soft color correlated phase space boundaries, the NNLO set of regions which enters eq. (5.16) is given by: It is convenient to re-organise the sum over regions in eq. (5.16) by introducing the subdivergence subtracted regionsC ijk ,Ŝ ij , andC ij , as certain subsets of U (2) . The region C ijk is defined as the set of all regions which contain C ijk . The regionŜ ij is defined as the set of all regions containing S ij apart from those also containing C ijk ; and the regionC ij includes C ij and its overlaps with the regions S i and S j . Using Θs we can also define these regions as follows: In terms of these combinations we can then write out the sum in eq. (5.16) as follows: An explicit representation for the pole part in terms of the different region approximants can then be written as follows: ..n+2 IŜ g k g l (s ij , a kl , a k , a l , t kl , t ik , t jk , t il , t jl ) Let us remark here that the poles of the observable O 1...n+2 do of course not depend on the parameters a i... and b ij... . To get a simpler expression independent of these parameters one can alternatively set all the parameters to unity, i.e., b ij... = 1, a i... = 1. However to explicitly verify the cancellation of these parameters in the pole parts constitutes a welcome cross-check for its validity for a given process.
Given the results of all the integrated counter-terms, for which we present simple integral representations in appendix C, one can assemble the functions IC g i g j g k and IŜ g i g j which make up the basic new building blocks needed at NNLO to construct the integrated counter-terms for arbitrary multiplicities. We give these functions as expansions in including terms up to O( 0 ); although the results of this paper allow to construct them to arbitrary order if needed. Since these functions are lengthy, due to the many parameters JHEP08(2018)006 on which they depend, we provide them in computer readable format with this publication. We nevertheless provide them here for the useful case where one sets for all i, j, k: and

The poles for the H → gggg phase space integral
A simple example which allows us to test the validity of eq. (5.41) is given by the quantity We consider the corresponding amplitude in the heavy quark effective theory where the Higgs boson couples to gluons directly via the effective Lagrangian: with C eff a Wilson coefficient, H the Higgs boson field and G µν a the gluon field strength tensor. We can evaluate the inclusive quantity O H→g 1 g 2 g 3 g 4 using IBP reduction and the JHEP08(2018)006 master integrals presented in [68] to get (setting Q 2 = 1): Using eq. (5.41) we can confirm the pole parts of this expression independently to obtain (using again the parameter choice of eq. (5.42)): The α i... and β ij... parameters thus cancel in the poles and, more importantly, reproduce the correct result. A proper integrand-level implementation of these counter-terms must therefore numerically cancel the finite log-dependent parts which make up the finite part of the integrated counter-term O Singular H→g 1 g 2 g 3 g 4 . Unfortunately it is infeasible to check the finite part using the slicing approach, since the hierarchies which we have assumed for the four different slicing parameters are too large to simulate numerically. To guarantee an accuracy of 1% one would have to take the smallest parameter to be around 10 −8 . Even with enormous computer resources this would prove to be challenging. Further more one would have to find an explicit decomposition of the squared amplitude which respects the distribution of θ-functions which we have fixed in the various soft limits.

JHEP08(2018)006
In order to accomplish this task we will therefore have to resort to the subtraction method to which one should be able to promote the slicing method. In the subtraction method the hierarchy between slicing parameters can be relaxed and an accurate stable numerical evaluation should be feasible. Since it requires working out suitable phase space parameterisations or momentum mappings for all the 26 different counter-terms this is a however a formidable task in itself and will be left therefore to future work.

Conclusions
In this work we introduced a new scheme for the subtraction of IR divergences in real radiation phase space integrals, which is based on a particular Feynman diagram dependent slicing observable.
To construct this observable we introduced slicing parameters or normal coordinates for each soft or collinear singular region. We then exploited this freedom by demanding a particular strict hierarchy in the size of these parameters. This hierarchy lead to simple counter-terms from the point of view of analytic integration but it also leads to numerical difficulties as a slicing scheme. By promoting the slicing scheme to a subtraction scheme the hierarchy of the slicing parameters can subsequently be dropped and any physical value of the slicing parameters can be used; in this way good numerical convergence can be obtained. Here this was demonstrated explicitely for a simple NLO example, but not yet at NNLO.
Based on the geometric properties of the observable we established a subtraction formula which summarises the combinatorics of the various counter-terms for single and double real emissions and conjecture its general form for an arbitrary number of unresolved emissions.
We applied the formalism to final state real radiation at NLO and NNLO in Yang Mills theory and integrated all the required counter-terms. We employed reverse unitarity and IBP reduction to simplify the calculation of the most complicated counter-terms. We showed in particular how the master integrals required for these counter-terms can be extracted from existing calculations of unrelated quantities. We were thereby able to compute or extract all required counter-terms in terms of Γ and p F q hypergeometric functions to all orders in the dimensional regulator.
We tested the integrated counter-terms by reproducing the poles of the purely double real emission contribution to the gluonic Higgs decay in the heavy quark effective theory.
There exist many possible directions to extend this work in the future. The most important step will be to show that the scheme can indeed be employed to build local counter-terms at the level of the integrand. The next logical step would be to extend the scheme to include also initial states; this would open up a new path for computing LHC observables. Another step is to extend the scheme to real-virtual corrections; one can foresee that this should be a straight forward application of the techniques presented here for the case of real radiation at NLO. Beyond one can also imagine to use the scheme for N 3 LO calculations and/or to include massive quarks into the formalism.

JHEP08(2018)006
Given the simplicity of the integrated counter-terms, the potential locality of the counter-terms and their well defined combinatorial properties, we believe that the proposed scheme may well become an important method for performing higher order calculations in perturbative QCD in the future.

Acknowledgments
This work has been supported by the European Research Council (ERC) advanced grant 320389 and the NWO Vidi grant 680-47-551 "Decoding Singularities of Feynman graphs". I am indebted to S. Badger and T. Melia with whom I collaborated on the project during its early stages and who were influential in establishing some of the basic ideas. I would also like to thank them for their continued enthusiasm and the countless conversations we had on the topic since then. I would further like to thank M. Borinsky, H. Kissler, D. Kreimer, E. Laenen, and W. Waalewijn for useful discussions and E. Laenen, T. Melia and V. Del Duca for comments on the manuscript.

A Region cancellations
In the following we will derive several identities among different sets of overlapping regions.
• Let us start with the identity lim A||B,A||D where A, B and D are non-intersecting sets of momenta. We consider here the special case where the region C AB ∩C AD corresponds to a collinear momentum configuration only. It is sufficient to show that s ABD ≤ b ABD given that s AB ≤ b AB and s AD ≤ b AD . Now since A||B and A||D it follows that A||B||D, which in turn implies that which is in accord with the ordering given in eq. (4.39) and guarantees eq. (4.29).
• We proceed with the identity lim A,B→0 where again A and B are two non-intersecting sets of momenta. In order to prove this identity we have to specify the hard momenta which enter the constraint of the soft slicing parameter. Two different choices will be relevant for us. The derivation is easiest in the case when the hard momenta of the different slicing parameters are chosen identically as say k and l, such that which is automatically satisfied given that a A a AB . Taking the limits in the other order leads to a similar conclusion and we conclude that (S A ∩ S B ) ⊂ S AB . However there exists another case of interest where the hard momenta appearing in the different soft order parameters are not identical but are nested in the following sense: Again we obtain different results for S AB depending on whether we take first the limit in a A or a B . If for instance we take first B → 0 the conclusions are as in the case before. If however we first take the limit A → 0 we have to ask whether s B(kl) < a AB s kl is guaranteed by s B(il) < a B s il . For this purpose it is useful to write where we have used that s i(kl) is much smaller than s kl and y B;kl ∈ (0, 1]. It thus follows that s B(kl) a AB s kl , given a B a AB . Thus for the cases of interest eq. (4.28) is fulfilled, given the ordering of eq. (4.39).
Using similar arguments one can show that a more general identity is also true.
• Let us now consider the following 4-term cancellation identity:

JHEP08(2018)006
where A is a set of momenta not containing the single momenta i and j. The identity derives from the fact that the overlap region C Ai ∩ C Bi contains singular momentum configurations of two types: or their overlap. While for the second case A||i||j we can relie on the identity given in eq. (4.29), we must show that it holds also for the case A → 0. To accomplish this it is sufficient to show that since the other two terms in eq. (A.15) are contained in a sub-region of these regions and must thus cancel by the same mechanism. The soft region is given by the bound Using the constraints s Ai < b Ai and s Aj < b Aj we find Thus, assuming b Ai ∼ b Aj we must fulfil the bound: Now since the momenta i and j are not allowed to be collinear to the momenta in A, we can write this as b Ai ≤ a A b Aij 2 .
(A. 20) This bound corresponds to the worst scenario, since it may be that A||i||j may not be allowed by the measurement function; nevertheless this inequality is consistent with the ordering suggested in eq. (4.39).

B Reverse unitarity with Heaviside step functions
In this section we will briefly review the method of reverse unitarity and show how it can be used to find IBP relations among the integrated counter-terms. The central idea of reverse unitarity is always to treat Dirac delta functions as cut propagators by making the replacement: Using reverse unitarity one can then find relations among integrals of the kind I(s), whose integrands |M| 2 can be related by IBP relations. Note that since here we are interested in the limit s → 0 the final integral over the variable s is always trivial. This would not be the case away from the limit. Similarly one can repeat the trick above for several step-functions.

C Integrated counter-terms
In the following we provide expressions for the counter-terms associated to all the required regions.