Momentum mappings for subtractions at higher orders in QCD

Subtraction schemes provide a systematic way to compute fully-differential cross sections beyond the leading order in the strong coupling constant. These methods make singular real-emission corrections integrable in phase space by the addition of suitable counterterms. Such counterterms may be defined using momentum mappings, which are parametrisations of the phase space that factorise the variables that describe the particles becoming unresolved in some infrared or collinear limit from the variables that describe an on-shell phase space for the resolved particles. In this work, we review existing momentum mappings in a unified framework and introduce new ones for final-collinear and soft counterterms. The new mappings work in the presence of massive particles and with an arbitrary number of soft particles or of clusters of collinear particles, making them fit for subtraction methods at any order in perturbation theory. The new mapping for final-collinear counterterms is also used to elucidate relations among existing final-collinear mappings.


Introduction
The discovery of the Higgs boson [1,2] at the CERN Large Hadron Collider (LHC) and the subsequent study of its properties has been the crowning achievement of the Standard Model (SM). These results fixed the value of the last unknown parameter of the SM and confirmed many predictions of the properties of the Higgs boson, such as the value of its coupling strengths to the SM fermions. While the LHC, as a hadron collider, was designed as a discovery machine, it is now also our best tool to improve the precision of our knowledge of the Standard Model, both in the Higgs sector, where we have much to learn, and in the other aspects of the theory. As statistics are steadily increasing and experimental techniques improving, the experimental uncertainties are progressively shrinking and projections for the end of the LHC programme indicate that it will provide precision results in the Higgs sector as well as on the strong and electroweak interactions. In many cases, in particular for Higgs physics, it is expected that the leading source of uncertainty on SM observables will be due to theoretical predictions. As a result, it is crucial that we improve the precision of theoretical calculations, which can often be achieved by going to higher orders in perturbation theory. The next-to-leading order (NLO) both in the Quantum Chromodynamics (QCD) and in the electroweak (EW) couplings is now routinely achievable for most processes thanks to the development of highly efficient automated tools. Since in many cases of interest NLO predictions are not sufficient to match the experimental precision, the next target are next-to-next-to-leading order (NNLO) calculations in the QCD coupling for which no fully automated method exists yet.
At the NLO, exploration of this IRC behavior of amplitudes has led to the establishment of subtraction methods [29,30] as the standard approach, in which real-emission corrections, which exhibit IRC divergences when integrated over phase space, are made integrable by the addition of suitable counterterms. These counterterms exhibit the same divergent behaviour as real-emission matrix elements, but are typically much simpler functions of the kinematics. This permits to integrate them in dimensional regularisation, expose their singularities as poles in the regulator and cancel them against the IRC poles of the virtual loop corrections -as predicted by the KLN theorem.
All the subtraction schemes mentioned above are based on devising a set of counterterms which approximate the matrix element in the limits where they become singular, such that the difference be computable in four dimensions and the counterterms can be integrated in d-dimensions as a Laurent series in the regulator. Defining a scheme amounts to specifying how these counterterms are obtained for a generic process, or at least a generic class of processes. There is of course a lot of freedom in the selection of counterterms, as is exhibited by the number of well-established schemes and by the continued appearance of new approaches [63][64][65].
All the existing schemes exploit the factorisation properties of amplitudes in the singular limits, which are the points in phase space where some massless partons are collinear to one another or soft. The leading behaviour of the amplitude close to the singular surfaces in phase space is known, and in order to achieve a working subtraction it must be the case that the sum of all counterterms features the same leading behaviour. Counterterms, however, are functions of phase space as a whole (possibly equal to zero in parts of it in the case of sector-based approaches) and one must define them away from the singular limits. This is where schemes differ from one another, since the definition of their singular limits is far from enough to make the counterterms unique.
One aspect of the freedom of choice in the definition of counterterms for subtraction schemes is momentum mappings. Momentum mappings are parametrisations of the phase space where the variables that describe the particles becoming unresolved in some infrared or collinear limit are factorised from the variables that describe an on-shell phase space for the resolved particles. This factorisation property is key to make a subtraction method general using the following procedure: • for a given limit, we choose a mapping that factorises phase space into a lowermultiplicity phase space and "unresolved variables"; • we write our counterterm as a universal singular factor (e.g. an Altarelli-Parisi kernel) multiplied by a squared amplitude with the appropriate reduced multiplicity; • we choose the momenta of the lower-multiplicity phase space as the arguments of the lower-multiplicity squared amplitude and we express the universal singular factor in terms of the "unresolved variables" and possibly of the lower-multiplicity phase-space momenta. The factorisation of both the integrand and the parametrisation allows then the ddimensional integration of the singular factor over the "unresolved variables" for fixed lowermultiplicity kinematics. Thanks to this procedure, the integral of the real-emission contribution yields a Laurent series in the dimensional regulator, whose poles cancel the poles of the corresponding virtual squared amplitudes locally in the lower-multiplicity phase space.
Momentum mappings were introduced at NLO where the original dipole subtraction [30] uses a factorising parametrisation for massless partons, of which that of the CoL-oRFul subtraction [66] is a variation. A variant of the mapping of Ref. [30] was introduced to handle massive particles [67] and another solution was proposed subsequently by Nagy and Soper [68]. As we shall see in section 2, different schemes use different mappings at NNLO.
Because momentum mappings are always a means to an end, their properties and defining features have been little studied so far. However, the current activity in setting up subtraction schemes that are truly general calls for a transversal study of momentum mappings. The aim of the present paper is to make a step in this direction. More precisely, we will review a number of existing momentum mappings in a unified framework, introduce new ones and present important observations about their application.
The paper is organised as follows. Section 2 provides a more detailed discussion of what a momentum mapping is by setting up explicitly their definition at NLO for collinear and soft configurations and showing how they are realised in different existing schemes both at NLO and NNLO. Section 3 introduces a new momentum mapping for final-collinear coun-terterms which is shown to be a generalisation of existing mappings. The new momentum mapping works in the presence of massive particles and with an arbitrary number of clusters of collinear particles. It can be used to show the equivalence of the mappings of Refs. [67] and [68] in the case of massive particles with a single recoiler. Section 4 introduces a new momentum mapping for soft counterterms which works in the presence of massive particles, as well as in specific kinematic configurations where existing soft mappings fail. Section 5 presents ideas to allow subtraction schemes to be setup in a way where integrated counterterms can be computed once and for all independently of the choice of mappings, which we realise explicitly on the specific case of final state NLO collinear counterterms. Section 6 presents our conclusions.

