Automation of next-to-leading order computations in QCD: the FKS subtraction

We present the complete automation of the universal subtraction formalism proposed by Frixione, Kunszt, and Signer for the computation of any cross section at the next-to-leading order in QCD. Given a process, the only ingredient to be provided externally is the infrared- and ultraviolet-finite contribution of virtual origin. Our implementation, currently restricted to the case of e+e- collisions, is built upon and works in the same way as MadGraph. It is particularly suited to parallel computation, and it can deal with any physical process resulting from a theory implemented in MadGraph, thus including the Standard Model as well as Beyond the Standard Model theories. We give results for some sample processes that document the performances of the implementation, and show in particular how the number of subtraction terms has an extremely mild growth with final-state multiplicity.

Introduction
The vast majority of new-physics searches at the LHC, and some of those at the Tevatron, are based on final states signatures that feature the presence of a large number of jets (typically, from four to ten), generically denoted as multi-jet configurations. The fairly difficult problem of giving theoretical predictions for multi-jet processes has been completely solved in the past few years, but the solutions are only accurate to the leading order (LO) in QCD. All approaches are thus based on the computations of tree-level, multi-parton matrix elements, performed by dedicated computer programmes (see e.g. ref. [1] for a review and a list of references); the LO accuracy implies that the computations can be carried out in four dimensions, and that all divergences resulting from integration over the phase space are avoided by means of kinematic cuts. Through the computation of the matrix elements, one obtains final states composed of quarks, gluons, and other accompanying particles (e.g. W 's or Z's). At this point, one may choose to identify each light quark and gluon with a jet, thus making use of the local hadron-parton duality, and compare the theoretical results to data; we shall call these results matrix-element predictions. A more realistic approach, which is widely used by experimental collaborations, is that of feeding parton-level final states obtained from matrix element computations to event generators, which by showering them eventually result in hadron-level final states, where any jet may contain several tens of particles; we shall call these showered predictions. As is known, showered predictions are only meaningful in the context of a proper matching formalism, which avoids the double counting of configurations that can be obtained both from the matrix elements and from the showers. Several solutions are available for the problem of tree-level matching [2,3,4,5], which involve the simultaneous treatment of final states with different multiplicities.
While fairly successful phenomenologically, the approaches discussed above suffer from the usual limitations of LO computations. Scale uncertainties, which are typically large for processes with high-multiplicity final states, render it unreliable the predictions of the absolute values of the cross sections. A common strategy is therefore that of fixing the normalization equal to the data for a given jet multiplicity, and of predicting the cross sections for larger and smaller multiplicities. It is therefore desirable to extend the accuracy of multi-jet computations to the nextto-leading order (NLO) in QCD. It is in fact worth recalling that NLO results, so far available only for small-multiplicity jet final states, have been a very important ingredient in establishing QCD and the SM as the correct theories of strong and electroweak interactions at LEP, SLD, and the Tevatron, as well as to exclude the presence of Beyond the SM signals in the data so far. We also recall that NLO calculations, apart from resulting in smaller theoretical uncertainties that allow one to trust absolute predictions, induce non-trivial effects in particle spectra that cannot be obtained in general by a simple rescaling of LO results.
The problem of the computation of multi-jet cross sections at the NLO has attracted considerable attention in the recent past. There are two main obstacles that must be cleared. Firstly, one needs to compute the one-loop corrections. There are now several different solutions to this problem, that are reasonably automated, and are based on unitarity methods: BlackHat, CutTools, and Rocket (see refs. [6,7,8] for recent results, and for a more complete list of references). A more traditional approach, based on Feynman diagrammatics, is that of GOLEM [9]. It is reasonable to expect that most of these codes will attain larger flexibility and speed in the near future. Secondly, one must compute the real-emission contributions, and combine them with the one-loop ones in order to obtain the physical cross sections.
The aim of this paper is that of addressing the second of the problems mentioned above. We assume that ultraviolet-renormalized one-loop corrections are given, and we fully automate the following steps: d) output of the results on an event-by-event basis, in the form of a set of fourmomenta and a weight, that one can use to construct (infrared-safe) observables.
We stress that item b) involves the cancellation of infrared poles as prescribed by the KLN theorem. Even if the explicit form of the one-loop contributions is not available, this task can be carried out analytically, thanks to the fact that the general structure of the poles is known for any process; the residues are proportional to tree-level matrix elements that can be computed automatically.
There are two main ingredients in our work: a programme that computes treelevel matrix elements, and a universal formalism for the analytical cancellation of infrared singularities, that allows one to implement in a computer codes short-distance cross sections free of any singularities. For the former, we use MadGraph [10,11]. For the latter, we adopt the formalism of Frixione, Kunszt, and Signer, presented originally in refs. [12,13], which we shall refer to as FKS in the rest of this paper. Through the implementation of items a)-d) above, we basically achieve an NLO version of MadGraph/MadEvent, up to a missing (infrared-and ultraviolet-finite) piece of one-loop origin, in which the only operation required from a user is that of typing in the process to be computed. We call this programme MadFKS. The current implementation of MadFKS, and the results given in this paper, are relevant only to the case of e + e − collisions, which is sufficient to highlight the capabilities of the code.
On the other hand, the formalism is written in full generality, and can therefore be immediately applied to any other kind of collisions. We defer to a forthcoming paper the complete implementation of the cases where initial-state hadrons are present. We stress, however, that what is done here will apply without any modifications (except for the trivial multiplication by the relevant parton density functions) to those cases as well. Thus, the only new piece of implementation required by initial-state hadrons will concern the emission of a parton from an initial-state parton. We finally point out that the results of MadFKS are matrix-element predictions, but also give the necessary building blocks to arrive at NLO showered predictions in the context of matching techniques such as MC@NLO [14] or POWHEG [15,16].
This paper is organized as follows. In sect. 2 we give a general overview of the issues in NLO computations of multi-jet cross sections, and of the defining features of the FKS subtraction method. In sect. 3 our notations are introduced. In sect. 4 we summarize the most relevant formulae of the FKS subtraction for any type of collisions, including hadron-hadron ones. In sect. 5 we explain how these formulae can be implemented in a computer code. A few refinements which we use in MadFKS to optimize the performances of the code are discussed in sect. 6. In sect. 7 we present selected results of MadFKS, for the case of e + e − collisions. Section 8 reports our conclusions. Further technical information such as longer formulae, and possible variants of the implementation, are given in the appendices.

Computations of multi-jet cross sections at the NLO
Given a production process, the computation of its NLO corrections in QCD implies the evaluation of the corresponding one-loop and real-emission corrections (for a pedagogical introduction to these concepts, see e.g. sect. 4 of ref. [1]), which are eventually added up to get the physical cross section. Although from the principle point of view there is no difference between the computation of a small-multiplicity and that of a large-multiplicity jet cross section, in practice the situation is entirely different. One-loop computations based on straightforward Feynman diagram evaluation rapidly become too involved to be carried out 1 , and new techniques have been developed recently to bypass this problem (see e.g. refs. [18,19,20,21]). On the other hand, the computation of real-emission and Born (i.e., tree-level) diagrams poses no problem, and is limited only by CPU.
There are two issues when summing one-loop and real-emission contributions. The first is that of proving (analytically) the cancellation of the infrared singularities that arise in the intermediate steps of the computation, and to write the leftovers as finite contributions that can be implemented in computer codes. The second is the actual implementation of these finite terms, which crucially includes an integration over the phase space. The former issue has been fully solved [12,22,23,24] in the 90's, for any jet multiplicity (in other words, the jet multiplicity is simply a parameter in the resulting formulae), in the context of the so-called universal cancellation formalisms, which use either the (approximate) slicing method, or the (exact) subtraction method; it is nowadays acknowledged that the slicing method is unsuited for describing complicated final states, such as those in multi-jet production.
The latter of the two issues mentioned above has actually not even been considered, for multi-jet cross sections, until recently. There are now several proposals [25,26,27] that aim at constructing local counterterms for any given real-emission matrix elements (item a) above). What is outlined in items a)-d) is achieved in ref. [28] (restricted to massless SM particles), and very recently in ref. [29] (which includes the treatment of massive SM particles). All of these approaches have tackled the problem using the universal subtraction method originally introduced in ref. [22] (referred to as dipole subtraction henceforth).
As was discussed in the introduction, in this paper we use the FKS subtraction formalism. Although dipole subtraction is the most widely used method for NLO matrix-element computations, FKS is the only technique used so far in practice for performing NLO showered computations, in the MC@NLO and POWHEG frameworks. However, as we shall show in this paper, FKS subtraction has also excellent performances in the context of matrix-element computations. While dipole and FKS subtraction methods are formally equivalent as shown in ref. [30], differences between the two arise because of the way in which the building blocks for the local counterterms of the real-emission matrix elements are combined (these building blocks corresponds to the soft, collinear, and soft-collinear singularities). In the dipole method, one emphasizes the role of soft emissions. A structure arises where one sums over an "emitted" parton, and two colour partners that exchange the emitted parton; the sum over dipoles is thus a sum over three indices. In the FKS method, the emphasis is on collinear emissions, and the structure that emerges is therefore one where there is a sum over parton pairs, i.e. a sum over two indices. The consequence of this is that, in general, the number of independent subtraction terms is smaller in the FKS formalism than in the dipole formalism.
One of the key features of FKS subtraction is that, in the integration of the realemission matrix elements, one effectively defines partonic processes with at most one soft and one collinear singularities, which results in a much simpler subtraction structure than that of the original matrix elements. These processes are furthermore fully independent from each other, and can therefore be computed separately. This is a genuine parallelization, more efficient than that of computing the same production process several times with different seeds for random number generation, especially in view of the use of adaptive-integration routines (the latter type of parallelization can obviously still be set up in the context of FKS subtraction, if need be).
The original FKS papers [12,13] addressed the case of processes with final states composed only of massless quarks and gluons. The extension to generic processes, with the presence of both strongly-interacting massive particles and of non stronglyinteracting particles, is almost trivial. Several of these cases have in fact been implemented in MC@NLO using FKS subtraction (see e.g. ref. [31]). Here, we shall present a complete summary of the relevant formulae, which can also be applied to computing cross sections to NLO accuracy for the production of non-SM particles.
In summary, we advocate the use of FKS subtraction for matrix-element predictions for a variety of reasons: structure identical to that of collinear emissions, thus lending itself naturally to matching with parton shower Monte Carlos not based on colour dipoles; ease of importance sampling; small number of subtraction terms, and the modest growth of their number with the multiplicity, with beneficial effects on numerical stability; organization of the calculation in a way which is parallel in nature; and the availability of all results necessary for the treatment of fully-polarized processes, which allows one to performing sums over helicity states with Monte Carlo methods.
In the rest of this paper, all formulae for NLO cross sections will make use of FKS subtraction.

