Factorisation and Subtraction beyond NLO

We provide a general method to construct local infrared subtraction counterterms for unresolved radiative contributions to differential cross sections, to any order in perturbation theory. We start from the factorised structure of virtual corrections to scattering amplitudes, where soft and collinear divergences are organised in gauge-invariant matrix elements of fields and Wilson lines, and we define radiative eikonal form factors and jet functions which are fully differential in the radiation phase space, and can be shown to cancel virtual poles upon integration by using completeness relations and general theorems on the cancellation of infrared singularities. Our method reproduces known results at NLO and NNLO, and yields substantial simplifications in the organisation of the subtraction procedure, which will help in the construction of efficient subtraction algorithms at higher orders.


Introduction
Infrared divergences arising from exchanges of soft and collinear massless particles are well known to cancel in infrared-safe observable cross sections, where singularities in virtual corrections to scattering amplitudes are compensated by divergences arising from the phasespace integration of unresolved real radiation [1][2][3][4]. The concrete implementation of this cancellation in perturbative calculations for massless gauge theories is relatively straightforward for low-multiplicity final states and for highly inclusive cross sections, where the involved phase-space integrals and the structure of typical observables are sufficiently simple (witness, for example, the four-loop calculation of the total cross section for annihilation of electroweak gauge bosons into hadrons [5,6]). The situation is considerably more challenging for higher multiplicities and for typical collider observables, where real radiation is subject to intricate phase-space constraints, possibly involving non-trivial recursive jet algorithms. In these cases the phase-space integration must be performed numerically, and the cancellation of soft and collinear divergences is much more difficult to implement. Common approaches involve the definition of approximate real-radiation matrix elements with the correct singularity structure, which are then integrated analytically in order to achieve the required singularity cancellation before numerical tools are employed.
Any solution to the subtraction problem hinges upon our general understanding of infrared divergences in perturbation theory. In particular, the structure of soft and collinear singularities in virtual corrections to scattering amplitudes is very precisely understood [7][8][9][10][11][12][13][14][15][16][17]: divergent contributions to generic massless gauge theory amplitudes can be factorised JHEP12(2018)062 from the hard scattering in terms of a small set of universal functions, defined by gaugeinvariant operator matrix elements. Furthermore, these functions obey evolution equations that can be solved in terms of soft and collinear anomalous dimensions, which are completely known in the massless case up to three loops [18,19]. Differential information on real-radiation matrix elements is somewhat less detailed: the latter have been shown, in considerable generality, to factorise in soft and collinear limits into products of lowerpoint amplitudes multiplied times universal kernels [20][21][22]; all the relevant kernels needed for NNLO calculations are known [22][23][24][25], with partial information available at N 3 LO as well [26][27][28][29][30].
At NLO, such factorisation properties were first employed for the general cancellation of infrared singularities in the so-called 'slicing' approaches [31,32]: these involve isolating singular regions of phase space by means of a small resolution scale (the 'slicing parameter'), approximating real radiation matrix elements by the relevant infrared kernels below that scale, and integrating the latter in d dimensions, so as to explicitly cancel the infrared poles of virtual origin. This procedure yields a correct result up to powers of the slicing parameter, which then has to be taken as small as possible, compatibly with numerical stability. In order to avoid this parameter dependence, 'subtraction' algorithms [33][34][35], were later developed at NLO: in these schemes, one introduces local infrared counterterms containing the leading singular behaviour of the radiative amplitudes in all relevant regions of phase space. One then subtracts the local counterterms from the radiative amplitude, leaving behind an integrable remainder, and one adds back to the virtual correction the exact integral of the local counterterms over the radiation phase space, cancelling explicitly the virtual infrared singularities; the resulting finite cross section can safely be integrated numerically, and the whole procedure is exact, not involving any approximation. These NLO subtraction algorithms are currently implemented in efficient generators [36][37][38][39][40][41][42][43][44], and the handling of infrared singularities is not a bottleneck for phenomenological predictions at this accuracy.
At NNLO and beyond, the construction of general subtraction algorithms is the subject of intense current research. The technical difficulties are significant, due to the proliferation of overlapping singular regions when the number of unresolved particles is allowed to grow, and due to the increasing complexity of the soft and collinear splitting kernels at higher orders. Several schemes have been proposed to address the NNLO problem, belonging either to the slicing [45][46][47][48][49][50][51][52] or to the subtraction [53][54][55][56][57][58][59][60][61][62][63][64][65][66][67] families. Novel ideas are also being introduced [68,69], and the first studies of simple N 3 LO processes have recently appeared [70][71][72]. The variety of NNLO methods developed so far underscores both the phenomenological interest and the technical difficulty of the problem, which so far has not been solved in full generality. It is clear that in the near future it will become phenomenologically relevant, and theoretically interesting, to extend the application of NNLO methods to more complicated processes, and to devise subtraction algorithms at higher orders. Such extensions will require a high degree of optimisation of existing procedures, and possibly the implementation of new methods and theoretical ideas.
In this paper, we propose a theoretical framework to systematically analyse the structure of soft and collinear local subtraction counterterms to any order in perturbation theory. Our guiding principle is the well-understood structure of infrared divergences in virtual cor-JHEP12(2018)062 rections to scattering amplitudes. We note that the detailed structure of virtual factorisation must be reflected in the organisation of local counterterms: this implies significant simplifications, in particular for overlapping soft and collinear singularities, which are straightforwardly handled in the virtual case. Furthermore, we note that explicit high-order calculations of soft anomalous dimensions have shown that many kinematic and colour structures which could potentially contribute to infrared divergences are in fact absent or highly constrained, a feature that must also be reflected in the form of the real-radiation counterterms. Finally, we note that virtual corrections to infrared singularities exponentiate non-trivially, providing connections between low-order and high-order contributions. These interesting and well-understood properties have not so far been fully exploited for the analysis of realradiation subtraction counterterms, and we hope that our discussion in this paper will lead to progress in this direction. Indeed, our central result is a set of definitions for local soft and collinear counterterms, written in terms of gauge-invariant matrix elements of fields and Wilson lines, and valid to all orders in perturbation theory, which can be shown to cancel all virtual and mixed real-virtual singularities on the basis of general cancellation theorems [2,3], and of simple completeness relations. These definitions can easily be shown to reproduce known results at NLO and NNLO, and provide the basis for a first-principle calculation of higher-order universal infrared kernels. Applying this technology at NNLO, we find a simple and physically transparent organisation of soft and collinear subtractions, including in particular the treatment of double counting of the soft-collinear regions.
The paper is organised as follows: in section 2, we briefly review the infrared factorisation of multi-parton scattering amplitudes for massless gauge theories; then, in section 3, we present a basic outline of the subtraction problem at NLO and NNLO: a companion paper [73] is devoted to a detailed construction of a full subtraction algorithm for final-state singularities; in sections 4 and 5, we present our definitions for soft and collinear local counterterms, valid to all to all orders in perturbation theory; in section 6, we briefly illustrate the definitions by showing how they reconstruct the well-understood structure of final-state infrared subtraction at NLO; in section 7 we apply our general results to the problem of NNLO subtraction, and we provide precise expressions for all the local counterterms required for hadronic massless final states; finally, we discuss future developments in section 8.