Review of momentum mappings
This section is dedicated to introducing momentum mappings with explicit examples, setting up definitions and providing a resource where the different approaches used in the literature and their properties are described. We will first review the example of CoL-oRFul subtraction for final-state NLO singularities, which is built on both a soft and a collinear mapping, allowing us to illustrate both important aspects and possible issues relating to mappings. We then move on to discuss the three main different types of mappings and their realisations in different subtraction schemes. Finally we discuss how these elementary mappings can be combined to handle counterterms where disjoint sets of particles become unresolved.

An example at NLO: CoLoRFul subtraction
As a first example, let us describe how CoLoRFul subtraction handles the regularisation of the real-emission corrections to the process e + e − → qq where q is a massless quark, i.e. e + e − → qq g. The matrix element for this process diverges in three limits: when the gluon momentum becomes zero (p g → 0), or when it becomes collinear to the momentum of either the quark (p g p q ) or the anti-quark (p g pq). The well-known collinear and soft factorisation formulae for the squared amplitude multiplied by an IRC-safe observable O read as follows, where zp q = (1 − z)p g and z pq = (1 − z )p g in the respective collinear limits, µ is the regularisation scale, = (d − 4)/2 is the dimensional regulator, α s is the strong coupling and C F is the fundamental Casimir of the strong interaction gauge group. The amplitudes M qq and M qqg are respectively the amplitudes for the Born (e + e − → qq) and real-emission (e + e − → qq g) processes. Note that amplitudes are unambiguous only when their arguments are momentum conserving and on-shell, meaning that both |M qq (p q , pq)| 2 and |M qq (p q + p g , pq)| 2 are well-defined only in the appropriate limits, as momentum conservation between the initial and final state would otherwise not be respected in the former and the on-shell condition for the quark would be violated in the latter. As we discussed in the introduction, a momentum mapping resolves this ambiguity by defining mapped momentã p q andpq as functions of p q , pq and p g such that: (i) they are on-shell and momentum conserving outside of the strict limit; (ii) the integration measure is factorised. Next, we illustrate how this works.

Soft mapping
The issue with using the right-hand side of eq. (2.1) to define a soft counterterm is that the momenta p q and pq do not add up to the initial state total momentum Q = p e + + p e − , such that the matrix element is not well-defined outside of the strict limit. The solution proposed in the CoLoRFul scheme is to use instead the momenta, is the ratio of centre-of-mass energy of the qq system to the one of the qqg system, and Λ[λQ, Q − p g ] is a Lorentz transformation that maps Q − p g to λQ, given in eq. (2.30). The role of Λ is easy to understand in the rest frame of Q, where it ensures thatp q +pq has a null total 3-momentum. Replacing p q/q withp q/q does not affect the leading behaviour of the right-hand side of eq. (2.1), since in the soft-gluon limit p g → 0 we have λ → 1 and thereforep q/q → p q/q . As a result, a valid soft counterterm is Note that the momenta in the eikonal factor are taken to be the mapped momenta. This is an arbitrary choice, since choosing to keep the original quark momenta in the eikonal factor would again yield the same leading behaviour in the soft limit. This counterterm provides an appropriate regulation for the soft divergence of the matrix element, as it is a well defined function of the same variables as the matrix element and reproduces its behaviour in the soft limit. In the computation of the local counterterm required for integrating the realemission matrix element over the 4-dimensional phase space, property (i) is exploited: the arguments of the reduced matrix element |M qq (p q ,pq)| 2 are on-shell momenta that verify momentum conservation, which allows it to be unequivocally defined and efficiently derived. Property (ii) relates to the d-dimensional integration of the counterterm, which is performed to expose the poles in the analytic regulator at the local level in the Born phase space. In the phase-space integral that defines the cross section, the CoLoRFul soft mapping provides a change of variables that allows us to factorise the integration measure as follows, Using this expression, we can now write the d-dimensional integral of the counterterm, The universal integrated counterterm is the integral over the unresolved gluon momentum of the eikonal factor, Note that the dependence on Q indicated on the left-hand side is not explicit in the formula, but is generated by λ defined through eq. (2.5). The integral in the last equation can be done analytically once and for all in d-dimensions and it exposes the implicit soft singularities of the real-emission matrix element as explicit poles in the dimensional regulator. This leaves the integral over the Born phase space, which is usually done numerically as the integrand is now integrable in the limit d → 4 and the observable function can be arbitrarily complicated, Furthermore, the singular eikonal factor for a soft gluon emitted from a quark-antiquark dipole is universal and the phase-space factorisation used above is generalisable to arbitrary final states as we will illustrate below. As a result, the integrated singular factor S qq g will appear unchanged for arbitrarily complicated processes that feature this type of singularity.