Partonic processes
For any given partonic process, we shall denote by n the number of final-state particles at the Born level. The contributions to the NLO partonic cross sections will therefore either have a 2 → n or a 2 → (n + 1) kinematics, which we shall call Born (or n-body) or real-emission (or (n + 1)-body) kinematics respectively 2 . They will be denoted as follows: Born kinematics ; (3.1) We shall write the corresponding phase spaces as respectively. We shall denote by n ∅ the number of final-state particles which are not strongly interacting (e.g. leptons); n H will denote the number of massive, stronglyinteracting particles (e.g. heavy quarks). Depending on the process, some or all of the n ∅ + n H particles will not belong to the Standard Model. n-and (n + 1)-body processes will be characterized by the same values of n ∅ and of n H ; the identities of the corresponding particles will also be the same in the two classes of processes. On the other hand, the number of light quarks and gluons in the (n + 1)-body processes (denoted by n (R) L ) will be equal to that in n-body processes (denoted by n (B) L ), plus one. We have therefore: If 1 ≤ k ≤ n + 2 or 1 ≤ k ≤ n + 3 is an index which runs over all particles, we shall adopt the following labeling scheme: massless quarks and gluons ; non strongly-interacting particles . (3.10) In order to deal with the various types of colliding particles, it will also be useful to introduce the quantity n I , which can take the following values n I = 1 ⇐⇒ HH collisions , (3.11) n I = 2 ⇐⇒ eH collisions , (3.12) where e and H generically denote a non-hadronic and a hadronic particle respectively; n I is therefore the smallest value that a label can assume when running over light quarks and gluons (without loss of generality, we have conventionally given label 1 to the non-hadronic particle in an eH collision). The identity of particle k will be denoted by I k ; it is also convenient to denote the anti-particle of I k by Ik, although I k can be used as well. A process will be unambiguously identified by giving the list of the identities of the particles involved, and we shall denote it by r = (I 1 , . . . I n+2 ) or r = (I 1 , . . . I n+3 ) . (3.14) The sets of all n-and (n + 1)-body processes will be denoted by R n and by R n+1 respectively 3 . In the following, we shall make use of the fact that a given n-body 3 If two processes in the form of eq. (3.14) are related by a permutation over the identities of final-state particles, only one of them will be included in R n or R n+1 . process can be related to at least one (n+1)-body process. We therefore introduce the following notations, where i is a final-state light quark or a gluon (3 ≤ i ≤ n (R) L + 2), and j is a strongly-interacting particle (n I ≤ j ≤ n (R) L + n H + 2): We have defined Thus, process r i \ in eq. (3.17) is constructed by simply removing parton i from the original list. It is clear that such an operation results in a physical process (i.e., which exists) only when i is a gluon. Process r j⊕i,i \ in eq. (3.18) is constructed by removing parton i from the original list, and by replacing particle j with one whose identity is that of the parton entering the (I j⊕i , I j , I i ) QCD vertex (thus, if I i = g and I j = g, then I j⊕i = g; if I i = g and I j = q, then I j⊕i = q, and so forth). If this QCD vertex does not exist (as e.g. in the case I i = u and I j =d), process r j⊕i,i \ is non physical, i.e. it does not exist 4 . The definitions of r i \ and of r j⊕i,i \ are obviously motivated by the relevant soft and collinear limits respectively, but will be used in several different contexts 5 .

Matrix elements
NLO cross sections receive contributions from both tree-level and one-loop amplitudes, which we shall denote as follows:

19)
A (n,0) (r) −→ Born tree amplitude (r ∈ R n ) , (3.20) The amplitudes above are typically relevant to a given spin and colour configuration; they include all the coupling constant factors. Starting from amplitudes, one needs to construct several quantities which will enter short-distance cross sections. We denote them as follows: In these equations, s = (k 1 +k 2 ) 2 = 2k 1 ·k 2 , and ω(I) is the product of spin and colour degrees of freedom for particle I. In 4 − 2ǫ dimensions, we have ω(q) = 2N c ≡ 6 and . These average factors serve the sole purpose of fully specifying the divergent part of the one-loop contribution in the CDR scheme, which is not used in numerical calculations (see app. B for more details). We stress that the formulae which give the physical cross sections in FKS, and that are implemented in computer codes, are finite in four dimensions. Note that symmetry factors accounting for the presence of identical particles in the final state are not included in the matrix elements defined above; they will be inserted later in the expressions for the shortdistance cross sections. As the notation suggests, r ∈ R n+1 in eq. (3.22), while r ∈ R n in eqs. (3.23)- (3.25). The quantities defined in eqs. (3.22)-(3.25) are sufficient to perform all of our computations, with the exception of an azimuthal contribution to the collinear limit of M (n+1,0) , which requires a reduced matrix element that we shall denote by M (n,0) . The precise form of this quantity is irrelevant in this section; it will be used in sect. 5.3, and its definition will be given in app. D. The colour operators Q(I) that enter the definition of the so-called colour-linked Born's, eq. (3.24), give a representation of the colour algebra, and can be defined as follows: with t a and T a the SU(3) generators in the fundamental and adjoint representations respectively. Clearly, eq. (3.26) can be extended to higher-dimensional representations if need be. It is trivial to show that where Note that M (n,0) kl is an interference term related to the exchange of a soft gluon between particles k and l, which appears in the soft limit of the real matrix elements. This term has a coefficient proportional to an eikonal factor, which therefore vanishes identically when k = l and k is massless (see app. A).
The final formulae of the FKS subtraction method, to be given below, will only deal with non-divergent quantities. The interested reader can find in the original paper (ref. [12]) the proof of the cancellation of the infrared and collinear singularities that arise in the intermediate steps of the computation. As in all computations at the NLO, the finiteness of the partonic short-distance cross sections is in general a consequence of imposing kinematic cuts on final-state particles. Without loss of generality, we can always assume these cuts to be equivalent to the request of having either n (B) L or n (B) L + 1 jets in the final state, the jets being reconstructed with an arbitrary algorithm. This condition will be symbolically denoted by (3.36) The condition in J n (B) L is sufficient to prevent the appearance of phase-space singularities 6 in n-body quantities such as M (n,0) , M (n,0) kl , and M (n,1) . On the other hand, the (n + 1)-body tree-level matrix elements will still diverge in some regions of the phase space; these divergences are subtracted (in all subtraction methods) by means of suitable counterterms. In order to classify these divergences and to eventually 6 It is understood that the condition in J n (B) L also prevents all divergences due to final-state leptons and photons. This is necessary if the subtraction of the corresponding QED singularities is not carried out. subtract them, for any given process r ∈ R n+1 we introduce the following set of ordered pairs (called the set of FKS pairs) In words, a pair of particles belongs to the set of FKS pairs if they induce soft or collinear singularities (or both) in the (n + 1)-body matrix elements, even in the presence of jet cuts. The first element of the pair will be called the FKS parton, and the second element will be its sister 7 . Conversely, by definition P FKS takes into account all phase-space singularities of the (n + 1)-body tree-level matrix elements after jet cuts. We note that the condition k j → 0 is simply not relevant if j = 1, 2, since in such a case the momentum k j is fixed. In the computation of an NLO cross section according to the FKS formalism, each pair belonging to P FKS will correspond to a set of subtractions of soft and collinear singularities which, when combined with the real-emission matrix element, will result into a finite contribution to physical observables. Each of these contributions is separately finite, and can therefore be computed independently from the others, and integrated with numerical methods. It is immediately obvious that some of the pairs (i, j) with i, j ≤ n (R) L + n H + 2 will not belong to P FKS for any physical process. This is the case e.g. when i and j are both massive and strongly interacting; or when i and j are a quark and an antiquark of different flavours; or when one of the members of the pair is massive and strongly interacting, and the other member is not a gluon. In general, however, for large multiplicities the number of elements in P FKS will scale approximately as n (R) L (n (R) L + n H ). On the other hand, depending on the identities of the particles in r, some of the contributions due to the various FKS pairs will actually be identical, because of the symmetry properties of matrix elements and phase spaces. This implies a drastic reduction of the number of independent terms actually needed in the computation. We shall return to this issue in sect. 6.1.

Cross sections
We write an NLO cross section in the collision of two particles P 1 and P 2 with momenta K 1 and K 2 using the factorization theorem (4.1) As the notation suggests, we have (see eqs. (3.1) and (3.2)) will be a parton density function (PDF). Otherwise, it may simply be equal to δ(1−x α ) (e.g. when describing an electron beam at a fixed energy), or to a more complicated function that effectively describes the energy loss when the incoming particle P α turns into the particle I α that enters the hard reaction (e.g., f (Pα) Iα may be the Weizsäcker-Williams function with P α = e and I α = γ). These details are in any case irrelevant in what follows, where we shall deal only with the short-distance partonic cross sections dσ (n) , dσ (n+1) , and dσ (n+1) ; if necessary, we shall refer to functions f (Pα) Iα generically as PDFs. In this paper, we shall denote by µ the common value of the factorization and renormalization scales, µ = µ F = µ R , which simplifies the writing of the formulae. At the end, one is always able to recover the separate dependence upon µ F and µ R , by exploiting the renormalization group invariance with respect to these two scales. All the relevant formulae are given in app. C.

n-body contributions
We decompose the n-body cross section into four terms: (4. 3) The first term on the r.h.s. of eq. (4.3) is the Born contribution (r ∈ R n ) 4) with N (r) the symmetry factor that takes into account the presence of identical particles in the final state. The remaining three terms on the r.h.s. of eq. (4.3) have, roughly speaking, a collinear, soft, and one-loop origin respectively; they are separately finite, and we understand the kinematics given in eq. (3.1). We have: where Note that, in this equation, only light quarks and gluons are involved. We have s = 2k 1 · k 2 , and E k is the energy of parton k in the c.m. frame of the incoming particles I 1 and I 2 (i.e., the partonic c.m. frame in the case of hadronic collisions). The Casimir's C(I) have been defined in eqs. (3.29) and (3.30); the other colour factors are: Equation (4.6) contains two free parameters, ξ cut and δ O . The same parameters will be used in the subtraction of the soft and final-state collinear divergences, respectively, that affect the real-emission contribution, as it will be described below. Thus, although eq. (4.5) does depend on ξ cut and δ O , the physical cross section does not. This constitutes a powerful check of the correctness of the implementation of the formalism; we shall return to this point in sect. 7. In eq. (4.6) we have denoted by Q the Ellis-Sexton scale [32], which is an arbitrary mass scale that one may use to write the one-loop results; more details on its role will be given in app. B and in app. C. Here we just stress that the NLO cross section is exactly independent of Q. Finally, we remind the reader that µ = µ F = µ R . Note that the contribution in round brackets in the first line on the r.h.s. of eq. (4.6) is non trivial only if at least one of the incoming particles is a hadron. Still, the form in eq. (4.6) applies to any kind of collisions, provided that one defines The n-body contribution of soft origin reads: The quantities E Finally, the n-body contribution of one-loop origin is: (4.14) The quantity V (n,1) F IN is the non-divergent part of M (n,1) introduced in eq. (3.25). Its definition is therefore unambiguous only in a definite scheme, and the one to be used in eq. (4.14) is the Conventional Dimensional Regularization (CDR) scheme. The extraction of V (n,1) F IN from a complete one-loop computation, and its relation to the popular Dimensional Reduction scheme, is discussed in app. B.