Infrared factorisation for virtual corrections
We begin by describing the simple multiplicative structure of infrared poles that emerges from the factorisation of fixed-angle multi-particle gauge-theory amplitudes, in order to illustrate the potential simplification that might follow for real soft and collinear radiation. Infrared singularities in these amplitudes factorise in a way which is reminiscent of the renormalisation of ultraviolet divergences: for an amplitude involving n massless particles with momenta p i , the result takes the form [14][15][16] (2.1) In this compact notation, the amplitude A n and the finite coefficient function F n are vectors in the finite-dimensional space of colour configurations, and the divergent factor Z n is a JHEP12(2018)062 colour operator. Soft-collinear factorisation implies evolution equations, which lead to the exponentiation of infrared poles in terms of a finite infrared anomalous dimension matrix Γ n . One may write where all infrared singularities are generated by the integration of the d-dimensional running coupling over the scale λ, extending to λ = 0 [74]. The infrared anomalous dimension matrix Γ n is strongly constrained in the massless case by the factorisation of soft and collinear poles (see eq. (2.6) below). In full generality, one writes where s ij = 2p i · p j , Γ dip n contains only two-particle correlations, and ∆ n is constructed out of quadrupole correlations, starting at three loops [18,19], and constrained to depend on momenta only through the conformal-invariant cross ratios Up to two loops, only the dipole part of the infrared anomalous dimension matrix is relevant. It can be written as [13][14][15][16] where γ i is a collinear anomalous dimension, dependent on particle spin and related to the corresponding field anomalous dimension. The operators T i act as 'gluon insertion' operators, in a manner dependent on the colour representation of the hard particle i, as discussed in [34,75]. The coefficient of the logarithmic term is extracted from the light-like cusp anomalous dimension for colour representation r, γ r K (α s ), assuming that γ r K (α s ) = C r γ K (α s ), and dropping the quadratic Casimir eigenvalue C r : this assumption ('Casimir scaling') is known to be valid up to three loops, while there is solid numerical evidence that it breaks down at four loops, due to the presence of fourth-order Casimir invariants [76,77].
Eqs. (2.3) and (2.5) highlight several remarkable simplifications in the general structure of infrared poles: first of all, exponentiation ties together different orders in perturbation theory in a non-trivial way; furthermore, one observes that correlations involving three coloured particles are absent at NNLO at the level of the soft anomalous dimension, and can only arise in amplitudes through the mixing of one-and two-loop effects upon expanding the exponential; finally, to all orders in perturbation theory, non-dipole corrections are severely constrained to depend on momenta only through the variables in eq. (2.4). We expect these simplifying features to be reflected in the detailed structure of real radiation, and our goal is JHEP12(2018)062 to set up tools to uncover and implement these simplification. In order to proceed, we note that the compact expression in eq. (2.2) is not sufficiently detailed to extract information relevant to the subtraction problem, where it is important to distinguish the contributions of soft and collinear configurations, and to understand the issue of double counting of softcollinear poles. It is therefore necessary to take a step back to the full factorisation formula underlying eq. (2.2), which can be written as [7][8][9][10][11][12][13][14][15][16][17] where for simplicity we suppressed the dependence on the renormalised coupling α s (µ 2 ) and on the regulator . In eq. (2.6), the colour vector H n is a finite remainder (related, but not equal to F n in eq. (2.1)). For each hard massless particle with momentum p i , we introduced a four-velocity vector β i , β 2 i = 0, obtained by rescaling p i by an arbitrary hard scale, say β i = p i /µ, and a 'factorisation vector' n i , n 2 i = 0, responsible for isolating the collinear region for particle i, and in order to enforce the gauge invariance of the collinear factors. For each hard particle, the jet function J i collects all collinear singularities associated with the direction defined by p i . The jet functions are spin dependent, and defined in terms of gauge-invariant matrix elements of fields and Wilson lines. For outgoing quarks with momentum p and spin polarisation s one defines where the Wilson line operator is For (outgoing) gluons with momentum k and polarisation λ, the definition is more delicate, due to the requirement of gauge invariance: a straightforward substitution of a gluon field for the quark field in eq. (2.7) is not satisfactory, due to the non-homogeneous term in the gluon gauge transformation. The issue has been well understood for a long time, initially in the context of giving operator definitions of parton distribution functions for gluons [78].
In that case, the requirement is to find a gauge invariant quantity reducing to a gluon number operator in a physical gauge; a possible solution is to use a particular projection of a field strength operator in place of the gluon field in the equivalent of eq. (2.7): the homogeneous gauge transformation of the field strength can then be compensated by the Wilson line insertion. At amplitude level, an elegant proposal was put forward in the context of SCET in [79,80], and we will use it in what follows. We define where we have not displayed colour indices, the covariant derivative D µ = ∂ µ − ig s A µ is evaluated at x = 0, and the extra power of g s on the left-hand side compensates for the effect of differentiating the Wilson line.