Collinear mapping
The discussion of the CoLoRFul approach for the subtraction of the NLO soft divergences already outlined the main advantages of using momentum mappings for subtraction: they allow for the definition of reduced kinematics which completely separate the variables which control a singular limit of the matrix element from the variables that describe the reduced process, i.e. the variables that enter the matrix element and the observable. This in turns permits the analytic integration of the singular factor independently of the process and the observable.
This is of course also true for collinear singularities and their associated mappings. Let us consider the limit g q in the example at hand; the issue is that the momentum of the parent quark p q + p g which appears on the right-hand side of eq. (2.2) is on-shell only in the exact limit. The collinear mapping used in CoLoRFul defines a valid parent momentum via a "democratic" shift proportional to the total momentum Q of the process, and adjusts all other (massless) momenta through rescaling, As in the case of the soft mapping, an essential feature of this mapping is that the transformation between p qg , pq andp qg ,pq becomes trivial in the collinear limit. Indeed in this limit y gq → 0, therefore α → 0 which in turn impliesp qg → p qg andpq → pq. One can then define a local counterterm that is valid over all of phase space and reproduces the collinear limit of the matrix element as The real-emission phase space can be rewritten as a convolution over s qg , and using the mapping one finds where J s qg , Φ qgq is the Jacobian of the change of variable (p qg , pq) → (p qg ,pq) at fixed s qg . As in the soft case, we can then separate the integral over the reduced phase-space variables Φ qgq that enter the matrix element and observable in the counterterm from the integral of the singular factor taken over the decay phase space of the off-shell parent quark with momentum p 2 qg = s qg to the quark-gluon final state, Φ qg . This decay phase space depends on the momentum of the parent particle p qg , which is itself a function of s qg and Φ qgq through the mapping, which we explicit in the integral measure above. This phasespace factorisation yields the integrated counterterm, which, as in the case of the soft counterterm, features the universal integrated singular factor, The fact that this integral only depends onp qg and Q is specific of the mapping adopted for this scheme, and in general the dependence can be over the whole reduced phase space. This integral can be performed once and for all to expose the related IRC phase-space singularities as poles in the dimensional regulator, leaving the reduced phase-space integral, to be done numerically for arbitrary observables after it has been combined with the virtual contribution. Now that we have illustrated the role that mappings play in subtractions, let us review the different choices that have been made in the existing methods for NLO and NNLO subtractions. We will first focus on mappings used to define counterterms that regulate final-collinear limits, then those associated to initial-collinear limits and finally those used to regulate soft limits. 1

Final-collinear mappings
The basis of all the final-collinear mappings we present in this section is the well-known factorisation of arbitrary phase spaces into the production of a parent particle and its decay, already used in eq. (2.14). Let us consider n = m + k − 1 final-state particles with momenta {p} n = {p 1 , . . . , p n }. Their phase space may be factorised as follows, What we call a final-collinear mapping is a change of variables In practice, in order to regulate the divergences associated to k final state particles becoming collinear, we are interested in keeping m i unchanged for i = 1, . . . , m − 1 and settingm K = 0 so that p K describes the momentum of an on-shell massless parton. The existence of a singularity in the (m . . . n)-collinear limit also depends on having m m = · · · = m n = 0.
In this section we will always consider that all momenta not involved in the collinear splitting are mapped. It is however of course always possible to leave some momenta unchanged and to affect only a subset through the mapping. When discussing this possibility, also in the context of initial-collinear and soft mappings, we will refer to the momenta that are affected by the mapping as recoilers while we will refer to those that are not changed as spectators.

Dipole mapping
In the dipole subtraction scheme for NLO subtractions [30], the following mapping is used to define an on-shell reduced matrix element for the (s, t)-collinear limit: we specify a massless recoiler momentum, which we take here to be p r , and definẽ with p a 1 ...an = p a 1 + . . . + p an . It is easy to see that momentum conservation is obtained since we have p r + p s + p t =p r +p st . One can easily generalise this mapping to k unresolved The phase space for n = m + k − 1 final-state particles is written explicitly as

Rescaling mapping
The rescaling mapping, which is used in the CoLoRFul scheme, can be seen as a generalisation of the dipole mapping where the pair of collinear momenta (s, t) recoils against all other final-state particles instead of a single recoiler r. In order for the rescaling mapping to be applicable all recoilers need to be massless. Indeed, each of the resolved momenta is rescaled by an appropriate factor to restore momentum conservation as shown in eq. (2.11) and below,p It is again easy to generalise this to multiple unresolved momenta t 1 , . . . , t k [44] as follows, Incidentally, we note that in the centre-of-mass frame α t 1 .
The relation in eq. (2.24) is invertible and allows us to express α t 1 .
We can then write the following phase space factorisation formula for n = m + k − 1 final-state particles, In this case, a parametrisation that leads to a simpler expression for the phase space is the one already adopted in the CoLoRFul scheme,