(n + 1)-body contributions
The short distance cross sections dσ (n+1) in eq. (4.1) are due the contributions of the real-emission matrix elements M (n+1,0) , with their phase-space singularities suitably subtracted. As is known, the integration of the subtracted matrix elements is the most involved from the numerical viewpoint, and the difficulty increases with the number of external legs, because of the proliferation of singularities. In the FKS formalism, this problem is simplified by effectively partitioning the phase space, in such a way that in each of the kinematic regions resulting from the partition at most one soft and one collinear singularity are present. This partition is achieved by introducing a set of positive-definite functions 8 where the argument r ∈ R n+1 stresses the fact that different choices for S can be made for different processes; this property is useful to reduce the computing time.
In the following, we shall often understand the argument r in the S functions. There is ample freedom in the definition of the S functions, but the following constraint must be obeyed: (4.16) Furthermore, these functions must have the following behaviours in the kinematic configurations associated with the singularities of M (n+1,0) : S ij = 0 ∀k with I k = g and ∃l with (k, l) ∈ P FKS or (l, k) ∈ P FKS . (4.20) In words, S ij goes to zero in all regions of phase space where the real emission matrix elements diverge, except if this involves particle i being soft, or particles i and j being collinear. Note that eqs. (4.17) and/or (4.19) need not be satisfied if one of the partons in the relevant FKS pair is massive; in such cases, the limits can assume arbitrary values. Also, eq. (4.20) implies in particular that the limit of S ij is zero when j is soft; thus, although according to eq. (3.37) there may be soft singularities associated with FKS sisters, they are damped by using S ij . The functions h ij (z) introduced in eq. (4.17) are defined in 0 ≤ z ≤ 1, and have the following properties: The physical cross section is independent of the choice of h(z); we shall return to this point in sect. 7. As stressed above, the particular form of the S functions is not important in what follows; we shall give an explicit construction in sect. 5.2. As a final remark, we point out that eqs. (4.21) and (4.22) imply that the notation h ij is redundant here; it will however become useful in the context of the optimization we shall carry out in sect. 6. Using eq. (4.16), one proceeds by writing Thanks to eqs. (4.19) and (4.20), each term on the r.h.s. of eq. (4.24) will only be singular when particle i is soft, and/or particles i and j are collinear 9 . Therefore, each term can actually be regarded as describing a production process with the simplest possible structure of phase-space divergences. Furthermore, these contributions are fully independent from each other.
We shall now give the expressions of the subtracted real-emission cross sections. In order to do so, we work in the c.m. frame of the incoming particles: In this frame, for each pair (i, j) ∈ P FKS we introduce the variables ξ i and y ij , where In other words, ξ i is the rescaled energy of the FKS parton, and y ij is the cosine of the angle between the FKS parton and its sister. Using these variables, the soft and collinear singularities of S ij M (n+1,0) correspond to ξ i = 0 and to y ij = 1 respectively (note, however, that the matrix element may not be singular in either or both of these limits, depending on the identities of particles i and j). We have (see ref. [12]) The variable ϕ i parametrizes the azimuthal direction of the FKS parton, but its precise definition is irrelevant here. The quantity d φ ij n is related to the n-body phase space, and is implicitly defined as follows: The actual form for d φ ij n depends on the specific parametrization adopted for dφ n+1 , but the following equations must always hold: The distributions entering eq. (4.29) are defined as follows 10 : In eq. (4.33) and (4.34) ξ cut and δ are free parameters, that can be chosen in the ranges It would be possible to choose different values of ξ cut and δ for each of the cross sections dσ (n+1) ij contributing to eq. (4.28). In practice, we shall choose the same value of ξ cut for all (i, j) pairs, the same value of δ for all pairs (i, j) with j ≥ 3 (which we shall denote by δ O ), and the same value of δ for all pairs (i, j) with j ≤ 2 (which we shall denote by δ I ). The parameters ξ cut , δ O , and δ I are therefore associated with soft, final-state collinear, and initial-state collinear singularities respectively. As was discussed in sect. 4.1, the physical cross section on the l.h.s. of eq. (4.1) is strictly independent of ξ cut , δ O , and δ I , while this is not the case for the short-distance cross sections dσ (n) , dσ (n+1) , and dσ (n+1) .
For the reader unfamiliar with the use of plus distributions in the context of cross section calculations, it is useful to explicitly show the structure of eq. (4.29). We start by rewriting such an equation with the following shorthand notation: We can now expand the plus distributions using their definitions, eqs. (4.33) and (4.34): (4.38) There are four terms in the integrand in eq. (4.38). We shall call them "event" (Σ ij (ξ i , y ij )), "collinear counterevent" (Σ ij (ξ i , 1)), "soft counterevent" (Σ ij (0, y ij )), and "soft-collinear counterevent" (Σ ij (0, 1)). In general, they may have four different kinematic configurations, but one can actually reduce them to two, as discussed in sect. 5.1. Note that, in the case of a massive sister (m j = 0), the collinear and softcollinear counterevents vanish identically, because of the damping factor (1 − y ij ) in eq. (4.29).
Finally, consider the case m i = 0. According to eq. (3.37), (i, j) ∈ P FKS only if I j = g. The only possible singularity arising from the FKS pair (i, j) is the soft one k 0 j → 0; this is however damped by S ij . It follows that all counterevents in eq. (4.38) vanish, and the distributions in eq. (4.37) coincide with ordinary functions. Therefore, using eq. (4.30), one simply gets Since there are no subtractions involved, one is not bound to use the variables ξ i , y ij , and ϕ i in the parametrization of dφ n+1 in this equation. On the other hand, there is nothing that prevents one from doing so; in such a case, note that the smallest value ξ i can assume is 2m i / √ s, and we understand in eq.
The fact that when m i = 0 one does not need to subtract any singularities suggests that from the viewpoint of the efficiency of the calculation it is not convenient to define the quantity dσ (n+1) ij in eq. (4.39) as an independent contribution; we shall discuss this issue in sect. 6.

Degenerate (n + 1)-body contributions
These contributions, denoted by dσ (n+1) in eq. (4.1), are the finite remainders left by the subtraction of initial-state collinear counterterms. We write The condition (i, j) ∈ P FKS (with j = 1, 2) in eq. (4.40) implies that the realemission matrix element is singular when k i k j . This typically means that, if P j is not a hadron, the corresponding contribution dσ (n+1) ij need not be included. A notable exception is I i = q and I j = γ, when the matrix element has a (QED) singularity, which is cancelled by the counterterm due to the inhomogeneous part of the (Altarelli-Parisi) evolution equation for the hadronic photon. In any case, by writing (i, 1) ∈ P FKS and (i, 2) ∈ P FKS in eq. (4.40) we can use this equation for any kind of collisions. As shown in ref. [12], one obtains where, in analogy with eq. (4.33), we have introduced the distribution We have also defined the functions i.e. the unregularized Altarelli-Parisi splitting functions times a (1 − z) factor that damps the z → 1 soft singularity. The explicit forms of P (0) ab and P (1) ab are given in app. D. The functions K ab are related to the choice of the PDF scheme; for all practical purposes they are trivial, since in the MS scheme we have K ab = 0. We point out that in the special case alluded to at the beginning of this section, i.e. I i = q and I j = γ, the relevant dσ (n+1) ij must be proportional to α em . We can therefore still use eqs. (4.41) and (4.42), provided that we include a factor α em /α S in the kernels P qγ and K qγ (see ref. [13]).
Consider dσ is fully analogous). The fact that (i, 1) ∈ P FKS implies that r 1⊕ī,i \ exists and belongs to R n , and that I 1⊕ī is either a light quark or a gluon. The notation embodies the fact that the matrix element and phase space must be computed with 2 → n kinematic configurations, where incoming partons have momenta We conclude by stressing that the symmetry factor in eqs. (4.41) and (4.42) is N (r), and not N (r 1⊕ī,i \ ) or N (r 2⊕ī,i \ ). This can be easily understood as follows. Suppose that there are p identical particles in the final state of r, with labels i 1 , . . . i p ; N (r) will thus contain a factor p!. On the other hand, we shall also have dσ (4.46) Apart from allowing one to reduce the computing time by a factor of p, this equation also shows that symmetry factors N (r 1⊕ī 1 ,i \ 1 ) and N (r 2⊕ī 1 ,i \ 1 ) will emerge naturally, since .
This equation obviously holds also for p = 1.

Implementation
The most involved part in the computation of an NLO cross section is the integration of the subtracted real-emission (i.e., (n+1)-body) matrix elements. In this section we shall therefore mainly present technical details on the ingredients that will enter such an integration in our implementation. In particular, we shall discuss the kinematics relevant to the generation of counterevents, give the formulae relevant to the evaluation of the damped matrix elements in the singular limits, and present an explicit construction of the S functions. We stress that the latter is fully general, and easy to automate, but nothing prevents one from using alternative forms, provided they fulfill the conditions given in eqs. (4.16)- (4.20). Although the practical applications in this paper are limited to the case of e + e − collisions, the formulae we give in the following can be applied without modifications to any kind of incoming particles.

Kinematics
The integration of the subtracted real-emission matrix elements over the (n + 1)body phase space entails the generation of 3n − 1 random numbers, of which three correspond to the integration variables ξ i , y ij , and ϕ i , that have been singled out since they are directly associated with the FKS parton. We denote the remaining 3n − 4 random numbers by x α ; we can write a generic 2 → n + 1 kinematic configuration as follows: As implied by eq. (4.38), this configuration will be that of the event. Starting from it, we construct the configurations corresponding to the soft, collinear, and soft-collinear counterevents 11 : with ϕ i and x α in eqs. (5.2)-(5.4) the same as those used in the generation of the kinematics of the event, eq. (5.1). For a given choice of the 3n − 1 independent variables, the momenta k l depend on the phase-space parametrization adopted. We have considered three different such parametrizations, in order to test the convergence behaviour and numerical stability of the numerical integration of short-distance cross sections. All obey the following equations: These properties trivially follow from the observation that the kinematic configurations of the counterevents are degenerate, and effectively correspond to a 2 → n configuration. This configuration, whose momenta we denote [14,16] byk, can simply be defined with a relabeling of the soft counterevent kinematics: with σ denoting the relabeling. We stress that, although eqs. (5.5) and (5.6) need not be fulfilled by any phasespace parametrization, and are not mandatory for the implementation of FKS subtraction, they are quite intuitive from the physical point of view, and help improve the numerical stability of the results for differential distributions, while also reducing the computing time.
We finally point out that the results of this section apply to the case j ≥ 3, i.e. to final-state emissions. The situation is more involved in the case of j = 1, 2, i.e. of initial-state emissions, but similar simplifications are possible (see e.g. app. A.4 of ref. [14]). We postpone a discussion on this issue to a forthcoming publication.