JHEP12(2018)062
We note that jet functions are single-particle quantities and do not carry any colour correlations from the full amplitude: the fact that collinear poles have this property is a highly non-trivial consequence of gauge invariance and diagrammatic power counting. Colour-correlated singularities arise only from soft gluons, which, at leading power in their momentum, cannot transfer energy between hard particles, but induce long-range colour mixing. The soft factor S n is therefore a colour operator, defined in terms of semi-infinite light-like Wilson lines radiating out of the hard collision, each along the classical trajectory of one of the hard particles. One defines where β i is the dimensionless four-velocity of the i-th hard particle, and where, for simplicity, we do not display the color indices of the Wilson lines.
Gluons that are both soft and collinear to one of the hard coloured particles are present both in the jet functions, eq. (2.7) and eq. (2.9), and in the soft matrix, eq. (2.10), and are therefore counted twice. It is however straightforward to subtract this double counting, since the soft approximation of the jet function is simply given by the eikonal jet [12] and soft poles cancel in the ratio of the full jet to the eikonal jet, separately for each hard particle. This simple pattern of cancellation for soft-collinear regions (which in particular does not contain any colour correlations) will be reflected in the structure of local counterterms for real radiation.
We conclude this section with two technical remarks. First, we note that the requirement that n 2 i = 0 for all jet and eikonal jet functions is designed in order to avoid the presence of spurious collinear divergences associated with emissions from the n i Wilson lines. In practical calculations, however, it is highly economical to take the n 2 i → 0 limit, provided one can precisely control the contributions of spurious poles. 1 Finally, we note that in dimensional regularisation all correlators of (semi-)infinite Wilson lines are computed in perturbation theory in terms of scaleless integrals, which vanish in dimensional regularisation, so that the bare soft matrix and eikonal jets equal unity. One can therefore extract the infrared poles of these matrix elements by computing their ultraviolet divergences, which allows to make use of standard renormalisation group arguments. In practice, calculations can be performed with auxiliary regulators for soft and collinear poles: one may for example tilt the β i Wilson lines off the light cone, and introduce a suppression for gluon emission at large distances, as done for example in [82,83]. General theorems then guarantee [84][85][86] that the resulting anomalous dimensions are independent of the chosen collinear and soft regulators.

JHEP12(2018)062 3 Subtraction procedures at NLO and NNLO
We now provide a brief description of a subtraction procedure at NLO and NNLO, for the case of massless coloured particles in the final state, identifying the local counterterms required in this case. Our goal here is to present the general structure of the procedure, which is sufficient for the purposes of the present paper: a detailed construction of a complete subtraction algorithm for this case is presented in [73].
Let us begin by establishing some notation. Given a scattering amplitude with n massless particles in the final state, we write is the Born amplitude for the process at hand (which may of course already contain powers of the strong coupling), while A (k) n (p i ) is the k-loop correction. Given an infrared-safe observable X, one can then construct the perturbative expansion for the differential distribution of X, as At each non-trivial order in perturbation theory, the differential distribution contains contributions with different numbers of final state particles, and the cancellation of infrared singularities takes place upon integration over the phase space of unresolved radiation. Denoting with dΦ m the Lorentz-invariant phase space measure for m massless final state particles, and assuming that the observable involves n particles at Born level, one can write in more detail where δ m (X) ≡ δ(X − X m ) fixes X m , the expression for the observable appropriate for an m-particle configuration, to the prescribed value X. The integrands of the various terms can be expressed in terms of the squared scattering amplitudes involving n, n + 1 and n + 2 particles as