Lorentz mapping
An issue with the rescaling mapping defined above is that the resolved momenta need to be massless so that the rescaling does not change their mass. This holds both for the parent of the collinear particles and for the other final-state recoilers. A quick fix to apply the rescaling mapping in the presence of massive final-state particles would be to recoil only against the massless ones. Depending on the process this is however not always possible, as in the case of the real-emission contribution e + e − → tt g g in the limit where the two gluons become collinear. An option that is applicable in this case, proposed by Nagy and Soper [68], is to restore momentum conservation in the reduced phase space using a Lorentz transformation. In the case of two unresolved massless momenta s and t it takes the form, and Λ [Q −p st , Q − p st ] is a Lorentz transform that maps the 4-vector Q − p st to Q −p st , ensuring momentum conservation. 3 It is possible to define such a transformation independently of the space-time dimension as Such a Lorentz transformation only exists if the two momenta have the same non-vanishing invariant mass K 2 =K 2 . The expression is therefore valid unless there is only one massless final-state recoiler. When this situation arises, as e.g. in the case of dijet production, another form of the Lorentz transformation Λ must be adopted, as for instance the one described in appendix A. It is easy to verify that indeed eq. (2.30) gives Λ[K, K]K =K and Λ · Λ T = 1.
The mapping (2.29) is again generalisable to k unresolved massless momenta p t 1 , . . . , p t k , (2.31) The phase space factorisation for n = m + k − 1 final-state particles is Note that, as shown in appendix B, this mapping can be generalised to the case where some of the unresolved momenta, as well as the parent momentump t 1 ...t k , are massive, which makes it suitable for quasi-collinear counterterms in processes with massive coloured particles.

Initial-collinear mappings
Initial-collinear mappings are used to subtract divergences that occur when a set of finalstate particles become collinear to an initial-state particle. There are two main differences compared to final-collinear mappings, • initial-collinear mappings generate a convolution on Bjorken momentum fractions that cannot be cast into a factorised form; • the unresolved phase space contains one particle less than the number of final-state particles that are not resolved. In this section, we provide a single example of mapping, which was used at NLO in the original dipole subtraction [30], in antenna subtraction [35], in CoLoRFul subtraction [48] as well as in Frixione-Kunszt-Signer (FKS) subtraction [29]. Let us consider the factorisation of the amplitude for a process p a , p b → p 1 , . . . , p m+1 when p m+1 p a , for any reference 4-vector n such that n · p a = 0. We aim to parametrise the (m + 1)-particle phase space with total momentum Q = p a + p b in terms of variables which describe the (m + 1)-th emission and an m-particle phase space with total momentumQ = xp a + p b . Contrary to the case of final-collinear mappings, we thus not only have to change the finalstate momenta but also the initial-state ones. We set n = Q and realise the mapping as follows,p It is easy to observe that this mapping does indeed achieve its intended goal: in the p m+1 p a limit,p a = p a − p m+1 , soQ = Q − p m+1 and therefore alsop i = p i . Note that all final-state particles need to be shifted for the mapping to work. The phase space can be reparametrised in terms of the new momenta, where x in the second line is given by eq. (2.34) and is therefore a function of p a , p m+1 and n = Q. Using this mapping we do not completely factorise the phase space, but obtain a convolution where ξ entangles the energy of the emitted unresolved particle with the resolved initial state momentum. As a result, counterterms integrated over the unresolved degrees of freedom [dp m+1 ] will feature a dependence on ξ.
As in the case of final-collinear mappings, this transformation can be generalised without effort to k particles becoming collinear to the initial momentum p a [35], yielding the phase space convolution, [dp m+1 . . . dp m+k ](ξ). (2.38)

Soft mappings for massless partons
The soft mapping of section 2.1.1 is trivially extended to m massless partons of momenta p 1 , . . . , p m and one soft gluon of momentum p s . The m mapped momenta are defined by first rescaling all the hard momenta by a factor 1/λ and then Lorentz-transforming all of the rescaled momenta [66] p µ i = Λ µ ν Q, Here Λ is given in eq. (2.30) withK = Q and K = (Q − p s )/λ. The constraint K 2 =K 2 implies that λ = 1 − y sQ . (2.40) The phase space of eq. (2.7) is generalised to m hard partons, where the m momenta in the first factor on the right-hand side are those of eq. (2.39).
It is straightforward to extend the single-soft mapping of eq. (2.39) to a multiple-soft mapping for m massless hard partons of momenta p 1 , . . . , p m and r soft partons of momenta p s 1 , . . . , p sr , with n = m + r. The m momentap 1 , . . . ,p m are given bỹ where Λ is given by eq. In the case of two soft partons, the mapping of eq. (2.42) was used in [44].

Generalised rescaling mapping
In this section, we introduce a transformation which reparametrises the m-particle phase space of the momenta {p} m with masses {m} m in terms of m momenta {p} m with different masses {m} m and the same total momentum Q. We propose a novel application of this transformation as a final-collinear momentum mapping. The main application is of course the subtraction of genuine IRC singularities, where it is used to replace momenta of sets of massless final-state particles going collinear to each other with on-shell momenta for their massless parents. We foresee that other applications will be relevant, such as the stabilisation of quasi-collinear singularities in processes with massive coloured particles. We introduce the transformation in section 3.1, then we outline its usage as a mapping and derive the corresponding phase-space factorisation in section 3.2. We highlight important properties of this mapping in section 3.3, give an explicit application in section 3.4 and comment on counterterm integration in section 3.5. Finally, in section 3.6 we point out several special cases in which the transformation is significantly simpler to formulate or reduces to one of the mappings defined in the previous section.