Choice of S functions
As already mentioned in sect. 4.2, in ref. [31] a smooth form of the S function was adopted, which was seen to improve the numerical integration with respect to the form of refs. [12,13], where Heaviside Θ's had been used. In ref. [31], the S functions were constructed using invariants k k ·k l , but given that FKS subtraction uses energy and angles as integration variables, in a further generalization (see e.g. ref. [16]) one can also define: where the functions h ij have been introduced in eqs. (4.21)-(4.23), and The parameters a S and b S introduced in eq. (5.11) are real, positive, and arbitrary. The fact that the physical cross section is independent of their choices constitutes another check of the correctness of the implementation of FKS subtraction. We shall return to this point in sect. 7. As shown in eqs. (4.21) and (4.22), the functions h kl are given in terms of a function h(z), that we have chosen as follows: In practice, the specific form of h(z) has only a modest impact on numerical computations; we shall briefly comment on this point in sect. 7. The role of h kl in eq. (5.10) is simply that of avoiding to count twice the contribution of 1/d kl when both (k, l) and (l, k) belong to P FKS . In fact, in this case one can immediately obtain having used eqs. (4.23) and (5.9). Alternatively, one could have dropped the h kl factor in eq. (5.10), and summed over unordered pairs. We point out that the limits of the S functions given in eqs. (4.17)-(4.20) are related to the vanishing of d ij in the collinear (θ ij → 0 with m i = m j = 0) and soft (E i → 0) limits. When parton i is not a gluon, the soft limit is simply not relevant (the matrix element is non singular), and the S functions can assume arbitrary values. This implies that, in such a case, we may replace 2E i / √ s with 1 in the definition of d ij . Along the same lines, one may always replace β i and/or β j with 1 in d ij when particle i and/or j is massive; in doing so, d ij will not be proportional any longer to k i ·k j when a S = b S , but this is perfectly acceptable. We have tested these options, and will comment on this point in sect. 7.
The S functions defined in eq. (5.8) are suitable for numerical implementation, since they are well behaved in the whole (n + 1)-body phase space, and in particular in the soft and collinear limits. In fact, after trivial algebra (also using eq. (5.14)) one gets The E i dependence can be analytically removed in the second term on the r.h.s of eq. (5.15), and therefore Dd ij can be computed numerically both in the soft and collinear limits, which are needed for the evaluations of the soft, collinear, and softcollinear counterevents.

(n + 1)-body matrix elements
In the previous section we have seen how to compute the S functions in eq. (4.29) in the various singular limits. The other major ingredient for the evaluation of eq. (4.29) is the damped matrix element We shall now give explicit, finite, expressions for M (n+1,0) ij in the soft and collinear limits. Using the notation introduced in eq. (5.7), we have where σ is the relabeling of the indices introduced in eq. (5.7). Note that the eikonals on the r.h.s. of eq. (5.17) are finite at ξ i = 0, since the factor ξ 2 i in the numerator is cancelled by that resulting from E 2 i in the denominator (see eq. (4.26)). Furthermore, in the case m j = 0, the collinear limit of eq. (5.17) is also finite, since the factor 1 − y ij in the numerator cancels an identical factor in the denominator, when p = j or l = j (see eq. (4.27)); this collinear limit can thus be used to construct the softcollinear counterterm. The cancellations of the ξ 2 i and 1 − y ij factors are carried out analytically, and eq. (5.17) is therefore well behaved numerically in the whole (n + 1)-body phase space.
In the case m i = m j = 0, we also need to compute the collinear limit. We have (see ref. [12]) where z ji has been defined in eq. (5.9), and the functions Q ab ⋆ (z) and reduced matrix elements M (n,0) are given in app. D. The second term on the r.h.s. of eq. (5.18) vanishes when the integration over the azimuthal angle ϕ i is carried out. However, it is in general non zero pointwise in the phase space and, as eq. (5.18) shows, it actually has the same behaviour for y ij → 1 as the first term; it must therefore be included when constructing a local collinear counterterm. Using eqs. (4.26) and (4.27), eq. (5.18) becomes: where, in analogy to eq. (4.44), we have defined It is apparent that the r.h.s. of eq. (5.19) can be computed numerically also in the soft limit, E i → 0 (which is equivalent to z ji → 1), owing to the regularizing factors (1 − z) of eqs. (4.44) and (5.20). In the cases (I i , I j ) = (g, g) and (I i , I j ) = (q, g) (where q is a massless quark) the Altarelli-Parisi kernels have a soft singularity for z ji → 0. However, this singularity is eventually damped, in the computation of the cross section, by the factor h ij (z ij ) in eq. (5.8), and thus need not be subtracted.

Optimization
What we have done so far allows one to implement the FKS subtraction method in a computer code. We shall now tackle the problem of how to make this straightforward implementation more efficient. As discussed in sect. 3.2, the number of contributions to the sum on the r.h.s. of eq. (4.28) scales approximately as the square of the number of light partons in the final state. We shall show in the following that within the FKS formalism it is trivial to reduce this scaling property from a quadratic form to a linear function at most, and actually to a constant if one only increases the number of final-state gluons. We shall also discuss in this section the integration of the n-body contribution, and the implementation of multi-channel integration.

Reduction of the number of independent contributions
Since the number of contributions to eq. (4.28) coincides with the number of elements in P FKS , we start by observing that the definition given in eq. (3.37) is sufficient to include all singularities to be eventually subtracted, but is redundant. As an example, consider the case of I i = q, I j =q, with q any massless quark. According to eq. (3.37), both (i, j) and (j, i) belong to P FKS . On the other hand, the pairs (i, j) and (j, i) are responsible for the same divergence of the matrix element, the collinear one at k i k j . This divergence is subtracted twice in eq. (4.28), in the terms dσ , but is not double counted thanks to the presence of h ij (z ij ) and of h ji (z ji ) in the relevant S functions (see eq. (5.14)). Furthermore, as discussed in sect. 4.2, when the FKS parton is massive, all of the contributions dσ (n+1) ij are finite without subtraction; it is therefore not particularly advantageous to have a massive FKS parton. Several other examples can be given, the upshot of which is that it is more convenient to define P FKS as follows, rather than as in eq. (3.37): (6.1) Equation (6.1) differs from eq. (3.37) in the upper limit for the range of i, which implies that FKS partons are now restricted to be massless (this, in turn, has the consequence that the the condition k j → 0 need not be used any longer); and because of the presence of non-redundancy conditions, which we specify as follows: 2) At this point, the sum on the r.h.s. of eq. (4.28) still scales quadratically with the number of light partons; however, that sum obviously contains several identical contributions. Suppose for example that one has m final-state gluons i 1 . . . i m . There will be m(m − 1) gluon pairs in P FKS , which will all contribute to dσ (n+1) . However, it is clear that dσ iγ i δ , because of the symmetries of the matrix elements, S functions, phase space, and (IR-safe) observables. Likewise, if j is not a gluon, and is such that (i 1 , j) ∈ P FKS , then (i α , j) ∈ P FKS , for 2 ≤ α ≤ m, and dσ In general, the determination of which elements of P FKS give identical contributions depends on the process considered, and on the underlying theory (e.g., QED only or full EW interactions), and can be implemented with relative ease in a computer code. It is clear, however, that if several identical particles are present in the final state, one and only one of them (which one in particular, is irrelevant) will be a member of the FKS pairs. We shall denote by the subset of P FKS whose elements give non-identical contributions to the sum in eq. (4.28). We shall then compute with ς (n+1) ij an integer symmetry factor, equal to the number of identical contributions (for a given (i, j)) to the sum in eq. (4.28). The number of elements in P FKS (r) and P FKS (r) will be denoted by N FKS (r) and N FKS (r) respectively; the value of N FKS for several production processes will be reported in sect. 7.

n-body matrix elements
The integration of the n-body contribution, eq. (4.3), is straightforward, and is usually carried out independently of that of the (n + 1)-body contribution. In this section, we give the formulae that allow one to integrate these two contributions at the same time. For instance, this technique is the default in MC@NLO. We start by observing that if r α ∈ R n =⇒ ∃r β ∈ R n+1 such that r α = r i \ β with I i = g . (6.8) This states the fact that, given a set of particles that corresponds to a physical 2 → n process, by adding one gluon in the final state one gets a physically meaningful 2 → n + 1 process. It therefore follows that Note that, for an arbitrary 2 → n + 1 process r β , the corresponding 2 → n process r i \ β does not necessarily exist (e.g. if I i = q). In such a case, we have dσ (n) (r i \ β ) = 0. On the other hand, if one fixes i, if r i \ β exists then it is unique, and hence eq. (6.9) holds.
Let us now introduce the shorthand notation 10) The precise form of B can be read from eqs. (4.4), (4.5), (4.12), and (4.14), but is not relevant in what follows. Fixing i and assuming r i \ is a physical 2 → n process, we have, using eq. (4.31), Using eq. (4.18) we can also write This equation holds for any given i with I i = g. It is convenient to give a prescription for fixing i given the process. In order to do that, we observe that there is only one gluon among the first elements of the pairs that belong to P FKS , and this gluon is therefore a natural candidate for fixing i. We rewrite eq. (6.12) as follows: where on the l.h.s. of eq. (6.13) the use of a 2 → n + 1 process as an argument for dσ (n) is justified by the fact that there is a one-to-one correspondence between r and r i \ if one considers P FKS .
After these manipulations, the physical cross section of eq. (4.1) is written as follows: with dσ (n+1) given in eq. (4.40), and By inspection of eqs. (6.14) and (4.38), we see that the second term on the r.h.s. of eq. (6.16) has the very same structure as the soft counterevent of the first term. When cast in this form, it thus becomes obvious that the (n + 1)-and n-body contributions to an NLO cross section can be integrated together.
It is also easy to show that the degenerate (n + 1)-body contribution dσ (n+1) can be integrated together with the other two contributions. This is the default strategy in MC@NLO. We refrain from reporting the relevant derivation here, and postpone it to a forthcoming paper, where the implementation of the FKS subtraction will be extended to the case of hadronic collisions.
We conclude this section by observing that the logic behind eqs. (6.15) and (6.16) is that of getting the underlying Born process given a real-emission process. This structure matches quite smoothly that of MadGraph, which is the reason why we have chosen it in the first version of MadFKS. On the other hand, one may also consider the inverse logic, namely that of getting the real-emission process(es) given a Born process, which is especially convenient in the context of the matching of NLO computations with parton showers. We shall describe this option in app. E.