JHEP12(2018)062
where unobserved quantum numbers (such as colour) not affecting the observable X have been implicitly summed over. As briefly discussed in the Introduction, the problem of subtraction arises because the expressions X m for typical observables in the m-particle phase space, as well as the corresponding matrix elements, are very intricate, requiring numerical integrations of the real emission contributions. It is then often necessary to perform the cancellation of infrared poles analytically, before turning to numerical tools. The subtraction approach proceeds by seeking approximations to the real-radiation matrix elements which must be accurate at leading power in the appropriate variables (for instance, energies or transverse momenta) in all singular regions. To be more precise, let us first consider the NLO distribution. In that case, we seek a local counterterm function K n+1 in the (n + 1)-particle phase space, with the requirement that it reproduces the singular behaviour of the real-radiation transition probability R n+1 in all infrared limits, and, in our approach, with the further requirement that it should have a minimal degree of complexity, in order to allow for a direct analytic integration. Given such a function, we define the integrated NLO counterterm as where we introduced the single-particle phase space measure dΦ rad = dΦ n+1 /dΦ n . We can now subtract the local counterterm K n+1 from the real-emission probability R n+1 , obtaining an integrable function in the (n + 1)-particle phase space, and then add back to the distribution the integrated counterterm I n , which must cancel the explicit poles of the NLO virtual correction V n . The result is Note that no approximation has been introduced in passing from the second line of eq. (3.3) to eq. (3.6). Thanks to the infrared safety of the observable X, the integrand in the second line of eq. (3.6) is now integrable everywhere in the (n + 1)-particle phase space, and, at the same time, the first line is free of infrared poles. The differential distribution in this form is therefore amenable to a direct numerical evaluation. At NNLO, the cancellation pattern is considerably more intricate, but an exact subtraction procedure can still be constructed. At this order, infrared singularities arise in three different configurations: in the double-radiation transition probability RR n+2 , either one or two emitted particles can become unresolved, and in the real-virtual transition probability RV n+1 the single emitted particle can similarly become unresolved. It is therefore necessary to define three local counterterms: a function K n+2 in the (n + 2)-particle phase space, approximating RR n+2 in all singular regions with two unresolved particles, a function K (1) n+2 in the (n + 2)-particle phase space, approximating RR n+2 in all singular regions with one unresolved particle, and a function K (RV) n+1 , in the (n + 1)-particle phase space, approximating RV n+1 in all singular regions where the radiated particle becomes JHEP12(2018)062 unresolved. It is furthermore appropriate to separate the double-unresolved counterterm as where the first term collects all double-unresolved limits which are reached hierarchically, with the first particle becoming unresolved at a faster rate than the second one, while the second term contains all remaining double-unresolved contributions, where the two particles become unresolved at the same rate (for a detailed discussion of how to achieve this separation, see [73]). One may then define the respective radiation phase spaces as and introduce the integrated counterterms as n+2 , As was the case at NLO, in eq. (3.5), also in eq. (3.9) the subscripts indicate the number of particles whose phase space still needs to be integrated. Specifically, I n+1 depends on the phase space variables of (n + 1) particles, and has explicit infrared poles cancelling those of the real-virtual transition probability RV n+1 ; the resulting finite combination, however, can still have singular limits when the radiated particle becomes unresolved: those singular limits must be subtracted by combining K (RV) n+1 with I (12) n+1 , in order to cancel the respective explicit poles. Our final expression for the subtracted NNLO distribution is therefore One verifies that no approximation has been made in going from the third line of eq. (3.3) to eq. (3.10). Furthermore, each line in eq. (3.10) is both finite in four dimensions, and integrable in the respective phase spaces. Clearly, eq. (3.10) is only the starting point in the construction of a full-fledged subtraction algorithm: the next crucial step is the explicit definition of the necessary local counterterms, which must properly organise all soft, collinear and soft-collinear regions avoiding double counting; in the process, it is necessary to construct precise phase-space mappings in order to exactly factorise radiative from non-radiative phase spaces; finally, the local counterterms must be analytically integrated in the respective radiation phase spaces. In the remainder of this paper, we discuss a systematic construction of the local JHEP12(2018)062 counterterms, which we will carry out explicitly up to NNLO, but which is applicable in principle at any perturbative order. A detailed algorithmic implementation of eq. (3.10) for final-state massless partons has been presented in [73]. In what follows, our main concern is not the calculation of NNLO kernels, which have been known for a long time [22][23][24][25]: rather, we plan to show how information from the factorisation of virtual corrections allows to organise and simplify the NNLO subtraction procedure, pointing to possible future extensions to higher perturbative orders.