Definition
We begin by defining how the mapping acts on the {p} m phase space. Although the mapping can be formulated in a manifestly covariant form, as we shall see in eq. (3.10), for the sake of simplicity we work in the rest frame of the total momentum Q and use nonexplicitly Lorentz-covariant notation. In the considered frame, the 3-momenta involved in the mapping sum to zero, Spatial momentum conservation therefore remains valid if all 3-momenta are rescaled by a common arbitrary factor κ, Energies can then be set by imposing the mass-shell conditions, and the parameter κ finally be fixed by requiring energy conservation, where we abuse notation by using Q to refer to Q 2 . Once the values for the target masses {m} m are given, this is an algebraic equation for the unknown κ. The left-hand side of eq. (3.4) is a monotonous function of κ which varies between im i and ∞ over κ ∈ R + , and therefore admits a unique valid solution as long as the physical condition im i ≤ Q is respected. 4 The solution for κ can be promptly written in closed form if only two or three momenta are involved. For a generic number m of momenta {p} m , one can also show that any solution of eq. (3.4) is one of the solutions of a polynomial equation 5 of degree 2 2m−1 , indicating that the general case requires a numerical solution. In any case, it is immediate to see that if the target masses do not change,m i = m i , one finds κ = 1 and the mapping reduces to the identity. We have so far provided a procedure to generate m momenta with masses {m} m from m momenta with masses {m} m such that their total momentum Q is left unchanged. Let us now show how this affects the phase-space measure. For now we formulate the problem as a change of variables in the m-particle phase space integration, and we will cast this general approach to specific cases of phase-space factorisations for subtraction in the next section.
Working in the rest frame of Q, the original m-particle phase space reads In order to rewrite it in terms of the mapped momenta, it is useful to insert the identity, Changing variables according to p i = κ p i then gives where we introduced the Jacobian of the transformation J. We have checked this result by numerically computing the phase-space volume obtained by integrating over the two phase-space parametrisations for arbitrary choices of masses and up to m = 7 particles. This result is formulated in terms of non-manifestly covariant quantities, but since we worked in the rest frame of Q we can write in order to restore explicit Lorentz covariance.

Mapping of multiple clusters of collinear particles
Let us now see how the transformation applies to clusters of particles becoming collinear to each other. For a single cluster of k massless momenta within an n-particle phase space, By virtue of our ability to map several massive momenta to massless momenta and the easy generalisation of eq. (2.18) to several splittings, we can provide a phase-space factorisation that suits the subtraction of N clusters of k 1 , . . . , k N particles becoming collinear as depicted in fig. 1. The expression is as follows, As a specific example, this allows us to express the qq g g g phase space Φ(p q , pq, p g 1 , p g 2 , p g 3 ) in a factorised way of the form dΦ(qqg 3 ) × dΦ(q → q g 1 ) × dΦ(q →q g 2 ) in order to define counterterms for the "double-collinear" limit where gluon g 1 becomes collinear to the quark and gluon g 2 becomes collinear to the antiquark. This case arises in double-real radiative corrections to the final state qq g at NNLO.