Multi-channel integration
In the previous sections we have shown how to define a number N FKS of finite and independent short-distance contributions, whose sum is the observable cross section. For large final-state multiplicities, the integration of each of these contributions may prove difficult from the numerical point of view.
Let us start by summarizing the strategy adopted in a tree-level computation by MadGraph/MadEvent [33]. If the tree-level matrix elements are those relevant to a 2 → n + 1 process, one can use eq. (4.1), set dσ (n+1) = 0 and dσ (n) = 0 there, and define where J n (B) L +1 denotes that one reconstructs exactly n (B) L + 1 jets in the final state, which is sufficient to cut off all phase-space singularities of M (n+1,0) . Let us now denote by d α (r), α = 1, . . . N d (r) (6.18) the subset of Feynman diagrams that contribute to M (n+1,0) , and that cannot be obtained from each other with a permutation of four-momenta. In other words, MadGraph then splits the integral of dσ (n+1) in eq. (6.17) into N d (r) independent contributions 12 , called integration channels, and defined as follows: where and apply the subtraction procedure to this quantity. Unfortunately, this is not going to work. In the soft limit E i → 0, all diagrams that would induce a divergence in the square of the amplitude in eq. (6.19) vanish if squared individually, as in eq. (6.22), owing to the presence of a self-eikonal of a massless particle (an exception is therefore that of diagrams in which a gluon is emitted by a final-state massive particle). It follows that the cross sections dσ (n+1) α that are equal to zero are actually those that one would like to have the largest contributions to the sum in eq. (6.20).
From the physics viewpoint, this is easy to understand. The squares of individual Feynman diagrams work well in tree-level computations because these computations are basically classical physics. On the other hand, soft singularities are inherently a quantum effect, since they are due to the interference of diagrams. Thus, multichannel sampling based on squares of diagrams cannot possibly describe all complications due to interference.
We can bypass this problem in the following way. We first introduce a monotonically decreasing, smooth function f (x) with the following properties: We then modify eq. (4.24) to read The factor in square brackets in the second term on the r.h.s. of eq. (6.26) damps the only soft and collinear singularities that survive the damping of S ij . That second term is therefore finite over the whole (n + 1)-body phase space, and is thus effectively a tree-level computation that can be dealt with in the same way as in standard MadGraph. The effect of the function f in the first term on the r.h.s. of eq. (6.26) is that of suppressing hard and large-angle emissions of the (i, j) system. We can therefore assume that the typical kinematic configuration will always be well approximated by the degenerate configuration obtained by letting particles i and j become collinear, or when particle i becomes soft. As we have seen in sect. 5.1, these two limits result in the same four-momentum configurations, and we can choose one or the other depending on which of the two will result in a physical 2 → n process as far as particle identities are concerned. In practice, if I i = g we shall use the soft limit, while if I i = g the collinear limit will be used. To avoid unnecessarily complicating the notation, we shall denote both cases by r i \ in this discussion. Starting from eq. (6.26), we get (6.28) We stress again that, in spite of the presence of J n (B) L rather than that of J n (B) L +1 , the cross section in eq. (6.28) is finite, and needs no subtraction. We have also defined where dσ (n+1) ij has been given in eq. (4.29), and we applied in eq. (6.27) the same symmetry factor ς (n+1) ij as in eq. (6.7) for obvious reasons. As discussed above, the kinematic configuration used to evaluate the prefactor |d α | 2 /D in eq. (6.29) is that of either the soft or collinear counterevents, which coincide. Therefore, this prefactor will have the same value in the case of the event and of all the counterevents. On the other hand, f (2k i ·k j /s) will in general be different when computed with event and counterevent kinematics, but the smoothness of this function guarantees that the cancellation of phase-space singularities will still take place.
Note that the formulae above can be applied without modifications even if one sets f (x) ≡ 1. This means that n-body Feynman diagrams are used for multichannel sampling in the whole (n + 1)-body phase space. This may imply a loss of efficiency when hard-emission regions are integrated over. However, as we shall see in sect. 7, this is not the case for any of the processes we have considered in this paper.
From eq. (6.27), we can read the total number of integration channels of our implementation of an NLO cross section, for a given process r: Note that the second term on the r.h.s. of this equation is roughly proportional to N FKS (it would be exactly so if all N d (r i \ ) were equal). So while N FKS is directly related to the performance of the FKS formalism in keeping the number of independent subtraction terms as small as possible, the ratio N ch /N FKS can be used to measure the increase of the number of integration channels due to the inherent complexity of the underlying Feynman diagram structure. Note that, by setting f (x) ≡ 1, eq. (6.30) becomes N ch (r) =