Local counterterms for soft real radiation
Our general strategy to define local counterterms is to construct eikonal form factors and radiative jet functions including real radiation: these functions, when integrated over the final-state phase space and combined with their virtual counterparts using completeness relations, build up eikonal and collinear total cross sections, which are finite by the general theorems of refs. [1][2][3][4]. Let us begin with the case of purely soft final state radiation (which of course includes soft-collinear particles as well). Considering n hard particles, represented by Wilson lines in the soft approximation, radiating m soft gluons, we define the eikonal form factor where in the second line we have defined multiple soft gluon currents J µ 1 ...µm S , in the third line we have introduced the perturbative expansion of the form factors, and we are not displaying colour indices to simplify the notation. A well known property of the soft approximation at leading power in the soft momenta is spin-independence: thus the multiple soft gluon currents are independent of the gluon polarisations λ i , and the definition easily generalises to the emission of final state soft fermions. Note that at this stage the form factor contains loop corrections to all orders in perturbation theory.
Our underlying assumption is that the exact amplitude for the emission of m soft gluons (which may in turn radiate soft quark-antiquark pairs) from n hard coloured particles obeys, to all orders, the factorisation where the remainder R n, m is finite in four dimensions, and integrable in the soft particle phase space. After renormalisation, the amplitude A n, m is ultraviolet finite, and all virtual soft poles, as well as all contributions that are non-integrable in the soft particle phase space, are contained in the soft form factor S n, m . Eq. (4.2) is proven to all orders for m = 0, and it is consistent with all known perturbative results, in particular with the JHEP12(2018)062 arguments of [22,23,25]; a formal all-order proof has however not yet been provided: we treat it as a working assumption, which is known to be correct at NNLO. Squaring eq. (4.2), and performing the trivial helicity sum, one finds, at leading-power in the soft momenta where we introduced the eikonal transition probability Eq. (4.5), up to simple modifications, 2 can be interpreted as an eikonal total cross section. When all coloured particles are in the final state, such a cross section is finite to all orders by the standard cancellation theorems (which can be verified by explicit power counting); with initial state colour, the eikonal cross section is affected by collinear divergences which can be treated by conventional collinear factorisation [87]: indeed, in our framework, these collinear divergences are included in eikonal jet factors to be discussed in section 5. As far as soft divergences are concerned, we conclude that the kernels S n, m provide completely local soft approximations to the relevant squared matrix element, valid at leading power in the soft momenta, and they cancel the virtual soft poles order by order in perturbation theory: this identifies them as candidate counterterms for subtraction in the soft sector.
Let us now illustrate this general framework with simple examples, recovering known results at low orders. A classic case in point is single-gluon emission from a multi-particle configuration at tree level. Eq. (4.2) for m = 1 and at lowest order reads with the definition * (λ) (k) · J (0) (4.7)

JHEP12(2018)062
Explicit calculation expanding the Wilson-line operators in powers of the coupling, or directly with eikonal Feynman rules, easily yields the well-known result for the tree-level soft-gluon emission current [22,75] Squaring the tree-level amplitude one finds the leading-power transition probability where we used the fact that at tree level there is no need to distinguish between H (0) n and A (0) n ; we recognise the colour-correlated Born probability, multiplied times the standard eikonal prefactor. Multiple soft-particle radiation at tree level is similarly easy to compute: for the case of two gluons, one directly recovers the result of [22] with the last line representing uncorrelated emission from two different hard particles, and the first two lines collecting terms arising from double emission from a single hard particle. Currents corresponding to the radiation of soft quark-antiquark pairs, or for emissions with higher multiplicity, can similarly be computed directly in Feynman gauge in a straightforward manner. At loop level, the organisation of counterterms becomes more interesting. Let us for example consider single-gluon emission at one loop: expanding eq. (4.2) for m = 1 to first non-trivial order we find The first term corresponds to a tree-level soft-gluon emission multiplying the finite part of the one-loop correction to the Born process; in the second term the soft function is evaluated at one-loop, and therefore has both explicit soft poles and singular factors from single soft real radiation: it multiplies the Born amplitude. In this case, the proposed factorisation appears to differ from the one proposed in [25], which reads

JHEP12(2018)062
Here the Catani-Grazzini soft current J CG (k, β i ) multiplies the full n-particle amplitude, including loop corrections containing infrared poles, whereas in eq. (4.2) for m = 1 the hard function H n (p i ) is finite. It is, however, easy to map the two calculations, using eq. (4.2) for m = 0, and solving for the one-loop hard part H (1) n (p i ). One finds where we normalised S (0) n to the identity operator in colour space. This leads to an expression for the Catani-Grazzini one-loop soft-gluon current in terms of eikonal form factors, as * (λ) (k) · J (1) (4.14) Comparing eq. (4.14) with the calculation in [25], one easily recognises that the same combination of Feynman diagrams is involved, and one recovers the known result Phrasing the calculation in terms of eikonal form factors allows for a straightforward and systematic generalisation to higher orders. For example, expanding eq. (4.2), for m = 1, to two loops, one finds The expression for H (1) n is given in eq. (4.13); furthermore, one can similarly derive an expression for H (2) n from the two-loop expansion of eq. (4.2) for m = 0, obtaining Substituting the expressions for the hard parts into eq. (4.16), and comparing with eq. (4.12), one finds the two-loop soft-gluon current * (λ) (k) · J (2) (4.18) Note that in expressions such as eq. (4.18) the ordering of factors is important, since the form factors S are colour operators. Note also that all terms in eq. (4.18), except the first one, are already known for general massless n-point Born processes. The two-loop soft-gluon current was computed for n = 2 by extracting it from known two-loop matrix elements in refs. [27,28,88]. Eq. (4.18) provides a precise framework for the calculation for generic processes with n coloured particles at Born level. Clearly, it is not difficult to derive expression similar to eq. (4.18) for the case of multiple soft-gluon radiation at the desired loop level.