Commutativity and associativity
Processes where multiple disjoint clusters of massless particles can become collinear at the same time feature multiple singular limits that require subtraction. As mentioned above, one singular kinematic configuration is the one where all the children particles become collinear to their respective parents at once. However, the limits where only one or some  children clusters become collinear are also divergent and need to be subtracted, and a pattern of cancellation between the different counterterms is required for the subtraction to work. As a result, it is useful to have a mapping that ensures that the different counterterms for these sub-limits yield reduced phase-space points that match under appropriate conditions. This is guaranteed to happen if the properties of associativity and commutativity are respected by the mapping.
• Commutativity is realised if mapping a process with multiple separate simultaneous splittings sequentially (i.e. splitting by splitting) yields the same reduced phase-space point independently of the order chosen, as illustrated in fig. 2a.
• Associativity is realised if mapping a single splitting with multiple children in one step or sequentially merging subsets of the children yields the same reduced phasespace point, as illustrated in fig. 2b. Let us see how these properties are realised in the case of the generalised rescaling mapping.
Commutativity is easy to prove. Take N clusters of particles whose sums of constituent momenta {p 1 , . . . , p N } are mapped one after the other to on-shell parent momenta recoiling against m other momenta {q 1 , . . . , q m }. Let σ be the permutation of 1, . . . , N that specifies the order in which clusters are merged into their parents. Each step leaves the direction of the clusters' or parent 3-momenta unchanged and rescales them by a parameter κ σ i , where i labels the step . Whatever the order σ of the sequential cluster merging, the final mapped phase space {p 1 , . . . ,p N ,q 1 , . . . ,q m } verifies the following energy-conservation equation, As we already argued, the left hand side is a monotonous function of κ σ = κ σ 1 . . . κ σ N over R + , so that there is a unique physical solution for the final phase-space point independently of the order of the iterated mappings σ. By the same argument, the result is also independent of whether multiple clusters are merged simultaneously or one after the other. The proof for associativity follows along the same line. Without loss of generality, let us consider for simplicity the case of three momenta {p 1 , p 2 , p 3 } with masses {m 1 , m 2 , m 3 } being combined into a momentump 123 with mass m 123 . If we perform the mapping in one step, the vector direction p 123 is kept unchanged and we need to solve for κ 123 in (3.14) If we first map the momenta p 2 , p 3 into an intermediate momentump 23 with mass m 23 , and then mapp 23 withp 1 , we have where κ 1,23 κ 23 must verify the same energy conservation condition as κ 123 and the final parent momentum spatial direction is still that of p 123 . As a result, p 123 = p 123 . Furthermore, any other momentum q j in the process is mapped by having its spatial components rescaled, yielding the same final momentum as well. The result of the mapping is thus independent of whether the particles of a cluster are merged all at once or one after the other.
3.4 Application to e + e − → qq q q g Associativity and commutativity make the subtraction of iterated limits in schemes without sectors considerably more straightforward. Let us illustrate this with the example of a double-unresolved limit of e + e − → qq q q g, where q and q are quarks of different flavour, andq andq are their respective antiquarks. We number the particles as q 1 ,q 2 , q 3 ,q 4 , g 5 . The squared amplitude |M| 2 features, amongst others, two singularities in the limits where either quark-antiquark pair becomes collinear and the other has generic kinematics, C 12 and C 34 . These divergences need to be regulated if integration is to be performed in 4 space-time dimensions. As discussed above, this can be achieved using the factorisation properties of the squared amplitude and a momentum mapping to build counterterms as follows, P 12 |M| 2 (p 12 ,p 3 ,p 4 ,p 5 ) approximates |M| 2 (p 1 , p 2 , p 3 , p 4 , p 5 ) in C 12 , P 34 |M| 2 (p 1 ,p 2 ,p 34 ,p 5 ) approximates |M| 2 (p 1 , p 2 , p 3 , p 4 , p 5 ) in C 34 , where P ij is the splitting kernel for the appropriate limit and we omit spin-correlation indices. Tilded and hatted momenta indicate the mapped momenta in the mappings for the limits C 12 and C 34 respectively.
In the limit C 12,34 where both quark-antiquark pairs are collinear to each other, the counterterms designed for the collinear configurations C 12 and C 34 both approximate the matrix element thus leading to over-subtraction. Moreover, each of these two counterterms for singular single-unresolved configurations in turn presents a divergence in the region of phase space where the opposite mapped quark-antiquark pair goes collinear, i.e.p 3 p 4 andp 1 p 2 . We denote the limits which approach these singular kinematics with C34 and C12. Note that in general the loci of the limits C34 and C12 and those of C 34 and C 12 do not coincide. In order to regulate the divergences associated to these double-unresolved configurations, a counterterm for the limit C 12,34 and counter-counterterms for the limits C12C 34 and C34C 12 need to be introduced,        P 12 P 34 |M| 2 (p 12 ,p 34 ,p 5 ) approximates |M| 2 (p 1 , p 2 , p 3 , p 4 , p 5 ) in C 12,34 , P 12 P34|M| 2 (p 12 ,p 34 ,p 5 ) approximates P 12 |M| 2 (p 12 ,p 3 ,p 4 ,p 5 ) in C34, P12P 34 |M| 2 (p 12 ,p 34 ,p 5 ) approximates P 34 |M| 2 (p 1 ,p 2 ,p 34 ,p 5 ) in C12.

(3.18)
Starting from the singularities which correspond to the two single-unresolved configurations C 12 and C 34 , we obtained an integrand which contains six terms: the original squared matrix element and the five counterterms of eqs. (3.17) and (3.18). This situation is illustrated in fig. 3a. In the limit C 12,34 both quark-antiquark pairs are collinear, and as a consequencep 3/4 → p 3/4 andp 1,2 → p 1/2 . Therefore we also asymptotically havep 3 p 4 andp 1 p 2 , all mappings reduce to the identity and in the exact limit we have The three terms in eq. (3.18) have matrix elements that are evaluated for the same phasespace point which makes it possible for simple cancellation patterns to take place. However in one of the single-unresolved limit, say C 12 , we only findp 3/4 → p 3/4 but in general p 1,2 = p 1/2 . The mapping C 34 need not reduce to the identity and neither does C34. This means that the three terms in eq. (3.18) contain matrix elements (and in general measurement functions) that are evaluated at different phase-space points and it is highly non-trivial for cancellations to occur. This observation alone does not exclude that there might be a way to subtract all overlaps with a clever choice of counterterms, but it is clear that non-commutativity makes engineering iterative counter-counterterms a highly non-trivial task. A much simpler situation can be achieved using a commutative mapping. Indeed, commutativity ensures that the iterated counter-counterterms for the C 12 and C 34 limits have the same reduced kinematics as the C 12,34 counterterm for any phase-space point, as illustrated in fig. 3b. This is for example exploited in CoLoRFul subtraction [43][44][45][46][47][48][49][50][51][52], antenna subtraction [34][35][36][37][38][39][40][41][42] and local analytic sector subtraction [63,64]. On the other hand, the Lorentz mapping, which one could hope to use as an alternative for the rescaling mapping of CoLoRFul for massive final states, is neither commutative nor associative when there are more than two particles in the Born final state.