Results
In this section we present MadFKS results for several production processes in e + e − or µ + µ − collisions. These correspond to the full NLO cross section, eqs. (4.1) or (6.15), where we set V (n,1) in eq. (4.14). This is equivalent to assuming that the one-loop contribution is a pure-pole term in the CDR scheme (see app. B), which implies that the results given in this section are non physical. Indeed, our aim is not that of carrying out a phenomenological study, but to document the convergence properties of the FKS subtraction in MadFKS (i.e., the statistical errors obtained with a given number of integration points), and to make sure that our implementation is correct. We point out that the optimizations implemented so far are restricted to the FKS formalism; from the point of view of MadGraph, the current version of MadFKS must be considered as a benchmark upon which we shall build a fully efficient program. In particular, there are two main lines of development for the code. Firstly, the results we present here have been obtained by running all integration channels with the same number of points. This procedure will be improved following the same strategy as is used in MadGraph, namely that of performing preliminary low-statistics runs to determine which channels give significant contributions, in order to further the runs with larger statistics only for those. Secondly, in the current version of MadFKS all sums over helicities are performed exactly. A much more efficient procedure is that of performing such a sum using Monte Carlo techniques, as is typically done in the context of tree-level computations. This is also possible at the NLO, provided one is able to write the subtraction formalism for any given helicity configuration. This is straightforward in the context of the FKS subtraction method, since indeed all necessary ingredients are already available in the literature [12,34]; we shall further comment on this point, and give all the relevant details, in app. F. We have refrained from implementing these improvements (and several other minor ones) in the first version of MadFKS, since this is the best way to gauge in a clear manner the capabilities of the code; for example, by treating all channels on the same footing, one is able to check that they are all computed correctly, regardless of whether some are numerically relevant or not.
As explained in sect. 6, to obtain a cross section we find it convenient to sum over (n + 1)-body processes, according to eq. (6.15). In this section, therefore, when considering a given partonic process, it will always be understood as an (n + 1)body one. The corresponding n-body contributions will be added according to the procedure described in sect. 6.2.
We have set up two different types of numerical tests to check the correctness of the implementation. The first type of tests make sure that the local behaviour of our main ingredients in the computation of the real-emission contribution, namely the matrix elements and the S functions, behave as we expect them to do in the collinear and soft limits. For all real-emission processes r, and for all pairs (i, j) ∈ P FKS , we check numerically that eqs. (5.17), (5.18), and (4.16)-(4.20) are satisfied. This is done automatically, and if one of these tests fails, the programme proceeds no further.
The second type of tests rely on the fact that several free parameters are introduced in the FKS short-distance cross sections; the final results must be independent of these parameters, but partial results must not. To give one example, let us consider the case of ξ cut , introduced in sect. 4. This parameter enters analytically the n-body cross sections (see e.g. eq. (4.6) and the results for the integrated eikonals in app. A), and numerically the (n + 1)-body cross section through the definition of the soft counterterm, according to eq. (4.33). Therefore, the quantities dσ (n+1) and dσ (n) will separately depend on ξ cut , but their sum will not. The free parameters whose variations we have considered here are ξ cut , δ O , a S , b S . We have also investigated if the final results are independent of whether one uses the factor 2E i / √ s in the definition of d ij (this option is denoted by useenergy=.true.), or one uses 1 instead (this option is denoted by useenergy=.false.); see sect. 5.2 for more details. Each of these parameter choices corresponds to a set of (parallel) runs. Therefore, to be definite, we have used the (n + 1)-body process as a benchmark; this is a contribution to the four-jet cross section at the NLO. We have run at √ s = 100 GeV, and set the factorization, renormalization and Ellis-Sexton scales equal to the Z mass, µ 2 F = µ 2 R = Q 2 = m 2 Z . This process is indeed simple enough to run in a short amount of time, but has the required complexity to test the implementation in all its details. The condition in eq. (3.36), i.e. the hard cuts, has been imposed with the KTCLUS routine that implements the jet-finding algorithm of ref. [35], with Y cut = (10 GeV) 2 . The jet four-momenta are defined as the sum of the four-momenta of the partons in the jets; the lowest jet energy is therefore about 10 GeV. The cross sections and integration errors (as returned by Vegas [36]) for the process in eq. (7.2) are given in table 1. There are N FKS = 3 independent FKS pairs contributing to this process (these can be chosen to be (6, 3), (6,4), and (6, 7) according to the labels of eq. (7.2)), and five integration channels for each of them. We used 50000 integration points and 10 iterations per channel; only about 15% of these points result in kinematic configurations that pass the jetfinding conditions; at the level of phase space, the points are generated flat, i.e. no information on the hard cuts is available when one generates the kinematics. The results of table 1 are clearly independent 13 of the free parameters, typically within one standard deviation. As an even tighter consistency check, we have picked the largest and the smallest values among those in the table, and re-evaluated them with higher statistics. For a six-fold increase of the number of iterations, keeping the points per iteration to the same value as before, the corresponding cross sections change as follow: 3.5633 ± 0.0154 → 3.6086 ± 0.0051 , 3.6350 ± 0.0151 → 3.6007 ± 0.0053 . It is interesting to observe that, apart from the convergence of the mean values towards a common value, the errors in eq. (7.3) scale approximately as we would expect them to do if they would follow the Gaussian law typical of integrals of ordinary functions (while a subtracted cross section is actually a distribution). This gives us confidence on the fact that Vegas estimates correctly the integration errors.
We have also checked that the cross section is independent of the choice of the function h (see eqs. (6.5) and (4.23)); since the numerical effects are even smaller than those reported in table 1, we refrain from presenting the corresponding results here. We have performed the runs with It is remarkable that, with this choice, the contribution of dσ (n+1) 0,β (see eq. (6.28)) is smaller than the statistical errors reported in table 1, and is therefore not included in the results of the table.
As a final comment on the process in eq. (7.2), we point out a fairly pleasant feature that could have been anticipated. Namely, the results presented in table 1 have integration errors that are by far and large independent of the choices of the free parameters, thus including the case of "maximal" subtraction, ξ cut = ξ max and δ O = 2, when for each event one always subtracts the three counterevents. This is quite an useful property, since it implies that it will not be necessary to scan the space of the free parameters in order to find those values that maximize the numerical stability, all choices being almost equivalent. At the core of this property stands the fact that, for any given FKS pair (i.e., a finite and independent contribution), when scanning the (n + 1)-body phase space one finds at most one soft and one collinear singularities which, as discussed in sect. 5.1, are associated with a single kinematic configuration. Therefore, the behaviour of the subtracted cross section will not change much by varying ξ cut and δ. It is instructive to compare this with the case of dipole subtraction. There, by changing the parameter that is typically called α (the analogue of our ξ cut and δ), one can change drastically the number of subtraction terms included in the computation, and the tuning of α may become necessary. The capability of the FKS subtraction method to give numerically stable results also for large values of ξ cut and δ is especially promising in view of the fact that such values are usually the optimal choices when matching NLO computations with parton showers in MC@NLO or POWHEG.
The independence of the results of the parameters ξ cut and δ O allows one to check the correctness of the subtraction terms for the real-emission matrix elements, and of the logarithmic terms in the n-body contributions. In order to also check the non-logarithmic terms that enter the n-body contributions, we have computed the physical NLO cross sections for the processes e + e − → Z → 2 jets and H → 2 jets, where H is a SM Higgs with m H = 120 GeV, using the readily available relevant matrix elements for the virtual corrections. We have found complete agreement with the known results. We point out that, by doing so, we have tested all analytical formulae presented in this paper, except for eqs. e + e − → Z → uūggg (0.3601 ± 0.0014 (0.38%))×10 1 3 We have considered a variety of other processes, in order to prove the flexibility of the code, and to check the convergence of the numerical integration; the processes and the corresponding results are reported in table 2. We stress that, in order to obtain the results displayed in that table, we just had to write the relevant process in an input card for MadFKS; this is the same procedure as that done in MadGraph. All cross sections appearing in table 2 have been computed by setting ξ cut = 0.1, δ O = 0.6, a S = 1.5, b S = 1.5, and useenergy=.true.. We have run 10 Vegas iterations; for each of them, phase-space points have been generated in order to obtain about 10000 kinematic configurations per iteration that pass the jet-finding cuts, which corresponds to the statistics used to generate the process in eq. (7.2). For final states with only massless partons, we set √ s = 100 GeV; we have used the notation j to indicate the sum over partons (we have considered four flavours here); thus, e.g. e + e − → jjjjjj is, up to the missing finite one-loop contribution, the physical five-jet cross section. For top-quark production processes, we set √ s = 500 GeV.
The b quark is taken to be massive (m b = 4.7 GeV) in processes with ttbb + X final states, and it does not enter the jet-finding algorithm. Finally, for stop and gluino production we have used mt 1 = 400 GeV, mg = 400 GeV, and √ s = 1 TeV; we assume an effective Zgg axial vertex, and set the corresponding coupling equal to one. We have also studied the production of fully gluonic final states arising from the decay of an (off-shell) SM Higgs with m H = 120 GeV, that contribute to three-, four-, and five-jet observables. In order to keep this process on the same footing as the others that feature only massless final-state particles, we have set √ s = 100 GeV, but used a µ + µ − initial state in view of the smallness of the electron Yukawa coupling. As one can see from table 2, fully-gluonic Higgs decays have a single non-trivial contribution from FKS pairs. In the cases of processes involving massive particles we also tested the option of replacing β k and/or β l with 1 in the expressions of the relevant d kl (see sect. 5.2). We did not observe any significant differences with respect to the results obtained with the default choice. As one can see from the table, the numerical errors are fairly small even with the limited statistics we used. The largest of them (but still rather modest) are associated with the processes with a tt pair in the final state, plus either four gluons or a bb pair and extra gluons. It is clear that these processes feature several mass scales quite different from each other -the c.m. energy, the top mass, the b mass, and the minimum jet energy (∼ 10 GeV) -and therefore one expects the coefficients of the perturbative series to be plagued by several large logarithms, that inherently increase the complexity of the numerical calculation. We point out that, since we treat the b quark as massive, the bb pair is not an FKS pair, according to eq. (6.1). On the other hand, given the smallness of the b mass, it may be beneficial to treat such a pair as an FKS pair, which would also possibly imply the definition of a "quasicollinear" counterterm, whose analytical form has to tend to that in eq. (5.18) in the limit m b → 0. We leave the implementation of this option to a future work.
It is important to note that the growth of N FKS with the final-state multiplicity is always rather modest. From table 2, one can actually see that N FKS is a constant if the number of gluons is increased (N FKS for ttbbgg is one unit larger than in the case of ttbbg, since only in the former process one can form an FKS pair with two gluons). The increase of N FKS in the e + e − → (n + 1)j processes is mainly due to the corresponding increase in partonic subprocesses: in the case of four, five, and six particles in the final state of the (n + 1)-body processes, the number of contributing subprocesses is equal to 7, 7, and 17 respectively. Among individual contributions to the five-jet cross section, those with the largest N FKS are the four-quark, differentflavour processes such as e + e − → uūddgg, which have seven possible FKS pairs (since one distinguishes between quarks and antiquarks when pairing them with gluons, in The histograms for the latter two processes have been rescaled (by a factor of 9 and 250 respectively) in order to fit into the layout. We present thrust, C parameter, and the energy and polar angle of the leading jet.
order not to neglect possible charge asymmetries).
As discussed in sect. 6.2, we integrate the n-body matrix elements at the same time as the (n + 1)-body ones. On an event-by-event basis, we can therefore obtain both the NLO and the LO contributions. We have checked that the latter is, for all processes, fully consistent with the one predicted by standard MadGraph. If one switches off in MadGraph the optimizations relevant to the separate treatment of different integration channels, our LO computation has the same statistical accuracy as that in standard MadGraph. More importantly, if we only integrate the Born contributions to the processes listed in table 2 with the same number of points as that used for the NLO contributions (distributed equally among the possibly smaller number of integration channels), the resulting integration uncertainties are a relative factor 1.9 to 4.5 smaller than those relevant to the NLO results presented here. Exceptions are found for the (n + 1)-body processes ttbbg and ttgggg, whose LO contributions have integration accuracies better than the NLO ones by a factor 7 and 9 respectively. Overall, these figures give us another indication of the fact that the subtraction of singularities is achieved in very satisfactory manner from the numerical point of view.
We now turn to presenting differential distributions for some of the processes listed in table 2. As clarified in sect. 5 (see in particular sect. 5.1), for a given process and choice of (i, j) ∈ P FKS , to each random number there correspond two kinematic configurations (for the event and for the counterevents plus n-body contributions), with associated weights. These configurations give the complete information on all the four-momenta of the (resolved) final-state particles, and thus one can plot as many (infrared-safe) observables as one wants in the course of a run. The plots we present here are obtained with the same statistics (in fact, in the same runs) as that The histograms for the latter two processes have been rescaled (by a factor of 6 and 50 respectively) in order to fit into the layout. We present the energy and polar angle of the top quark and of the leading jet.
used for the total cross section results given in table 2. We point out that we do not use any smoothing procedure for the spectra we show: we limit ourselves to filling the histograms putting the weights in the bins determined by the kinematic configurations given by the code. This implies, in particular, that large weights for the event and counterevents may fall into different bins. We choose a relatively fine binning, in order to expose in a clear manner whether this mis-binning (which is unavoidable in any subtraction-based computation) is a severe problem or not.
In fig. 1 we show the distributions in the thrust, C-parameter, and energy and polar angle of the hardest jet, for the first three processes listed in table 2. In fig. 2, the same distributions are shown for three-, four-, and five-jet final states, and including the Z/γ interference effects -the plots in fig. 2 would therefore correspond to physical jet cross sections, had we included the proper one-loop contributions.
Finally, in fig. 3 we show the energies and polar angles of the top quark and of the hardest jet relevant to the e + e − → tt + (n + 1)g processes. All spectra are fairly smooth, except those relevant to the polar angles, where fluctuations are marginal for light-jet-only production processes, and more evident in the case of tt processes.

Conclusions
In this paper, we have presented a complete automation of the computation of any cross section at the next-to-leading order in QCD, with the exception of an infraredand ultraviolet-finite term of one-loop origin. This is achieved by embedding into the MadGraph framework the universal subtraction formalism of Frixione, Kunszt, and Signer (FKS); the resulting computer code has been named MadFKS. The only intervention needed to run the code is the definition of the inputs, among which is the production process one wants to study. The process is any physical process that results from a theory implemented in MadGraph, thus including the Standard Model and any user-defined Beyond the Standard Model theory. The current version of the implementation allows one to compute cross sections with colourless particles in the initial state, but there are no difficulties in extending it to other types of colliding particles; we defer this to future work. It is also worth mentioning that the formalism presented here requires only trivial modifications to be extended to the subtraction of QED infrared singularities. With minimal effort, one could therefore fully automate the FKS subtraction for the complete electro-weak corrections; we point out, however, that in order to actually perform such computations it will probably be necessary to implement a complex-mass scheme into MadGraph.
The key features of our implementation are a direct consequence of the FKS subtraction formalism. Through a dynamic partition of the phase space, one effectively defines partonic processes which are either non singular, or have at most one soft and one collinear singularity, thus rendering the subtraction of these singularities a fairly easy, and numerically efficient, task. These partonic processes are separately finite and fully independent of each other, and therefore the method is naturally suited to parallel computations. For each phase-space point, there are only two independent kinematic configurations; one corresponds to a fully-resolved configuration (all partons are hard and well separated from each other), and the other one is the associated partly unresolved configuration (where either one parton is soft, or two partons are collinear, or both), whose weight is the sum of all subtraction terms, of all finite remainders of the subtraction procedure, and of the Born; this structure is exact, i.e. does not imply any approximation. Being able to limit the independent kinematic configurations to two is quite beneficial for the numerical stability of the results, and in particular for the smoothness of the predictions for differential spectra. The latter point is easy to understand, since in FKS one renders it small the probability that weights large in absolute value and opposite in sign will not fall in the same histogram bin (thus, one reduces the mis-binning problem that affects all computations based on a subtraction method).
The FKS subtraction method allows one to introduce several free parameters in the calculation; partial results (e.g., the (n + 1)-body and the n-body contributions) depend on the values chosen for these parameters, but the physical cross section does not. We have exploited this property to check the correctness of our implementation. In addition, we have shown that the integration errors are largely independent of the choices of the free parameters, which is a direct consequence of having only two independent kinematic configurations for each phase space point. This is an important fact, since it allows one to choose the free parameters without having to worry about the numerical stability of the results.
We have presented predictions for several production processes in e + e − or µ + µ − collisions. The numerical integrations converge quickly, requiring relatively modest statistics to achieve accurate results. We have run in parallel, with each of the parallel jobs corresponding to one integration channel, and using the same statistics. This is certainly not the best way to perform a multi-channel computation. The most obvious improvement, which is straightforward to automate, is that of starting with a low statistics run for all channels, checking which contributions are largest in absolute value or/and have the largest fractional errors, and further the run with larger statistics only for those channels. We plan to release a public web interface to MadFKS as soon as the case of hadronic collisions will also be included, and some of the optimizations discussed in this paper will be implemented.
The current version of MadGraph, upon which MadFKS is constructed, is based on the evaluation of Feynman diagrams. We point out, however, that MadFKS can be rendered fully independent of Feynman diagram information, thus implying the possibility of using recursively-generated amplitudes. In fact, the only part in the present implementation which is strictly dependent on Feynman graphs is that relevant to multi-channel sampling. A different criterion for sampling, based e.g. on topological information or on the computation of the amplitudes, can be easily incorporated in MadFKS (where multi-channel sampling is the very last step of the computation).
If combined with one-loop computations, MadFKS will allow one to predict any NLO cross section, limited only by CPU time. The present bound is actually set by the capability of MadGraph to handle the large number of Feynman diagrams and colour configurations for very large final-state multiplicities; tree-level amplitude computations based on recursion relations will push this bound to larger multiplicities. The computer programme can now use any external routine that returns the finite one-loop contribution as defined in this paper, feeding it with four-momenta and information on the tree-level structure. It is clearly desirable to fully incorporate one-loop computations into the programme, which appears to be feasible given the recent progress made in the evaluation of virtual graphs. Either way, the implementation presented here is the first step towards a full automation of NLO computations matched with partons showers, as done in MC@NLO or POWHEG.