JHEP12(2018)062 5 Local counterterms for collinear real radiation
The strategy to define local collinear counterterms is very similar to the one adopted in the soft case. We begin by allowing for further final-state radiation in the operator matrix elements defining the jet functions in eq. (2.7) and eq. (2.9). This leads to the definition of radiative jet functions, which are universal, but distinguish whether the emitting hard parton is a quark or a gluon. In particular, let us consider first a final state with a hard quark carrying momentum p and spin s, and radiating m gluons. In this case we define where we extracted the quark wave function, so that J q, 0 coincides with the virtual quark jet defined in eq. (2.7), and is normalised to unity at tree level. Gluon polarisation vectors, on the other hand, are still included in the function J q, m , and could be extracted to define collinear currents in a manner analogous to what was done in eq. q, m proportional to g 2p+m s . Notice however that the gluon momenta in eq. (5.1) are unconstrained, and collinear limits must be explicitly taken at a later stage in the calculation.
At cross-section level, the definition of radiative jet functions is slightly more elaborate than was the case for soft functions, since one must allow for non-trivial momentum flow. This can be done in a standard way by shifting the position of the quark field in the complex conjugate amplitude, and then taking a Fourier transform in order to fix the total momentum flowing into the final state, setting l µ = p µ i + m i=1 k µ i . In the unpolarised case, one may sum over polarisations and define the cross-section-level radiative quark jet function as The perturbative coefficients J (p) q, m of the radiative jet function J q, m , computed in the collinear limit, provide natural candidates for collinear counterterms, to any order in perturbation theory, as will be illustrated below, in section 6 at NLO and in section 7 at NNLO. For gluon-induced processes, we can proceed in the same way, starting with eq. (2.9), and introducing the (amplitude-level) radiative gluon jet functions as where again we are not displaying colour indices, and polarisation vectors for the radiated gluons are included in the definition of J µν g, m . The definition (5.3) can be used to construct a cross-section-level radiative gluon jet function, as was done for the quark. It reads To illustrate the usefulness of radiative jet functions as collinear counterterms, let us focus, as an example, on the quark-induced jet function. In analogy to what was done in the soft sector, we note that summing over the number of radiated particles, and integrating over their phase space, by completeness one finds ∞ m=0 dΦ m+1 J q, m (k 1 , . . . , k m ; l, p, n) The r.h.s. of eq. (5.5) gives the imaginary part of a generalised two-point function, which is a finite quantity, since it is fully inclusive in the final state. The m = 0 contribution contains the virtual collinear poles associated with an outgoing quark of momentum p, and therefore the real radiation contributions for m = 0, given by eq. (5.2), must cancel those poles order by order in perturbation theory, as desired. Inclusive cross-section-level jet functions such as the integrated quantity in eq. (5.5) have been used in the context of threshold resummations for many years, starting with the seminal papers in refs. [89,90]. We can perform a simple test of the correctness of our method by computing the single-gluon radiative jet for an outgoing quark with momentum p µ . In the lowest perturbative order in the coupling constant receives contributions from three different diagrams, shown in figure 1, which give the result where p 2 = k 2 = 0, and up to corrections proportional to n 2 . It is easy to trace the contributions of the three diagrams in figure 1 in the axial gauge calculation of ref. [22]. Notice however that in eq. (5.6) the collinear limit for k, corresponding to l 2 → 0, has not been taken yet. This is easily achieved by introducing a Sudakov parametrisation for momenta p µ and k µ , and taking the k ⊥ → 0 limit, setting Due to the prefactor of order O (l 2 ⊥ −1 ], the leading behaviour in the l ⊥ → 0 limit is recovered by setting l ⊥ = 0 in the square bracket. This yields up to corrections of order l ⊥ . In the square bracket, as expected, we recognise the leading order unpolarised DGLAP splitting function P q→qg . It is interesting to perform the same check for the cross-section-level radiative gluon jet definition, which must reproduce the splitting kernel P µν g→gg when m = 1. The diagrammatic contributions, in Feynman gauge, are similar to those in figure 1, and are displayed in figure 2; in an axial gauge, n·A = 0, only the third graph, figure (2c), survives. In Feynman gauge, at amplitude level, the single-radiative jet function defined in eq. (5.3) gives which can be verified to be consistent with the computation performed in axial gauge.
Computing the single-radiative gluon jet function at cross-section level, we can use the Sudakov parametrisation