Jacobians
The phase-space factorisation of eq. (3.12) calls for a discussion of two potential challenges related to the Jacobian of the mapping J: |M| 2 (p 1 , p 2 , p 3 , p 4 , p 5 ) P 12 |M| 2 (p 12 ,p 3 ,p 4 ,p 5 ) P 34 |M| 2 (p 1 ,p 2 ,p 34 ,p 5 ) P 12 P 34 |M| 2 (p 12 ,p 34 ,p 5 ) P 12 P34|M| 2 (p 12 ,p 34 ,p 5 )  • J is a process-dependent function of the phase space. In fact, while for some other mappings such as the ones used in dipole subtraction J is only a function of a fixed number of momenta, here J is a function of all the momenta in the process, meaning that integrated counterterms would be process dependent.
• The Jacobian can be obtained as an explicit function of kinematics only for simple final states, since the degree of the equation that yields κ increases with multiplicity. This is an issue for analytic integration over s K , which would require knowing the full dependence of κ on this variable. The first issue is not noted here for the first time: it was already observed in the case of the rescaling mapping used in CoLoRFul subtraction for final-collinear limits. In that case, the Jacobian features an exponent which depends on the multiplicity, as can be seen in eq. (2.28). A simple but efficient solution was already proposed in Ref. [46] and exploits the fact that the Jacobian has to reduce to 1 in the corresponding collinear limit. As a result any working counterterm for that limit can be divided by J without spoiling the subtraction, since it yields the same leading behaviour. The same solution can be used for this mapping, and also solves the second issue raised above. We will come back to this point in section 5, where the implications of dividing the counterterms by the Jacobian are unfolded further.
3.6 Special cases and relations to other mappings 3.6.1 Rescaling mapping One can show that the generalised rescaling mapping reduces to the rescaling mapping presented in section 2.2.2 when all mapped momenta are massless, i.e.m i = 0 for all i. In that case, the rescaling parameter κ of the generalised rescaling mapping is given in closed form by where we have defined Note that α i is essentially the scalar product of p i with the unit light-like vector n = (1,p i ) and that it is zero for all momenta with m i = 0, so that the sum on the right-hand side of eq. (3.20) effectively runs over the massive parents that are mapped to massless momenta. The complete mapping then reads for massless recoilers. (3.22) Under the same assumptions the Jacobian (3.9) collapses to (3.23) The expressions presented above are slightly more general than those discussed in eq. (3.20) as they handle the case of multiple clusters of particles becoming collinear to each other simultaneously. It is easier to observe the correspondence with the existing literature when looking at specific examples, as discussed below.
Rescaling mapping for a single collinear cluster In the even simpler case of a single final-state collinear cluster of momenta p K = p m + · · · + p m+k−1 with multiple massless recoilers p i , we find This mapping was used for a collinear pair to formulate the CoLoRFul scheme at NLO in [66], and later to subtract triple-collinear limits in [44].
Rescaling mapping for two collinear pairs In the case of two collinear pairs {k 1 , k 2 } and {l 1 , l 2 } and only massless recoiling momenta p i , the rescaling mapping reduces tõ This is the expression used to handle two final-state collinear pairs of partons with an arbitrary number of massless recoilers r in CoLoRFul subtraction at NNLO [44].
Rescaling mapping for one collinear set and one recoiler In the case of a single collinear set K and a single massless recoiler r, we have Q = p r + p K and using p 2 r = 0 we find If the collinear set K is a pair of massless particles {i, j} the expression further simplifies to which is the mapping adopted for dipole subtraction at NLO [30].

Mapping to back-to-back kinematics
When the generalised rescaling mapping is applied to exactly two momenta p 1 and p 2 , conservation laws enforce that in their centre-of-mass frame p 1 + p 2 = 0 and p 1 + p 2 = 0. It is then easy to see that The rescaling parameter is just the positive solution of , (3.30) where λ indicates the Källen function. The full mapping reads and it is straightforward to see that eq. (3.9) reduces to In the case of one collinear cluster and one recoiler, m 2 =m 2 , this transformation corresponds to the momentum mapping used to subtract quasi-collinear singularities in the dipole formalism, when emitter and spectator are both final-state massive particles [67,70,71]. We have checked that eq. (3.32) agrees with Refs. [30] and [67] under this assumption. Finally, we note that there is a unique solution forp 1 andp 2 such that they are in the (p 1 , p 2 ) plane, which here is the same as the (p 1 , Q) or (p 2 , Q) planes. Observing that the Lorentz transformation (2.30) also ensures that the mapped momenta are in this plane, we conclude that for one collinear cluster and one recoiler the Lorentz mapping [68], the generalised rescaling mapping and the dipole mapping [67] are all identical.

Soft mappings with massive recoilers
The phase-space factorisation (2.41), which corresponds to the soft mapping of section 2.4 does not work in the presence of massive final-state particles. The reason is that the rescaling in eq. (2.39) does not let us trade p i forp i in δ(p 2 i − m 2 i ). There are possible workarounds for this issue. An example of such a fix would be to use the generalised rescaling transformation of section 3.1 to map all momenta onto the light cone, then apply the soft mapping of section 2.4 and finally use another generalised rescaling transformation to restore the appropriate on-shell conditions. However, the mappings (2.39) and (2.42) also cannot be used when there is a single massive resolved particle in the final state. Consider, for example, the real-emission corrections to inclusive Higgs production at hadron-hadron colliders via gluon fusion at NLO, where the final state is gluon plus Higgs, p g + p H . Then the reference vectors of the Lorentz transformation would be K = p H /λ andK =p H . The constraint K 2 =K 2 would imply thatp 2 H = p 2 H /λ 2 , which cannot be fulfilled by an on-shell Higgs. A solution that lifts both issues mentioned above consists in avoiding to rescale the hard momenta in the final state, and rescale the total momentum instead. Namely, we modify the single-soft mapping of section 2.4 to m partons of momenta p 1 , . . . , p m and masses m 1 , . . . , m m and one soft gluon of momentum p s , by transforming all of the recoilers' momenta as,p where Λ is given by eq. (2.30) withK = λ s Q and K = Q − p s . The constraint K 2 =K 2 is still given by eq. (2.40), however momentum conservation becomes Instead of a phase space factorisation as in eq. (2.41), we get a convolution, where Λ is given in eq. (2.30), withK = λ s 1 ...sr Q and K = Q − r j=1 p s j . The constraint K 2 =K 2 is given by eq. (2.43) which fixes λ s 1 ...sr . Momentum conservation becomes The phase space is given by the convolution (4.6)