Acknowledgments
SF would like to thank Zoltan Kunszt and Werner Vogelsang for discussions during the course of this work, and for a long-time collaboration on related subjects; and Eric Laenen, Patrick Motylinski, and Bryan Webber for work relevant to FKS subtraction with massive particles. FM is grateful to CERN for hospitality during the course of this work. RF and FM are partially funded by Technical and Cultural Affairs through the Interuniversity Attraction Pole P6/11. This work was supported in part by the National Science Foundation under award PHY-0757889.

A. Eikonal integrals
We define the integrals of the eikonal factors as follows: In eq. (A.1),Ê kl collects by definition all divergent terms (times an ǫ-dependent pre-factor, to be given explicitly below), while E kl is finite. The normalization of eq. (A.1) is conventional, and is equal to that used in ref. [12], times a factor 8π 2 (this factor has been compensated by a different normalization of the colour-linked Born's M (n,0) kl w.r.t. that of ref. [12]). The energy of parton i and the measure (in 3 − 2ǫ dimensions) over its angular variables dΩ i are defined in the c.m. frame of the colliding partons (see eq. (4.25)). In ref. [12] only the massless case m k = 0 and m l = 0 had been considered. The cases m k = 0, m l = 0, and m k = 0 with k = l, have been dealt with in ref. [31], but the analytical results for eq. (A.1) have not been published there. In this paper, we have computed the only missing ingredient, i.e. m k = 0 and m l = 0. We summarize here all the relevant results.
The definition of β l has been given in eq. (5.12).
where we have introduced the function (A.14) .  2) of ref. [12] to the case of arbitrary particle masses, and is consistent with the results given in ref. [37]. We stress again that the poles in eq. (B.2) have infrared origin, since the ultraviolet divergences are assumed to have been eliminated through renormalization. The Ellis-Sexton scale Q, originally introduced in ref. [32], is any scale that may be used in one-loop computations to express the arguments of all the logarithms appearing there as s ij /Q 2 rather than as s ij /s kl , where s ij and s kl are two invariants constructed with the four-momenta of the particles that enter the hard process (see e.g. sect. 6 of ref. [38]). Clearly, it is not mandatory to introduce the Ellis-Sexton scale in the computations of the virtual amplitudes, and fairly often one just uses the factorization or renormalization scale instead. On the other hand, if Q is kept distinct from both µ F and µ R , one sees that the NLO cross section is exactly independent of its choice; in other words, variations of Q do not lead to any estimate of missing higher orders in perturbation theory, but allow one to check the internal consistency of the implementation. Although indeed most of the available one-loop results do not use the Ellis-Sexton scale, it is always possible to introduce it; we show how this can be done in app. C.
Once the virtual amplitude is computed, one can interfere it with the Born amplitude and obtain M (n,1) . Thus, eqs. C. How to set µ F = µ R The formulae presented in this paper assume µ = µ F = µ R , and that another arbitrary mass scale, the Ellis-Sexton scale Q, is introduced in the computation of the one-loop corrections. In this appendix we explain how to relax the former condition, and how to insert the dependence on the Ellis-Sexton scale in an oneloop result that was computed without introducing it. We perform these tasks by imposing renormalization group invariance w.r.t. µ F and µ R . We start by giving our convention for the scales that appear in the various formulae of this paper: If µ F = µ R , all the formulae for the short-distance cross sections given in this paper must be computed with µ = µ F , except for the argument of α S , which must be set equal to µ R .
The statement above implies in particular that the PDFs are computed at µ = µ F , as customary. The convention given above cannot possibly allow us to recover terms proportional to logarithms of ratios of scales, which we now proceed to determine. In order to do so, we rewrite eq. (4.1) or eq. (6.15) symbolically as follows with C the unknown quantity that we need to determine. In order to proceed, we have to specify the power of α S that enters our cross section formulae; we denote by b an integer such that Note that in general we may have b = 0 (e.g. in the case of a purely-EW process at the Born level). We can now impose the invariance of the physical cross section w.r.t. to the renormalization scale where as usual this equation holds up to terms of O α b+2 S . Using the explicit forms of the short-distance cross sections that appear on the r.h.s. of eq. (C.1), and the convention stated above, we find that eq. (C.3) is equivalent to we obtain where dσ (B,n) is the Born cross section, given in eq. (4.4). This would be the end of the story if all short-distance cross sections could be computed according to the convention given above. This may not be the case for the one-loop contribution. We can however always manage to rewrite the one-loop contribution using our convention, regardless of its original form. To do this, we impose the condition that the physical cross section be independent of the factorization scale: which is the analogue of eq. (C.3). After some algebra (where one also uses the evolution equations for the PDFs), eq. (C.7) is found to be equivalent to the condition with C given in eq. (C.6). In order to proceed, let us thus consider eq. (4.14), and denote by t k all possible quantities with mass-squared dimensions that can be constructed with four-momenta. Following our convention, we must write where use has been made of the fact that, in an one-loop computation, there will be an explicit linear dependence on the logarithm of the renormalization scale; according to our convention, in the argument of this logarithm µ R must be replaced by µ F , hence eq. (C.10). On the other hand, the computation will in general contain terms of the kind log p t k /t l , with p = 1, 2. The procedure of Ellis and Sexton implies that these must be replaced by (log t k /Q 2 − log t l /Q 2 ) p , and b V in eq. (C.10) may thus include up to double logarithms of Q 2 ; this is however unimportant in what follows. Equation (C.8) can now be seen as a relation to solve for a V ; we obtain V (n,1) Note that the r.h.s. of this equation is indeed independent of µ R , since M (n,0) is proportional to α b S (µ 2 R ). Equations (C.9) and (C.11) can be used to insert the renormalization, factorization, and Ellis-Sexton scale dependences into a one-loop result originally computed with a single scale M. In other words, if we are given the quantity v(M 2 ) such that V (n,1) then we can construct the finite contribution of one-loop origin which complies with our convention as follows: .
(C. 13) We stress that the quantity v(Q 2 ) is proportional to α b S (Q 2 ), and therefore the whole r.h.s. of eq. (C.13) is proportional to α b S (µ 2 R ). In summary, we shall compute the cross section with the most general assignment of scales µ F = µ R = Q using (C.14) where the scales in the formulae for the short-distance cross sections that appear on the r.h.s. of this equation have to be chosen according to the convention given at the beginning of this appendix. When following this convention, the finite part of the one-loop contribution, eq. (4.14), can be derived using eq. (C.13) from a result v(M 2 ) obtained with a single scale M.

D. Azimuthal terms in collinear limits
When computing the collinear limit of the damped real matrix elements, we obtain eq. (5.18) (or the equivalent form eq. (5.19)). It is customary to write such a limit by only keeping the first term on the r.h.s. of eq. (5.18) but, as discussed in sect. 5.3, the second term is equally important if one wants to construct local counterterms. Appendix B of ref. [12] explains in detail how to compute the reduced matrix element M (n,0) , and the kernels Q ab ⋆ (z). These kernels are universal in the same sense as the Altarelli-Parisi ones. The notation suggests that already at the leading order (which is not the case for the Altarelli-Parisi kernels) their forms for timelike branchings are different from those for spacelike branchings. From ref. [12] we obtain M (n,0) r; k (C) = 1 2s The notation A (n,0) ± r j⊕i,i \ (D. 2) that appears in eq. (D.1) understands that parton j ⊕ i has helicity equal to ±. M (n,0) is therefore an interference between the plus and minus helicity states of parton j ⊕ i; the remaining helicities (i.e. excluding that of parton j ⊕ i) are summed over. The first term inside the curly brackets on the r.h.s. of eq. (D.1) is the ratio of two spinor products. We follow here the notation and conventions of Mangano and Parke, ref. [39]. We note that this term is a pure phase, and is therefore numerically well defined also in the soft limit (E i → 0); we have worked out its form analytically for the various phase-space parametrizations we have considered.
For final-state (i.e. timelike) branchings, the Q kernels read as follow [12]: The corresponding results for initial-state (i.e. spacelike) branchings are: The ⋆ symbol allows one to remember easily which parton is off-shell. We note that for branchings which involve an off-shell quark, the azimuthal term in the collinear limit vanishes identically. This is due to the conservation of the helicity along a light-quark line; see ref. [12] for details. We point out that in eq. (D.1) one must use the same conventions for the spinors in the computation of the amplitudes as in that of the pure-phase factor. As mentioned above, eq. (D.1) has been worked out using the conventions of ref. [39]. On the other hand, in the actual numerical implementation in MadGraph the amplitudes are computed following the conventions of HELAS [40]. Since eqs. (D.3)-(D.6) imply that we need to consider only the cases in which parton j ⊕ i is a gluon, we can write (D.11) Using the explicit representations of ref. [39] and of HELAS, it is a matter of algebra to arrive at ε µ (+, k; MP) = exp (−i ϕ k ) ε µ (+, k; HELAS) + αk , (D.12) where k is the gluon (light-like) four-momentum, which in the coordinate system chosen by HELAS has azimuthal angle equal to ϕ k , and α, β are two constants, which are irrelevant in what follows, since the corresponding contributions to eq. (D.11) vanish because of gauge invariance. Using these equations we obtain We conclude this section by reporting the 4 − 2ǫ dimensional forms of the unpolarized Altarelli-Parisi kernels for z < 1 [41]: According to eq. (4.44), the coefficients of the O(ǫ 0 ) and O(ǫ 1 ) terms in these equations are the kernels P (0) ab and P (1) ab respectively, which are used in several equations of this paper.