JHEP12(2018)062
To leading power in l ⊥ , and setting n 2 = 0, we end up with the expression The first two terms in the square bracket reproduce the expected splitting function; the third term, where the braces denote index symmetrisation, is proportional to either l µ or l ν : in the collinear limit, these corrections vanish when contracted with the factorised hard amplitude, which depends on the on-shell parent gluon momentum l. It is easy to check, by considering a final-state qq pair in eq. (5.3), that one may similarly recover the appropriate splitting function P µν g→qq ; kernels for double collinear emission can be reproduced with similar manipulations.
To complete our discussion, we note that the cross-section-level jet functions presented in eq. (5.2) generate all collinear singularities, including soft-collinear ones. These are therefore double counted, since they were already included in the soft region. In order to avoid this issue, following the logic suggested by the factorisation of virtual corrections in eq. (2.6), we may introduce radiative eikonal jet functions, defined by replacing the field ψ(0) in eq. (5.1) with a Wilson line (in the same colour representation), oriented along the hard parton direction β ν = p ν /µ. At cross-section level, this leads to the definition Notice that the radiative eikonal jet does not depend on the spin of the hard parton, so that eq. (5.12) applies to gluons as well; the Fourier transform fixes l µ to be the total momentum of the final state.
To test this definition, we compute the soft-collinear local counterterm for single radiation, and we easily find 2p · n p · k n · k . (5.13) In the limit of p µ collinear to k µ , we can employ the relations to obtain the explicit soft-collinear counterterm We note that the factor 2z in the numerator is necessary to enforce the commutation relation between soft and collinear limit at NLO: a basic feature that allows significant simplifications in the subtraction procedure [73].

JHEP12(2018)062
6 Constructing counterterms at NLO Our basic strategy for subtraction is to identify soft and collinear local counterterms starting from the known expressions for the poles of virtual corrections. We now proceed to illustrate how this works with the simple case of NLO massless final-states. Expanding eq. (2.6) to NLO, and using the fact that virtual jet functions are normalised to equal unity at tree level, we easily find Using eq. (6.1), it is straigthforward to construct the NLO virtual correction V n , entering NLO distributions as in eq. (3.3), and to express it in terms of the cross-section-level soft and jet virtual functions. One finds It is now a simple task to find local counterterms for these poles: one simply notices that the soft completeness relation in eq. (4.5), at NLO, implies the cancellation The local phase space integrands in eq. (6.3) and eq. (6.4), multiplied times the appropriate, finite, hard coefficients, must thus provide the necessary counterterms. In particular NLO soft poles are cancelled by integrating the combination note that, for gluons, the function J i, 1 is a spin matrix acting on the spin-correlated Born. The double subtraction of soft and collinear singularities overcounts the soft-collinear regions: one must therefore add back a local soft-collinear counterterm, given by i, E, 1 (k i ; l, β i , n i ) A (0) n (p 1 , . . . , p i−1 , l, p i+1 , . . . , p n ) 2 . (6.7) Using the tree-level results listed in section 4 and in section 5, it is easy to see that eq. (6.5) and eq. (6.6) reproduce standard results for NLO subtraction. One should however appreciate that the present approach provides a simple proof that the list of singular regions for real radiation considered here is exhaustive, and collinear regions for radiation from different outgoing hard particles do not interfere. While these facts are well-understood at NLO, their generalisations at higher orders are much less obvious. On the other hand, we note that these result do not yet constitute a subtraction algorithm at NLO: indeed, one can see that the tree-level matrix elements appearing in eq. (6.6) involve particles that are not on the mass-shell, except in the strict collinear limit, while momentum conservation is not properly implemented in eq. (6.5), except in the strict soft limit. A practical algorithm must provide a resolution of these issues, with the construction of suitable momentum mappings between the Born and the radiative configurations, either with global treatment of phase space, as done for example in [34], or with a decomposition into different singular regions, as done for example in [33] and in [73].

Constructing counterterms at NNLO
Extending the procedure of section 6 to higher orders is in principle straightforward, but it unveils and organises several non-trivial features of real radiation in singular regions of phase space. Let us begin by extending eq. (6.1) by computing the expansion of the virtual correction to the amplitude up to NNLO. The two-loop contributions can be written as Several comments are in order. We begin by noting that the first term on the first line is finite, being given by the action of the finite tree-level soft operator on the two-loop finite hard remainder. The second term contains two-loop soft and soft-collinear poles from the soft operator, giving singularities up to the maximum allowed degree, 1/ 4 . In the third term the one-loop soft operator acts on the one-loop finite hard remainder, giving a single soft pole and a double soft-collinear pole. The second line is the most interesting from the point of view of factorisation: it contains all double hard-collinear poles arising from twoloop virtual corrections associated with a single hard external leg, yielding singularities up to 1/ 2 . In particular, the second line does not generate any soft poles: indeed, while the function J (2) i (p i ) contains up to two soft poles, generated by gluons that are both soft and collinear to the i-th hard particle, the contributions in which both gluons are soft (on top of being collinear) are cancelled by the second term in square bracket, J  last term in the square bracket. Notice the factorised form of that last term: when one gluon is hard and the other one is soft, the soft gluon factorises from the matrix element in the usual way. This cancellation mechanism is illustrated, for a sample diagram, in figure 3. The last two lines in eq. (7.1) have a simpler interpretation: the third line contains single hard collinear poles arising simultaneously on two different hard legs, i and j; the fourth line contains single hard collinear poles on the i-th hard leg, accompanied by a soft single pole, or a soft-collinear double pole, or just multiplied times a finite correction.
The next step is to construct the virtual contributions to the squared amplitude at NNLO. In order for our procedure to work, these must in turn be expressed in terms of the cross-section-level virtual jet and soft functions, which is less than trivial since, at NNLO, all functions involved receive contributions both from the interference between the Born amplitude and the two-loop correction, and from the square of the one-loop amplitudes. For example, the two-loop cross-section-level virtual soft function is given by while the two-loop unpolarised cross-section-level radiative jet function for a quark emitting m gluons reads It is relatively simple to organise the virtual poles in the real-virtual contribution to the squared matrix element: this amounts essentially to a repetition of the NLO calculation, with n + 1 hard particles in the final state. One easily finds Double virtual poles, on the other hand, receive several non-trivial contributions, which can be organised as follows: n, i . (7.5)