Mapping independence of integrated counterterms
Multiple momentum mappings are suitable to define counterterms that cancel the divergences of real-emission processes. Different choice of mappings yield distinct phase-space factorisation and therefore, a priori, require redoing the integral over the unresolved degrees of freedom which yields the local cancellation of poles with virtual corrections. The poles themselves are universal, but the finite part of the integrals depend on the choice of mapping.
In this section, we argue that local counterterms may be defined in such a way that the analytic form of the integrated counterterms only depends on the mapping through the integration bounds, allowing them to be used for several choices of momentum mappings.
For the sake of simplicity and concreteness, we shall discuss how this can be achieved in the case of the counterterm subtracting the divergence from a single set of k final-state momenta {p} k = {p m , . . . , p n } becoming collinear in a n-particle phase space p 1 , . . . , p n with n = m + k − 1. We refer to their parent momentum as p K = p m + · · · + p n and to momenta before the splitting as {p} m = {p 1 , . . . , p m−1 , p K }. We exploit collinear factorisation of the squared amplitude to define a counterterm for this limit as follows, where P αβ is a splitting kernel, α and β are spin-correlation indices, and |M| 2 αβ is a shorthand for the spin-correlated reduced squared amplitude. As we already discussed in section 2.1, the reduced matrix element is only on shell exactly on the limit, where p 2 K = 0, so we can define our counterterm over the full phase space as where {p} m are mapped momenta and f is a function whose limit equals 1 when p m · · · p n .
The contribution of this counterterm to the total cross section is obtained by integrating over the real-emission phase space, which we first factorise as in eq. (2.18), where P αβ is the integrated kernel, J is the Jacobian of the change of variable ands 0 is the new bound of the virtuality integral after the mapping, which can a priori be a function of the mapped phase-space variables. This occurs, for example, in the case of the Lorentz mapping described in section 2.2.3. As anticipated in section 3.5, the choice f = J −1 appears to be particularly convenient. In this case the condition that f → 1 in the collinear limit is clearly respected, since the mapping has to become a trivial transformation in that region. The only mapping dependence is then contained in the expression of the upper bounds 0 of the virtuality integral, so that computing the integral P αβ analytically as a function ofs 0 permits the usage of the same integrated counterterm for different choices of mappings. Note that for this to be true, the integrand must be the same function of variables of dΦ k and dΦ m independently of the mapping, which might require adjusting the definition of splitting function variables in terms of un-mapped momenta in a mapping-dependent way.
We have verified this assertion on NLO final collinear splittings using a subtraction tool currently under development called MadNkLO. We implemented the NLO CoLoRFul subtraction scheme presented in [48] within this framework, and we slightly modified it to set f = J −1 . We validated our implementation by integrating the NLO real and virtual corrections to 3-jet production in e + e − collisions and comparing to MG5_aMC [72], and we found good agreement within statistical uncertainties both globally and differentially. We then defined a different subtraction scheme by changing the final-collinear mapping to the Lorentz mapping and setting f to the appropriate inverse Jacobian factor expressed in eq. (B.3). We kept the same functional form for the integrated counterterms as in our variation of the CoLoRFul scheme and inserted the corresponding expression ofs 0 . As mentioned in the paragraph above, this implies redefining the energy fraction variable z entering the splitting function such that the integrand has the same expression in terms of mapped momenta as the z used in CoLoRFul, where where α and v are variables of the factorized unresolved phase space. In the CoLoRFul mapping, z i,ij =z i,ij , but not in any other mapping. By choosingz i,ij as the energy fraction, which is a mapping dependent function of the un-mapped momenta, we ensure that the integrand of eq. (5.4) has a mapping independent expression. Both the integrals over the real and virtual corrections are individually affected by this modification, but as expected their sum is not altered by the change of mapping, within sub-percent statistical uncertainties.

Conclusions
In order to obtain more and more precise theoretical predictions for the LHC it is crucial to develop a fully-automated, efficient subtraction algorithm that can provide results at NNLO in QCD, and possibly beyond. Yet "an optimal subtraction method, able to efficiently deal with complex processes has yet to emerge" [73]. An intermediate goal would be to have a subtraction method which, up to the computation of the required two-loop amplitudes, works for every scattering process at NNLO accuracy.
While it is legitimate to push the existing subtraction methods to their maximum computational capabilities, it may be worth to dissect, analyse and compare them with the goal of eventually improving their features and components. The work presented here was inspired by the latter point of view.
In this paper we studied momentum mappings, which are parametrisations of the phase space that factorise the variables that describe the particles becoming unresolved in some infrared or collinear limit from the variables that describe an on-shell phase space for the resolved particles. In sections 3 and 4, we have introduced new momentum mappings for final-collinear counterterms and for soft counterterms. The new mappings work in the presence of particles of arbitrary mass and with an arbitrary number of soft particles or clusters of collinear particles, making them fit for subtraction methods at N k LO accuracy, with arbitrary k. In particular, the new mapping for final-collinear counterterms can also be used to show that at NLO the mappings of Refs. [67] and [68] for massive particles are equivalent in the case of a single recoiler.