E. Possible variants in the implementation
As we have discussed in the text, there are a large number of arbitrary parameters which enter our implementation of the FKS formalism into MadGraph, and any physical observable is independent of them. Each parameter choice can be considered as a different way to implementing the subtraction procedure in a computer code. In this appendix we report two variants of the implementation which are not parametric in nature, that we may want to consider in the future, and which aim at further reducing the number of independent contributions to the physical cross section. We also present a technique alternative to that discussed in sect. 6.2 to integrate together the n-and (n + 1)-body contributions, which is best suited for NLO computations matched with parton showers.
We start by considering the contributions to the (n + 1)-body cross sections due to (i, 1) and (i, 2) (for a given i, and assuming that both of these two pairs belong to P FKS ). Using the properties of the S functions, it is a matter of simple algebra to show that Here we have set which follows from the fact that we work in the c.m. frame of the incoming partons, and thus we can use eq. (4.25). Eq. (E.1) is the form used in refs. [12,13], and has the virtue of allowing one to subtract the two initial-state collinear singularities in one single integration. On the other hand, eq. (4.29) is closer to what MadGraph does when it integrates tree-level matrix elements. We postpone the comparison between these two possibilities to a forthcoming publication. We point out that in refs. [12,13] the notation was used S ij ≡ S ij denoted the S functions relevant to both i and j in the final state.
We now turn to considering the contributions to the (n + 1)-body cross sections due to all the pairs formed by a given FKS parton i and its massive sisters. Since m j = 0 for such sisters, the matrix elements are not singular in the collinear limits. This implies that the distribution in y ij defined in eq. (4.34) is just a regular function in eq. (4.29), and therefore there is a cancellation of the factors (1 − y ij ) between the numerator and the denominator. We thus immediately obtain On the r.h.s. of eq. (E.3) we have changed variables y ij → y i for all j. The variable y i can be equal to y ij for a given j, or can be defined as in eq. (E.2). The specific choice of y i should not matter much in terms of numerical performance, given that the collinear regions y ij = 1 are not associated with any peaks in the matrix elements (this may not be the case for particles with small masses, such as b quarks). Note also that multi-channel integration can be performed in eq. (E.3) following exactly the procedure that leads to eq. (6.29).
In the case of processes with many massive strongly-interacting particles, eq. (E.3) may be a more convenient alternative to eq. (4.29), since it will reduce the number of independent integrations. We did not explore this possibility in the present paper, and we plan to do so in the future.
We finally reconsider the problem of the simultaneous integration of the (n + 1)and n-body contributions. In sect. 6.2 we have shown how to split the n-body contribution into a sum of terms, each of which is then associated with a real-emission contribution due to a given partonic process and a given FKS pairs. We did so by unambiguously (but arbitrarily) "constructing" the Born-level processes by removing from the real processes the only gluon (in the cases in which there is a gluon) that is an FKS parton in P FKS . This procedure is clearly inspired by the connection between an (n + 1)-body process and its soft limits.
One can turn this logic around, and construct real-emission processes starting from a given Born process. In doing so, it is apparent that this construction can be achieved by considering a → bc branchings, with a a particle entering the Born process. These branchings are reminiscent of an underlying collinear kinematics, but this holds only at a formal level. For example, if one has a top quark in the Born process, by considering the branching t → tg ones gets a new final state with an extra gluon, which can be interpreted as a real-emission process, in spite of the fact that there is no collinear singularity associated with the tg pair.
At this point, one observes that for a given real-emission process, pairs of particles that lead to an underlying Born structure are in fact FKS pairs. For a given r B ∈ R n and a given r R ∈ R n+1 , after constructing P FKS (r R ) we introduce a generalized Kronecker symbol We remind the reader that two processes are considered identical if one of them can be obtained from the other with a permutation of the identities of final-state particles. It is then a matter of algebra to show that .

(E.6)
This equation suggests a definition analogue to that given in eq. (6.13). There, the n-body cross section dσ (n) was given an (n + 1)-body process as an argument. In this case, we give the (n + 1)-body cross section dσ (n+1) an n-body process as an argument: It is quite obvious that the same procedure can be applied to the degenerate (n + 1)body contributions of eqs. (4.41) and (4.42), also in view of the fact that they have a structure identical to that of the initial-state collinear limits of the (n + 1)-body cross sections. We thus define: dσ (n+1) (r B ) = dσ .
Equation (E.7) (or, which is equivalent for the sake of this discussions, eqs. (E.8) and (E.9)) may appear more complicated than its analogue, eq. (6.13). However, the reader should note that, owing to the presence of the Kronecker symbol, the sum over the (n + 1)-body processes in eq. (E.7) is almost trivial, and it amounts to simply counting the a → bc branchings allowed by particle identities. Given that one will have in any case to construct the sets of FKS pairs in order to be able to perform the computations of real-emission contributions, eq. (E.7) is simple to implement algorithmically -in particular, we stress that the current version of MadFKS does have all necessary ingredients to implement it.
We conclude this session by pointing out a couple of issues. Firstly, the use of eqs. (E.7)-(E.10) in place of those of sect. 6.2 implies that the n-body contributions are computed a smaller number of times. This may lead to a significant difference in CPU time when including the finite term of one-loop origin V (n,1) F IN into MadFKS, whose computation for large-multiplicity final states is quite laborious. Obviously, the procedure of sect. 6.2 can be improved in this respect (for example, by not calling the function that returns V (n,1) F IN every time the Born is computed), but the technique described in this section is perhaps a simpler alternative. Secondly, eqs. (E.7)-(E.10) are required for the matching with parton showers according to the POWHEG formalism. Although not strictly necessary, the same technique can also be used in the context of the generation of S events (which are Born-like) in MC@NLO [14].

F. Integration of matrix elements of given helicities
For large-multiplicity final states, the sum over all helicities is a time-consuming operation. A successful strategy at the tree level is that of replacing such an exact sum with one performed with Monte Carlo methods, which in turn requires the matrix elements be available for any possible fixed-helicity configuration. This is indeed the case in MadGraph, and the aim of this appendix is that of showing that the FKS formalism presented in this paper needs only trivial modifications to handle this situation.
We recall that we compute the physical cross section using eq. (6.15), and let us again consider for the moment the case of e + e − collisions only, which implies dσ (n+1) ≡ 0. The (n + 1)-body contribution is given in eqs. (4.28) and (4.29), while the n-body contribution is the sum (see eq. (4.3)) of the quantities given in eqs. (4.4), (4.5), (4.12), and (4.14). In these equations, the matrix elements appear that are defined in eqs. (3.22)-(3.25) (where in the latter one only the finite contribution is actually used, as specified in eq. (B.2)).
We now state the following fact: the equations we have referred to in the previous paragraph apply without any modifications to the case of a partonic process with fixed helicities, except for the fact that in eqs. In other words, the sum over helicities is now not understood in the matrix element. The (n + 1)-body contribution has the subtracted structure given in eq. (4.38). In eqs. (5.17) and (5.18) we have given the necessary ingredients to construct the soft and collinear counterterms respectively; these formulae however apply to the case in which the matrix elements are summed over helicities. When one considers the soft limit of the (n + 1)-body matrix element for a given helicity configuration, eq. (5.17) is basically unchanged. The colour-linked Born matrix elements which appear on the r.h.s. of that equation are computed with the same helicities as those entering the (n + 1)-body matrix element on the l.h.s., except for that of parton i, which is simply not present in the reduced process, as suggested by the notation r i \ . The eikonal pre-factors are also unchanged, but need be multiplied by a factor 1/2 in the case in which the helicity of parton i is kept fixed; if, on the other hand, such helicity is summed over, no modifications are needed in eq. (5.17).
The case of the collinear limit is slightly more complicated. Indeed, in the reduced process r j⊕i,i \ the helicity of the branching parton j ⊕i can in general assume any value, at fixed helicities h i and h j of partons i and j respectively. Equation The quantity M (n,0) has the same definition as in eq. (D.1), except for the sum over helicities, which is not performed here in accordance to the general rule of this section, given in eq. (F.1). We point out that in the second term on the r.h.s. of eq. (F.2) the sum over the helicity h of parton j ⊕ i is not performed, since that term results from the interference between the h = + and h = − states. This is also the reason why the kernels W do not depend on h. We finally stress that all helicities other than those of partons i and j have the same values in the (n + 1)-and n-body processes that enter the matrix elements on the two sides of eq. (F.2) respectively. The kernels T and W introduced in eq. (F.2) are the fully polarized versions of the (four-dimensional) Altarelli-Parisi and Q kernels respectively. They can be The T kernels relevant to the case of the branching parton with negative helicity are T −h b hc a→bc (z) = T +h bhc a→bc (z) (F. 22) for all types of branchings; we have again used the notationh b = −h b andh c = −h c . The Altarelli-Parisi kernels at fixed helicities in four dimensions can be easily obtained from the T kernels: . In order to see this, one has to make use of the identities that relate Altarelli-Parisi kernels of like and unlike helicities P ++ ab (z, ǫ) = P −− ab (z, ǫ) ≡ P ↑↑ ab (z, ǫ) , (F.26) P −+ ab (z, ǫ) = P +− ab (z, ǫ) ≡ P ↓↑ ab (z, ǫ) , (F.27) which, as the notation suggests, hold in general in 4 − 2ǫ dimensions.
Although not relevant to the applications presented in this paper, we conclude this section by discussing the generalizations of eqs. (4.41) and (4.42) (the degenerate (n + 1)-body contributions) to the case of processes with given helicities. These generalizations can be obtained from ref. [34], where the FKS subtraction method was extended to the case of polarized collisions. We have dσ (n+1) i1 (r; k 1 , k 2 ) = α S 2π h P (0)hh 1 The reduced matrix elements that appear in eqs. (F.28) and (F.29) are defined according to eq. (F.3); as the notation suggests, h is the helicity of parton 1 ⊕ī (or 2 ⊕ī). We are not aware of any results for the kernels K for polarized scattering, which are non zero in PDF schemes different from MS; therefore, when considering fixed helicities we shall limit ourselves to using MS PDFs. The damped Altarelli-Parisi kernels at fixed helicities, P can be obtained from eq. (F.23). That equation, however, cannot be used to obtain P (1)hah b ab , since the kernels T have been worked out only in four dimensions. To proceed, we make use of the standard definitions: P ↑↑ ab (z, ǫ) = 1 2 (P ab (z, ǫ) + ∆P ab (z, ǫ)) , (F.33) P ↓↑ ab (z, ǫ) = 1 2 (P ab (z, ǫ) − ∆P ab (z, ǫ)) . (F.34) The quantities ∆P ab are conveniently introduced in the study of collisions with incoming polarized beams. From refs. [41,42] we obtain ∆P gg (z, ǫ) = 2C A 1 1 − z − 2z + 1 + 2ǫ (1 − z) , (F.35) ∆P gq (z, ǫ) = C F (2 − z + 2ǫ (1 − z)) , (F. 38) after taking into account a finite scheme transformation which affects the O(ǫ) term of the qq kernel [42].