JHEP12(2018)062
We will now go through the various contributions to the r.h.s. of eq. (7.5), identifying in each case the real radiation counterterms that are needed to cancel the corresponding virtual poles. The double-soft virtual contribution (V V ) (2s) n , which, within our chosen organisation of the matrix element, contains soft-collinear poles as well, is given by where S n, 0 was given in eq. (7.2). To give a complete picture of the soft sector, at this point we also include in the discussion the single-soft virtual contribution (V V ) (1s) n , which is given by as well as the real-virtual soft poles in eq. (7.4). In order to cancel these poles, we need the completeness relation for the soft sector to NNLO, which reads It is natural at this point to identify three separate soft counterterms, characterised by their kinematic structure. The double-unresolved soft counterterm involves n-point hard kinematics, and double soft radiation; it is given by The single-unresolved soft conterterm involves (n + 1)-point hard kinematics, and single soft radiation; it is given by We now note that this procedure yields an expression for the complete double-unresolved soft counterterm K (2s) n+2 , but does not immediately distinguish between the two contributions defined in eq. (3.7). It is however easy, in this context, to identify the desired partition of the counterterm. Indeed, as discussed in ref. [73], the local counterterm K (12) n+2 is designed to be integrated in two stages: the first integration, in a single-particle phase space, yields the integrated counterterm I (12) n+1 , which must cancel the explicit poles of the real-virtual counterterm K (RV) n+1 , given entirely by the last term in eq. (7.11). From the point of view of factorisation, the desired function is then identified as follows. An explicit calculation of S (0) n, 2 from its definition in eq. (4.4) yields the sum of two distinct contributions: one in which the soft limits on the two radiated gluons are taken hierarchically, with one gluon being much softer than the other one, and one in which the two gluons have a comparable softness. This structure was identified in ref. [22], and is derived from eq. (4.10) by taking JHEP12(2018)062 the limit in which k 2 is much softer than k 1 , or viceversa. The hierarchical limit of K (2s) n+2 is constructed essentially by treating one of the two soft radiated particles temporarily as a hard one: it gives therefore precisely the desired function K (12, 2s) n+2 , which, upon integration, will cancel the explicit double-soft poles of the real-virtual local counterterm. A similar pattern can be replicated for the other double-unresolved local counterterms, in all cases in which a hierarchy between the two unresolved particles can be identified.
Turning to hard collinear poles, we first tackle the contribution with two hard collinear virtual gluons attached to the same hard outgoing leg. It is given by (7.12) In order to cancel the poles of the first two terms in eq. (7.12), we can use the NNLO expansion of eq. (5.5), which gives the finiteness condition i, 2 = finite , (7.13) and the analogous expression for eikonal jets. The third term of eq. (7.12) has a different structure, since it is a product of two one-loop functions. One can however cancel its poles with the same general approach, by using the fact that Once again, the contributions to different local counterterm functions can be identified by their phase space structure. We define while no single-unresolved counterterm in the (n + 1)-particle phase space is required in this case.

JHEP12(2018)062
are not paying special attention to the issue of Glauber gluons and possible factorisation violations: essentially, we are assuming that the formalism applies for sufficiently inclusive observables, such that Glauber gluons do not result in uncancelled infrared singularities. The issue is however very interesting from a theoretical point of view: infrared power counting shows that Glauber gluons do not contribute to infrared singularities for fixedangle scattering amplitudes (see, for example, [91] for a recent discussion of leading regions in coordinate space), but they can result in a breakdown of factorisation for insufficiently inclusive hadronic cross sections, when real collinear radiation is integrated over unresolved regions of phase space (see, for example, [92][93][94] for recent discussions). The tools developed in this paper, which allow for the study of real infrared radiation at the level of differential distributions, may in future help to shed light on the limits of factorisation theorems for less inclusive collider observables. 3 At the level of definitions of local counterterms, the extension to massive partons (which is of obvious phenomenological interest, in view of top-quark-related observables, and possibly b-quark mass effects) is not problematic: indeed, massive partons are not affected by collinear divergences (although it may be of interest to resum collinear logarithms of the quark mass), so that the structure of counterterms is in fact simpler when masses are present. In the massive case, on the other hand, more work is needed to properly define the phase space mappings associated with branchings involving massive partons [96], and to perform the corresponding integrations. On the other hand, the approach we have presented here is likely to have a significant impact in the organisation of future N 3 LO subtraction algorithms: indeed, at N 3 LO, the combinatorics of overlapping singular regions becomes considerably worse, and the impact of infrared exponentiation on subtraction is bound to become stronger. Work on a detailed extension of the present work to N 3 LO is in progress.
More generally, we hope that the present work will contribute to developing our knowledge of the infrared behaviour of real radiation at the differential level, to all orders in perturbation theory, bringing it to the same detailed level of understanding and control currently enjoyed by virtual corrections to fixed-angle scattering amplitudes and by inclusive cross sections.