Universal Decomposition of Phase-Space Integrands

One-loop integrands can be written in terms of a simple, process-independent basis. We show that a similar basis exists for integrands of phase-space integrals for the real-emission contribution at next-to-leading order. Our demonstration deploys techniques from computational algebraic geometry to partial-fraction integrands in a systematic way. This takes the first step towards a decomposition of phase-space integrals in terms of a basis of master integrals.


I. Introduction
Theoretical calculations of Standard-Model processes have played an important role in extracting information from experiments at both electron-positron and hadron colliders. This is especially true for calculations in perturbative QCD. These calculations are needed to measure Standard-Model parameters from data, and to confront data with expectations in the search for new physics. Nature has not been so kind to date as to reveal direct evidence of physics beyond the Standard Model at CERN's Large Hadron Collider (LHC), in spite of the tantalizing hints from the existence of dark matter and the plausibility of leptogenetic origins to the observed astrophysical baryon-antibaryon asymmetry.
Experiments are thus proceeding into an era of precision measurements, driven by the prospects of increasing luminosity at the LHC and the increasing success of experiments at overcoming issues arising from multi-event pileup [1]. This path calls for increasing precision in theoretical predictions, in order to minimize the part of theoretical uncertainties in the error budget of measurements and resulting constraints. Within a perturbative expansion, the large and quickly running QCD coupling α s means that efficiently predictions, especially for multi-jet observables, are only qualitative; that next-to-leading order (NLO) predictions are the first quantitative order; and that precision predictions require next-to-next-to-leading order (NNLO) predictions.
NNLO predictions require computations of two-loop amplitudes, a subject which has seen a great deal of progress in recent years. This progress has come on two fronts. Firstly, our understanding of how to compute the master integrals efficiently (e.g. with pentagon functions [2][3][4]) within the differential equations method [5][6][7] has improved. Secondly, finite-field methods [8,9] (see also recent work on extension to p-adic numbers [10]) have made it practical to obtain analytic forms for amplitudes using numerical approaches. In particular, this has facilitated the application of on-shell approaches, such as numerical unitarity [11][12][13], which build on earlier dramatic progress in the computation of one-loop amplitudes [14]. With these technical advances, the community has been able to compute numerous two-loop five-point amplitudes (see refs. [15,16] for recent cutting-edge five-point-

one-mass examples).
NNLO predictions also require the computation of real-emission integrals, and combining them with the two-loop virtual contributions. Each of these contributions is singular. The singularities are in practice always regulated dimensionally. They are manifest in the virtual contributions. In contrast, they arise entirely from the integration over phase-space in the double-emission contributions and partly from such integrations in the single-emission ones.
A direct computation of them appears intractable, and so various strategies have been proposed to isolate the singularities into simpler integrals which are tractable analytically.
These strategies fall into three main groups: slicing, subtraction, and hybrids of the two.
In a slicing approach [17], the integrals over the singular (soft and collinear) regions of phase space are performed analytically using the singular approximations to the matrix element. These integrals are universal because the soft and collinear approximations are universal. The integral over the remainder of phase space is done numerically. The definition of the singular regions can be done using pairs of invariants at the level of color-ordered amplitudes, or using other measures such as N -jettiness [18]. The slicing approach saw early use in NLO predictions of jet production. In a subtraction approach [19][20][21], one uses the universal functions describing the soft and collinear behavior of the matrix elements to build subtraction terms rendering the integrand finite everywhere in phase space. The subtracted integral is done numerically, while the integrals of the subtraction terms are performed analytically over all of phase space. The approach of ref. [19] in particular became the most widely used one at NLO.
A hybrid approach has been used extensively at NNLO, where a subtraction effectively controls the NLO matrix element, and N -jettiness slicing or q T -slicing [22] controls its use in the singular regions. There are also NNLO versions of subtraction approaches, including antenna subtraction [23], and approaches [24][25][26] based on generalizations of NLO strategies.
Deployment of these strategies at NNLO has proven far more difficult than might have been hoped, given the experience gleaned from their successful use and automation at NLO. This has motivated us to search for a new approach, which seeks to tackle the seemingly intractable problem of computing phase-space integrals directly, allowing for numerical computation of cross section while maintaining analytic control of singular terms. We are also motivated by the goal of implementing an approach that will allow computing a fully differential cross section at NLO, echoing the success of the 'matrix-element method at NLO' approach of refs. [27]. In this paper, we take the first steps towards this goal. We explore an approach based on decomposing phase-space integrands into simpler integrals, analogs of master integrands in loop computations. In the present paper, we will consider the problem at NLO.
At this order, we are interested in integrating over the emission of a lone gluon (or the splitting of a gluon into a quark-antiquark pair). Ultimately, we would seek to show that phase-space integrands from arbitrary processes at NLO can be reduced to a limited set, with no integrand term containing more than four denominator factors dependent on this lone gluon, and hence contributing to infrared or collinear singularities. In addition, we would seek to show reducibility of numerators. These simplifications are analogous to the well-known reductions of one-loop integrals to a superset of bubble, triangle, and box integrands.
The phase-space integrands are matrix elements, built at lowest order from squares of tree amplitudes. In this paper, we will take the very first step in this direction by showing that tree amplitudes of arbitrary multiplicity can be decomposed into sums of rational functions where each term contains no more than four denominators dependent on the singular gluon.
This shows that there is a process-independent finite basis of terms in both the tree amplitude and in the squared matrix element. We leave the question of reducing this basis to a minimal one to future work.
In the case of loop integrals, simple algebra and a consideration of appropriately constructed numerators suffice to demonstrate this reduction. The reduction for phase-space integrands is more difficult. We employ a masslessness-preserving remapping of partonic momenta along with techniques from computational algebraic geometry to carry it out.
In the next section, we present the setting for our investigation, that of jet observables in heavy-particle decay. The reader will find a more detailed roadmap to the remainder of the article at its end. The road starts in section III, where we review antenna factorization. In section IV, we discuss some subtleties concerning dimensional regularization and counting of degrees of freedom. In section V, we invert the antenna mapping. We use the inverse to building required changes of variables for later sections. In section VI, we review general aspects of reducing one-loop integrands. We review the standard reduction of integrands with constant numerators in section VII and recast this reduction in the language of algebraic geometry. We present a brief review of computational algebraic geometry in section VII B.
We give two examples of reducing contributions to tree-level amplitudes in section VIII.
In section IX, we introduce the notion of triskelia to classify all different contributions to an amplitude. We survey the reduction of all such different contributions in section X. In section XI, we consider the reduction of integrand numerators, recasting the loop-integral case in terms of computational algebraic geometry, and applying that recasting to tree-level amplitudes. We make a few concluding remarks in section XII.

II. Jet Observables
In collider applications of perturbative QCD, differential cross sections of various observables O built out of jets play an important role. Their prediction depends on a jet algorithm and on phase-space cuts on the jets and and other final-state objects. In this section, we use the decay of a colorless scalar φ as an example to set up the framework for our investigation into simplifying the real-emission contributions. We take it to couple to gluons via a f 0 φG 2 , where G µν is the gluon field strength and the coupling f 0 has dimensions of inverse energy.
The singly differential cross section in a jet observable O built out of n jet momenta J i has the form, The decaying scalar has four-momentum K 0 , and the D-dimensional Lorentz-invariant phase-space measure is, The jets are taken to be massless, as is typical of a modern jet algorithm using a massless recombination scheme such as the E T scheme. We take the jet algorithm to be infrared-and collinear-safe (IRCS). The function C J n imposes cuts on the protojets. The fully differential cross section with respect to all jets is d nD σ d D J 1 ⋯d D J n . It can be expanded order by order in perturbation theory, where the leading-order (LO) term is O(f 0 α n−2 ), the next-to-leading order (NLO) term is O(f 0 α n−1 ), the next-to-next-leading order (NNLO) term is O(f 0 α n ), and so on.
The most common type of jet algorithm is hierarchical. We may distinguish three phases in such algorithms: selecting the partons or protojets to be combined into a new level of protojets; combining the four-momenta of those protojets into the new protojet fourmomentum; and after reaching a stopping criterion, imposing cuts on the protojets to yield the set of final jet four-momenta. The selection procedure is typically implemented via a product of theta functions, while the combination procedure is implemented via a product of delta functions. (We retain the traditional name of 'recombination' for the latter.) In eq (2.1), it is the C J n function that imposes the required cuts on protojet energies and rapidities, while the two first phases are implicit in fully differential cross section d nD σ.
At LO, the fully differential cross section is, is the squared tree-level matrix element for the decay of the colorless scalar carrying K 0 into n partons carrying the given momenta k i . The jet function Rec represents the selection and recombination phases of the idealized jet algorithm, with weight localized on configurations yielding exactly n jets with jet four-momenta {J i } n i=1 , and 0 elsewhere. At LO, the protojets are simply the original partons, and the jet recombination algorithm is given by a product of delta functions δ (D) (J i − k i ). The integrals are then trivial, and the differential cross section is the same as the squared tree-level matrix element, Beyond LO, the jet algorithm will of course contain nontrivial recombinations of partonic momenta.
As an example, the k ⊥ algorithm (aka the Durham algorithm [28]) in the three-parton case would use a distance measure, with the recombination taking the following form at NLO, Here, y is the jet-size parameter, and bold symbols denote spatial vectors.
The observable in eq (2.1) is infrared-and collinear-safe by virtue of being a function of jet momenta resulting from an IRCS jet algorithm. More general IRCS observables, such as the distribution of jet energy radii, are certainly possible. The decompositions we will develop in this paper would apply to them as well.
The expressions in eqs. (2.1) and (2.4) contain objects, such as the fully differential jet cross section d nD σ d D J 1 ⋯d D J n and the recombination function , which are themselves delta distributions. In practical applications, each must be integrated over bins corresponding to those used in an experimental analysis.
One of our goals is to give a practical and analytically tractable definition of the fully differential cross section beyond leading order. This is similar to the goal in developments of the matrix-element method at NLO [27]. Our first task is to recast differential distributions in the form of eq (2.1). Let us first recall the conventional approach to higher-order corrections to distributions. At next-to-leading order (NLO), for example, we will have two contributions, virtual and real-emission. We can write the virtual contribution in a form paralleling the all-orders form (2.1). The virtual contribution is, (2.9) Here, M (1V) n is the interference of one-loop and tree-level matrix elements for the decay of the colorless scalar into n partons. The real-emission contribution is similarly given by, (2.11) This contribution again depends solely on tree-level matrix elements, but this time for the decay of the colorless scalar into n+1 partons. In general, both the virtual and real-emission contributions have −2 singularities (D = 4 − 2 as usual), which will cancel in any physical observable thanks to the KLN theorem [29]. If we are somehow able to calculate the integrals in the fully differential cross sections, then we could cancel singularities point by point in the jet phase space, which is one of our goals.
The integrals in eqs. (2.9) and (2.11) are taken over different partonic phase spaces, which makes cancellation of singularities difficult. We can, however, recast eq (2.11) in a form that makes it easier to discuss cancellations. We begin by changing variables in the real integration to one over n protojetsk i and one 'singular' momentum k r , (2.12) The protojet and singular momenta are all massless. The decomposition of the phase-space measure is exact, Here, d LIPS s indicates the singular phase space, and Jac is the Jacobian from the change of variables. The limits of integration on k r are determined by the original phase-space boundaries. The idea is that (within a slice of phase space) all singularities of an integral are associated to specific regions of the k r integration, without reference to integrations over the remainingk i . We will introduce a concrete form for a suitable remapping below.
We can think of the mapping from the original momenta to the protojets as an alternate, theory-friendly jet algorithm.
The change of variables is (in principle) invertible 1 , so that we can think of the original partonic momenta as functions of the protojet momentak i along with k r , (2.14) The protojet momentak i should be thought of as variables on a par with the partonic momenta at LO: massless jet momenta before the event is passed through jet cuts.
Introducing the protojet momenta, we rewrite the radiative contribution, The jet algorithm here is the theoretical idealization of the experimental jet algorithm.
We can now subtract and add back a function implementing a different jet algorithm, Rec T , to obtain, where, (2.17) The first term is finite, because any two IRCS jet algorithms will necessarily give identical results in the singular limit. So long as both approach the singular limit smoothly, this will give a convergence factor compensating the singularities, because all divergences in Yang-Mills theory are logarithmic. It is important to note that this finiteness holds pointby-point in the jet phase space so long as both algorithms yield massless jets. In contrast to subtraction algorithms [19,20], the subtraction term here is in the jet algorithm rather than the matrix element. It shares the goal of isolating divergences in a simpler analytic structure.
We set aside the first term, to be computed numerically after further remappings whose exploration we postpone to future work.
The divergences of the original expression (2.11) are now captured exactly by the second (R) term. Here, we choose the jet algorithm to be the remapping, with different remappings in different sectors of phase space (indexed by A). We can now perform thek i integrations, to obtain, Our choice of theory-friendly algorithm will also require a partitioning of phase space into sectors. That partitioning does not affect the decompositions we pursue in later sections.
For simplicity's sake, we suppress details of it in this article.
The integration in eq (2.19) will give rise to singularities in , but the KLN theorem tells us that they will cancel against the virtual contributions, so that the sum, is finite for any value of the jet momenta J i . (The J i arguments here can of course be treated as just dummy arguments.) The full NLO corrections require the addition of the finite term set aside above, This quantity could then be integrated over the jet phase space, multiplied by an appro-priate binning delta function, in order to obtain a prediction for the NLO corrections to any jet observable of interest, As in eq (2.1), the integration over the jet phase space could be performed in four dimensions, and could be performed numerically.
This is not the way that subtraction [19,20] or slicing [17]  The matrix-element method at NLO also defines a fully differential jet cross section from real emission. Approaches described in the literature [27], however, make use of either phasespace slicing or subtraction in order to isolate singular contributions, compute them, and cancel them against singularities in the virtual contributions.
Our goal is to lay out a different path. The integral in eq (2.19) is of course divergent, so it cannot be performed numerically. At the same time, it may appear analytically intractable.
Our aim in this paper is to take the first steps towards rendering it tractable.
We discuss in detail the first two keys essential to pursuing this path. The first is the choice of an analytically tractable jet algorithm. We will make use of a 3 → 2 clustering algorithm, based on antenna factorization [31] and harking back to the Arclus algorithm [32]. We will review the 3 → 2 clustering in the next section. In a real-world computation, this will leave us with an additional integral as in eq (2.17) quantifying the difference between this jet algorithm, and an experimental one (typically the anti-k T jet algorithm). As mentioned above, we leave a detailed discussion of the difference contribution to future work, limiting ourselves here to the comment that it is finite bin-by-bin in any distribution because of infrared safety.
The other key to a practical approach is an integrand or integral reduction to a basis.
Such bases are almost universally employed in modern loop amplitude calculations. The delineation of such a basis is one of the key results of this article. We will show that a finite basis exists, and consists of integrands arising from products of tree-amplitude terms with at most four distinct denominator factors. This is the first step to obtaining the analog of showing that an arbitrary one-loop integrand can be written in terms of integrands with five or fewer propagators. At one loop, one can in fact do better for integrals. We leave a refinement of our current results to future work.
The best approach requires making a special (if by now standard) choice of dimensional regulator. We will discuss this choice in section IV.

III. Antenna Factorization Revisited
The 3 → 2 clustering algorithm we will use is based on the antenna mapping [31]. While we cannot cluster two massless partons into a single massless object while maintaining overall momentum conservation in an event, we can cluster three massless partons into two massless objects, which we can think of as protojets. In the usual application, the trio of partons is color ordered, and the first and last are 'hard', while the middle parton is 'singular' in the sense that it is soft or collinear to one of the 'hard' partons. This aspect will require restriction to a slice in phase space, but we can cover all of phase space with appropriate slices. We can and will use the mapping without reference to a color ordering.
If we impose parity and fix an overall (D − 1)-dimensional rotation of the two jets by requiring them to follow the hard-parton momenta in all singular limits, the 3 → 2 mapping is unique up to a single function. We can map three massless parton momenta i, j, and k to two massless protojet momentaâ andb, along with a singular (radiated) momentum k r , We will refer to the momenta k i , k j and k k as the "recombining momenta".
The dimensionless function λ(i, j, k) can be chosen freely, subject only to sensible restrictions on its analytic structure and keeping the other dimensionless function appearing, ρ(i, j, k), real. It must be non-singular in soft and collinear limits, though it need not have a well-defined value in the soft limit. With appropriate special choices, we can recover the Catani-Seymour mapping [19]. The ρ(i, j, k) function is given by, In this equation, K 2 = sâb = s ijk is the total invariant mass of the three partons or equivalently the two protojets into which they are clustered. We obtain eq (3.1) by imposing on-shell conditions and momentum conservation to linear combinations of k i , k j , and k k . Keeping ρ(i, j, k) real implies that, If we want a momentum-independent limit on λ(i, j, k), we must insist that, We can see that eq (3.1) has the correct collinear limit independent of the details of λ(i, j, k), by examining the soft limit and collinear limits. In the soft limit, where k j → 0, and consequently s ij , s jk → 0, we see that, and hence, (3.7) In the collinear limit k i ∥ k j , s ij → 0, and, where K ij = k i + k j . Accordingly, again ρ(i, j, k) → 1, and hence, A similar simplification arises in the other collinear limit, where k j ∥ k k .
The explicit solutions in eq (3.1) contain square roots and rational functions. These are not compatible with the language of computational algebraic geometry we will employ in later sections. That language wants a reformulation in polynomial form. As an example of the class of reformulations we will use, let us re-express the constraints in purely polynomial form. Writing, we have, for the protojet on-shell conditions, and as momentum conservation conditions. We can simplify these equations by solving linear ones, Combined, eqs. (3.11) and (3.12) are 5 equations in 6 unknowns (the câ and cb). It is easy to see that the remaining degree of freedom can be identified with λ. Specifically, we can identify câ ,j = (1 + λ(i, j, k)) 2. This leaves us with two coefficients (câ ,i and câ ,k ) to be treated later as additional independent variables, alongside a pair of equations, (3.14) We note that the linear equation can also be solved, but only in an asymmetric manner; we avoid doing this. Substituting the solutions into eq (3.10), we obtain a compact polynomial representation for the mapping, We will also take the λ function as an additional independent variable.
Let us define two "remapping constraint functions" corresponding to eq (3.14), (3.16) Because the solutions for câ ,i and câ ,k are algebraic rather than polynomial in terms of the invariants, we do not want to solve for them explicitly. We could use them to impose the constraints of eq (3.14) in the form of equations,

IV. Dimensional Regularization
The expressions in section II make use of dimensional regularization. All integrations are in D = 4−2 dimensions. This is of course not suitable for the ultimate numerical integrations we must perform, so we must take the four-dimensional limit of appropriate intermediate expressions.
In the LO truncation of eq (2.1), the fully differential cross section is finite as → 0, so in that expression we can take the limit. The jet cuts eliminate regions where it has singularities, so we can also restrict the integration to four dimensions, and take → 0 in the measure as well. The extra integrations would contribute corrections of O( ) to physical observables; but at the end of a calculation, we can certainly take → 0, whereupon they will vanish. We can then calculate the tree-level matrix element in eq (2.5) using four-dimensional techniques, which are critical to efficiency.
At NLO, the situation is more subtle. The virtual matrix element contains explicit poles in , so we cannot take the limit there. While the -dimensional integrations over the corresponding components of the jet momenta are not needed in order to regulate singularities, they would shift the results for differential cross sections by finite amounts. In a physical quantity, this would be canceled by an opposite shift in the real-emission contributions.
The treatment in section II corresponds to the traditional, "conventional" (CDR) scheme in which all momenta and spins are continued uniformly to 4 − 2 dimensions. A different approach keeps observable momenta in four dimensions, and continues only unphysical momenta -virtual particles, soft or collinear real emissions -to 4 − 2 dimensions. This approach, common to the four-dimensional helicity (FDH), 't Hooft-Veltman (HV), and dimensional reduction (DR) schemes, makes integral bases smaller. It also allows the use of four-dimensional techniques in computing loop amplitudes. We shall denote these variants of dimensional regularization "four-dimensional" schemes below. They replace eq (2.8) by, What is the correct real-emission analog of this virtual differential cross section? It is clear that integrals over parton momenta other than k i , k j , and k k should be done in four dimensions; and clear that the integral over k r = k j should be done in 4 − 2 dimensions. We are left to determine the correct treatment of the protojets kâ and kb and of the partons k i and k k : are there additional integrals that must be included along with d LIPS D s ? Because of momentum conservation, it is clear that at least one of k i and k k must have -dimensional components; but are we obliged to integrate over them?
To gain insight into the correct prescription for matching the four-dimensional treatment of external momenta in the virtual contributions, we can examine the treatment in the slicing [17] and Catani-Seymour subtraction [19] approaches. We can also examine the potentially singular denominators that may arise. In both cases, there is effectively only one degree of freedom regulating the singularities. Both calculations, for example, can be embedded in five dimensions and do not require a six-dimensional embedding space.
Furthermore, if we examine the singular terms in the fully differential cross section d 4n σ V d 4 J 1 ⋯d 4 J n in the integrand of eq (4.1), we will find that they are proportional to the four-dimensional leading-order fully differential cross section rather than the D-dimensional one. Unitarity requires that we match this in the real-emission cross section after integration over the singular emission.
The different arguments above all consistently imply that the protojetsâ andb must be treated as four-dimensional, and must be integrated over four dimensions. Moreover, it is most natural to take the -dimensional components of k r = k j as the regulating integration variables. Momentum conservation means that the -dimensional components of k i and k k cannot vanish; but this degree-of-freedom restriction implies that they cannot be independent of k j either. Each component must be linearly dependent on the corresponding component of k j , so that the trio of partons can be embedded in five dimensions.
The mapping in eq (3.1) is valid when either all momenta are strictly four-dimensional, or when all are D-dimensional. Let us now see how to modify it for four-dimensional regulator schemes. Introduce a notation for four-and -dimensional components respectively, By design, µâ = 0 = µb; we will implement this via additional delta functions. This will give two linear equations relating µ i , µ k , and µ j . That will force both µ i and µ k to be proportional to µ j , and allows us to solve for the former, Putting in the solutions for câ ,{i,j,k} , we find, Equation (3.1) is unchanged, though it effectively takes the form, We can recast eq (4.3) as a pair of polynomial equations in scalar variables by contracting both sides of each equation with µ j , and clearing denominators. This gives us two additional remapping constraint functions, (4.6) (The latter is obtained from the sum of the two equations in eq (4.3).) These two functions will also give constraints in the form, (4.7)

V. Inverse Antenna Mapping
Unlike most reconstruction steps in jet algorithms, the antenna mapping is invertible.
Indeed, we need the inverse mapping, for example, in order to change variables from the partonic momenta to the protojet momenta, or equivalently from partonic invariants to protojet invariants. We can solve for k i , k j , and k k in terms of kâ, k r , and kb, as follows, where we define, and take τ (s 1 , s 2 ) to be, It is obvious in eq (5.1) that k i + k j + k k = kâ + kb and that k 2 j = 0; the reader may verify with a little bit of algebra that k 2 i = 0 = k 2 k as well. These statements hold independently of the precise form ofλ. The precise relation between the forward mapping's λ(i, j, k) and λ(s 1 , s 2 ) is complicated; but if we care about the inverse mapping generally (for example, inside an integral) rather than the exact inverse to eq (3.1), we can choose the dimensionless functionλ freely, subject to sensible constraints on its analytic behavior.
Let us verify the behavior in the soft and collinear limits. In the soft limit, k r , sâ r , s rb → 0; consequently, and the inverse mapping takes the form, In the kâ ∥ k r limit, we may define Kâ r = kâ − k r and then write, Here, sâ r → 0 while s rb → ζK 2 (1 + ζ), and accordingly, (1 +λ(sâ r , s rb )) , We then find for the inverse mapping, For the Jacobian for the change of variables from (k i , k j , k k ) to (kâ, k r , kb), we obtain (see appendix A for details), The only singularities appearing here are associated with τ , and this is also true when changing variables in invariants, (5.10) In these equations, x is a massless momentum other than (i, j, k), and so necessarily fourdimensional.

(5.15)
Its solution is eq (5.3). For later use, we define the right-hand side as a constraint function,

VI. One-Loop Decomposition Revisited
Our goal in the present paper is to show that the integrand arising from any squared tree- n+1 can be written in terms of a (small) set of master integrands, and to delineate such a set. We leave the determination of a minimal set to future work.
We will recast the reduction of one-loop integrands for an n-point process to simpler integrands in a technology which we can apply to the problem of reducing real-emission in-tegrands compatible with the theoretical jet algorithm introduced in section II and specified in sections III and V.
In modern one-loop amplitude calculations, one makes use of a standard integral basis, in which any Feynman integral in a process of interest can be written as a sum of basis or 'master' integrals with rational coefficients (rational in the Lorentz invariants and ).
Such a decomposition is relevant even if we do not make use of unitarity in computing the coefficients of the integrals in the expression for the amplitude.
We first review the initial steps in deriving this basis, focusing on the integrands rather than on the integrals. The well-known derivation relies on two techniques: partial fractioning, and use of dimensional identities. There are two reductions that we must consider, of integral, external momenta suffice to form a basis. This is the case on which we focus first, corresponding to Feynman integrals with five or more external legs. We take all internal lines to be massless, while external legs can be either massless or massive.
Consider first consider the reduction of 'scalar' n-point integrands (i.e. whose numerator is independent of the loop momentum) to simpler integrals. This is in a sense the most conceptually important reduction, as it takes us from a potentially infinite set of basis integrals to a finite set.
We will make use of Gram determinants, where, Empty sums are understood to vanish, and the index runs cyclicly mod 6. We can then rewrite, in which j is implicitly summed over 1 . . . 6.
In this expression, where the cross-out indicates an omitted momentum.
The inhomogeneous term in eq (7.3) is, This appears to be generally positive, as evaluated on randomly chosen configurations of momenta; but we offer no proof that it avoids vanishing. We will assume that it is nonzero.
Using eq (7.3), we can then write the following identity, where the cross-out indicates an omitted denominator. This gives a partial-fractioning (Replacing the five momenta in the lower argument with any other set of five external momenta will of course keep the equation the same.) However, one finds that the new coefficients ω ′ j are related, so that their ratios remain the same up to a term proportional to G(k 1 , k 2 , k 3 , k 4 , k 5 ) -which vanishes. This means that we only have one independent relation for the denominators in this case.
For larger n, we could of course iterate this reduction (or its straightforward generalization to integrals with external masses). Alternatively, we can write down all independent relations directly. We expect to find n − 5 independent relations. Modulo accidental relations, each linear relation should allow us to reduce the number of denominators by one. At the integrand level, without taking into account µ 2 terms (which yield integrals of O( )), we expect at best to reduce to terms with five denominator factors, as in the six-point example above.
We can understand the counting intuitively by imagining that we have rotated the D dimensional loop momentum into five dimensions; including any four independent external momenta, we expect to span a space of five independent quantities, which is exactly what happens in the above example.

A. Reduction Using Algebraic Geometry
In order to analyze real-emission integrands, we will need to make use of the tools of computational algebraic geometry. Before doing so, let us first set new lyrics to an old melody, recasting the reduction in the previous section into the language of algebraic geometry.
In the loop integral case, we take as variables all independent invariants built out of the loop momentum, a set we shall call W . (Momentum conservation gives relations between invariants.) We take as parameters all invariants built solely out of external momenta.
We will be interested in polynomial expressions in the variables, with coefficients that are rational functions of the parameters.
Let us explore the singularities of both sides of eq (7.7). Both sides will have poles when any D j vanishes. What about singularities were all D j to vanish simultaneously? In this case, the right-hand side would have no stronger singularity than having only five D j vanish simultaneously. The equation can only be consistent if the left-hand side, too, has no stronger singularity. This will be the case only if all D j cannot vanish simultaneously. At the same time, were the external momenta in arbitrary dimension rather than restricted to four dimensions, all denominators could in fact vanish simultaneously; so this impediment must depend on the momenta being restricted to four dimensions.
The simplest polynomial equation expressing this restriction is precisely the vanishing of the Gram determinant G 1 . We thus want to show that the simultaneous equations, D j = 0 (j = 1, . . . , 6) , G 1 = 0 , (7.8) have no solution. In the language of ideals, we want to show that the ideal built out of the D j and G 1 is the unit ideal. One way to do this, and computationally the most straightforward, is to compute the Gröbner basis of the ideal. If it is 1 (or a constant), the ideal is indeed a unit one. Moreover, there is always a way of expressing the elements of the Gröbner basis in terms of the original defining polynomials of the ideal, via the cofactor matrix. In the case of a unit ideal, this takes the following form, where the c i are polynomials in the variables W . (This is a polynomial analog to Bézout's identity.) For the hexagon integrand, the coefficients are all constants (rational functions of the parameters alone), In the next section, we review the basics of the computational tools we will use in later sections to partial-fraction real-emission integrands.

B. Algebraic Geometry Review
In this section, we provide a brief review of elements of computational algebraic geometry needed for the calculations in this paper. Many of these ideas have appeared earlier in the context of loop integral reduction, see e.g. refs. [34][35][36]. For a pedagogical introduction to these concepts we refer the reader to refs. [37][38][39].
A fundamental object of interest is the ring of polynomials in n variables x i over a field F, denoted by The field F is often a numeric field, such as the rational or complex numbers, but could also be a field of rational functions. In the context of this paper, the field F is taken to be the field of rational functions of the observable kinematics. Oftentimes we will work with a fixed, numeric phase-space point, in which case F reduces again to a numeric field. The tuple of variables {x 1 , . . . , x n } is taken to be the tuple {λ, τ, s i 1 r , s i 2 r , s i 3 ,r s i 4 r }, where the i k correspond to distinct partonic or protojet momenta. Elements of the polynomial ring R are expressed as linear combinations of monomials which we denote as where the α i ∈ Z ≥0 . It is useful to organize the set of monomials with a monomial ordering, which we denote by ≻. This is a total ordering of the exponents α of the monomials. Given a monomial ordering ≻, any polynomial p has a lead monomial LM ≻ (p), which is either 0 or largest monomial in p with respect to the ordering ≻. A useful class of orderings are the socalled "block orderings". We split the variables x i into two non-overlapping sets {y 1 , . . . , y s } and {z 1 , . . . z n−s }, and to each we associate a monomial order ≻ y and ≻ z . To compare two monomials x α and x β , one first compares the part of the monomial depending only on y variables using ≻ y and in the case of a tie compares the part of the monomials depending on the z variables using ≻ z . Unless otherwise specified, we use the degree reverse lexicographic ordering for variable blocks in this paper.
To answer this question, we introduce the concept of an ideal generated by the h i . This is defined as To answer this ideal membership question, we introduce the concept of the Gröbner basis.
Given J, a Gröbner basis is a different set of polynomials that also generate J. We will denote such a set as For brevity, we will say ∆ H,≻ (p) is the remainder of p modulo H and write this as What is special about a Gröbner basis is that having zero remainder modulo such a basis is in one-to-one correspondence with ideal membership. That is, More generally, we can consider ideals whose associated variety is not empty. It turns out that for such an ideal, one can define a geometric dimension which is derived from the associated variety. Quite naturally, when the ideal is associated to a finite set of points, both the ideal and variety are said to be "zero-dimensional". This geometric notion has a useful algebraic consequence, when considering remainders modulo zero-dimensional ideals. This is sometimes referred to in the physics literature as the "shape lemma" [35]. Let us regard the ring R as an infinite-dimensional F-vector space, spanned by monomials in the variables.
In such an interpretation, if we let J be an ideal and G be a Gröbner basis of J with respect to ≻, the operation of computing the polynomial remainder ∆ G,≻ is a linear map acting on R. For J whose dimension is non-zero, the image of ∆ G,≻ is an infinite-dimensional subspace of R. However, if J is a zero-dimensional ideal, and F is algebraically closed, it turns out that the image of ∆ G,≻ is a finite-dimensional subspace of R. Moreover, the dimension of this vector space is equal to the number of points in the associated variety, when counted with appropriate multiplicity.
c. Co-factors and Syzygies For our purposes, it is important to obtain an explicit set of c i such that the relation eq (7.13) holds. Given a Gröbner basis of J, and a polynomial h 0 in the ideal, polynomial reduction by the Gröbner basis explicitly constructs a set of The matrix A ij is known as the co-factor matrix. In the computer algebra system Singular, A ij can be computed by making use of the command lift. Together with eq (7.22), this gives us an explicit set of polynomials c j such that eq (7.13) holds given by i + d (2) i . Furthermore, if x is any element of R then it is clear that the tuple xd i satisfies eq (7.24). The solutions of eq (7.24) have the structure of a so-called module. It is somewhat useful to consider a module as analogous to a vector space, in both cases elements are frequently represented as tuples of elements of a "scalar" algebraic object.
We call this particular module the syzygy module of {h 1 , . . . , h m }. This syzygy module is a subset of R m and therefore a submodule of R m . Unlike vector spaces, modules do not always have bases. Nevertheless, they do have generating sets and there exist algorithms to compute a generating set of a syzygy modules. We denote a generating set of the syzygy module of {h 1 , . . . , h m } as Syz(h 1 , . . . , h m ).
In analogy to Gröbner bases in polynomial rings, one can also introduce Gröbner bases for modules. This will allow us to systematically move through and understand the set of solutions c i of eq (7.13). Specifically, we will to be able to find as well as elements with desirable properties such as low degrees. To this end, let us label the cartesian basis vectors of R m as {e 1 , . . . , e m }. Naturally any element of R m can be expressed linearly in these vectors with elements of R as coefficients. One can introduce a monomial ordering on R m by extending an ordering ≻ on R. We will make use of the "term over position" extension which we label as ≻ TOP . We say that x α e i ≻ TOP x β e j if x α ≻ x β or x α = x β and j > i. When extending the degree reverse lexicographic ordering, this can be implemented in the Singular Mathematica interface as {DegreeReverseLexicographic, ModuleDescending} [40].
With a module term ordering, one can once again define a reduction operation, which we will denote in the analogous manner to eq (7.17). Furthermore, one can also define a Gröbner bases with respect to ≻ TOP and such Gröbner bases can be computed in computer algebra systems such as Singular. Importantly, if the module ordering extends a degree ordering, then remainders will have the lowest possible degree.

VIII. Reduction of Numerator-Free Real-Emission Integrands
We are ultimately interested in squared real-emission matrix elements in Yang-Mills theory. In modern approaches, these would be computed at the amplitude level using spinor variables, and then squared analytically or numerically. We will follow the same logic as in the general one-loop reduction, and instead perform a gedanken calculation of the integrand.
Each amplitude can be written as a sum of cubic tree diagrams; we can replace four-point vertices with a connected pair of three-point vertices and a numerator factor canceling the additional denominator. To study the reduction of numerator-free integrands, it is then sufficient to study diagrams in φ 3 theory.
Any term in the squared matrix element arises from the product of two tree diagrams.
The analog to a scalar loop integral is given by the product of two massless cubic diagrams.
Let us begin by considering the reduction of single diagrams before taking the products that appear in the squared matrix element. The following discussion is sufficient to show that a finite number of collections of denominators can occur in any real emission integrand.
Given the scheme of dimensional regularization we employ (see Section IV), in realemission contributions, the post-remapping singular momentum k r can again be rotated into five dimensions; but it has only four independent components, unlike the loop momentum, because it is on shell. The counting of denominators is then reduced by one; we expect to reduce an expression with five independent denominators to a sum of expressions each with four denominators. We first find five denominators in an eight-point amplitude, so that is the analog of the six-point integral that we will study. Each cubic diagram gives rise to one denominator structure. Many different denominator structures can arise from the ensemble of diagrams, of course. We will classify them and enumerate them in Section IX.
The simplest diagram has all three recombining momenta adjacent in the ordering of legs.
For example, the diagram shown in fig. 1 gives rise to the following expression, 1 s j3 s j34 s j345 s j1345 s j12345 . (8.1) In this example, we take the momenta k 1 , . . . , k 5 to be massless; the recombining momenta k i,j,k are always massless.
Label the denominators, For this denominator, the following Gram determinant, gives rise to an equation analogous to eq (7.6), where, and, While G 2 vanishes on physical configurations because the momenta k 1 , . . . , k 5 are fourdimensional, it does not vanish manifestly when expressed in terms of invariants. In this way, it gives rise to a partial-fractioning identity in eq (8.4).
Other diagrams, however, give rise to denominators where no Gram determinant (nor even a set of Gram determinants involving the partons alone) suffices to obtain a partialfractioning identity. For such diagrams, we need to deploy the machinery of computational algebraic geometry reviewed in section VII B.

A. Reduction Using Algebraic Geometry
Let us now consider the general case of real emission integrand partial fractioning, where we cannot simply rearrange a Gram-determinant identity. An example is the term, 1 s i4 s i45 s j1 s j12 s j123 , (8.7) which arises from the diagram in fig. 2.
In the abstract, with T 1 = s i4 , T 2 = s j1 , T 3 = s i45 , T 4 = s j12 , and T 5 = s j123 , we seek an identity of the form, where the Z j are functions of the variables which vanish on physical configurations. These functions would implement the four-dimensionality of the protojets or the dependence of the protojets on the recombining momenta. As variables, we take a set of invariants dependent on the recombining momenta including the T i , along with the mapping coefficients câ ,i , câ ,k , and λ appearing in eq (3.14). Because we impose the four-dimensionality via Gram constraints, we need 20 invariants, the number independent in D dimensions, rather than the 14 independent in four dimensions.
More concretely, we seek an identity of the form, where the Gram determinants G 3 , . . . , G 6 are,  6). The Grams G 3,4 enforce the four-dimensionality of kâ ,b , and G 5,6 remove otherwise vanishing terms from the c j . The remapping constraints R 1,2 enforce the on-shell conditions for kâ ,b , and R 3,4 implement the appropriate variant of dimensional regularization. Unlike the case of eq (8.4), we cannot obtain such a decomposition if we insist that the c j be functions of the nonrecombining invariants alone. We would however expect to obtain an identity of the form of eq (8.8) if we allow the c j to be polynomials in the variables.
In order to find such a decomposition, we would need to show that the joint equations, 1, . . . , 4) , (8.11) have no solution. With presently available Gröbner-basis algorithms and implementations, this problem appears hopelessly intractable if attempted in a fully analytic way. We use the computational algebra system Singular, but we expect the performance of other available systems to be similar.
In order to accomplish our goal, we must simplify the calculation. We do so as follows: we switch to variables built entirely out of kâ, kb, and k r in addition to the partons k 1,...5 ; we solve the Gram-determinant constraints explicitly; we use finite-field numerical values for all momenta except k r ; and we use modular arithmetic in Singular, which yields vastly faster Gröbner-basis computation than with rational arithmetic.
In this approach, we should think of the recombining momenta k i , k j , k k purely as functions of kâ, kb, and k r via the inverse mapping eq (5.1). The singular momentum k r is the only variable momentum; all others are held fixed to specific finite-field values. That is, we do not choose values for the recombining partons k i , k j , and k k ; nor do we choose a specific form for theλ function in the inverse mapping. We will take as our variables a set of independent invariants built out of the D-dimensional momentum k r , along with the symbols τ andλ. These take the place of the set W in the one-loop integrand reduction.
Choosing numerical values for k 1,...,5 and kâ ,b satisfies the constraints G 3 = 0 = G 4 automatically; we can solve the constraint G 6 = 0 by choosing a set of four invariants, for example is D-dimensional, there is no further relation between these invariants arising from k 2 r = 0.) The vanishing of G 5 then follows as well.
The several R j are replaced by the lone constraint R τ = 0, with R τ defined in eq (5.16).
In addition, we must supply a polynomial constraint equation fixing the functional form of λ; a convenient choice is, Rλ ≡λ(s râ + s rb ) + s rb − s râ , (8.13) withλ defined by Rλ = 0. Other choices are possible, but we find that this one gives the simplest identities. It also simplifies the Jacobian (5.9). As variables for computing the Gröbner basis, we take the four variables in W r , along with τ andλ, V = W r ∪ {τ,λ} = {s râ , s rb , s r1 , s r2 , τ,λ} . (8.14) Our ideal is then ⟨B⟩, where, One might wonder if this setup is D-dimensional, as we only include variables with a fourdimensional dependence on k r . However, we also do not include the on-shell condition for k r in B. This is consistent as, due to the spherical symmetry beyond 4-dimensions, the dimensional components of k r can be used to explicitly solve the on-shell condition. Hence the setup is indeed D-dimensional.
Recalling the discussion of section VII B, we see that the first step to finding a partial fractions identity is to compute the Gröbner basis, GröbnerBasis B; V ; (8.16) in Singular, we use the slimgb function for GröbnerBasis . We use a block ordering with the four invariants forming one block, and {τ,λ} the other, and reverse degree lexicographic (degrevlex) within each block.
We find that the Gröbner basis is indeed {1}, which shows that the desired equations have no common solution, and that a partial-fractioning identity then exists. To find the desired identity, we make use of the fact that the cofactor matrix relates the original generating set to the Gröbner basis. We can obtain the cofactor matrix C using Singular's lift function. By definition, This offers us an identity of the form, This is of the desired form, and already allows us to reduce terms with five-factor denominators to a sum of terms with four-factor denominators (as the R τ and Rλ terms will ultimately vanish on physical configurations).
Nevertheless, this identity is not optimal as the C j (for j = 1, . . . , 5) are not necessarily independent of the variables in W r . This risks introducing powers of k r into the numerators.
To address this issue, we recall from section VII B that C is in fact only well-defined up to syzygies of B. If z is a syzygy of B, that is if z ⋅ B = 0, then, It is therefore clear that we can shift C to reduce or eliminate powers of k r in its first five components. To do this systematically, we again make use of Gröbner basis techniques.
Specifically, the remainder of polynomial division by a degree ordered Gröbner basis, is always of the lowest possible degree. Therefore, in order to obtain the "simplest" form for C, we will reduce it against a modified Gröbner basis of the syzygies of B. The remainder will be the desired simplest form.
We first compute the syzygy module of B, which we can do in Singular using the syz function, We now wish to use this syzygy module in a Gröbner basis computation to find an alternative choice for C, with lower polynomial degree. A subtlety in implementing this is that we are uninterested in the values of C 6 and C 7 as they multiply the constraint functions R τ and Rλ which are zero on physical configurations. Nevertheless, we find that the value of {C 1 , . . . , C 5 } with the lowest polynomial degree does not come together with the {C 6 , C 7 } with the lowest polynomial degree. For this reason, we are forced to instruct the Gröbner basis algorithm to throw away C 6 and C 7 . To this end, we modify the syzygy module by adding two syzygies which freely change the last two components of B, corresponding to R τ and Rλ, {(0, 0, 0, 0, 0, 1, 0), (0, 0, 0, 0, 0, 0, 1)} . (8.21) In any degree ordering, this has the effect that the components C 6 and C 7 will always reduce to zero, and have therefore been discarded within the computation. This procedure could be phrased more elegantly and formally in terms of quotient rings.
We now compute the Gröbner basis of the modified syzygy module, In computing this Gröbner basis, we make use of a term over position module ordering which extends a block ordering with the four invariants forming one block, and {τ,λ} the other.

IX. Triskelia
In the previous section, we discussed an algorithm to construct partial-fraction identities for a given denominator structure. We turn now to the task of systematically organizing the set of denominators that arise in a tree-level amplitude. We focus on an organization that will allow us to characterize all different ways a denominator can depend on the singular momentum when using the inverse antenna mapping (5.1) to express the recombining momenta in terms of the protojet and singular momenta. This will then allow to systematically collect partial-fraction identities for all possible denominators in a tree-level amplitude. This dependence can be graphically organized by combinatoric devices that we dub triskelia. The best-known triskelion in the wider culture is probably the Isle of Man's.
Our starting point is to consider the complete set of tree-level Feynman diagrams for an n-point amplitude in a cubic theory. As discussed earlier, this suffices for Yang-Mills. The external legs are the recombining momenta k i,j,k as well as n − 3 numbered momenta. Using momentum conservation, each propagator denominator can be rearranged to depend on at most one of the three recombining momenta. The denominator factors independent of them are simply constants, and we can pull them out front as overall prefactors.
We can re-express the recombining momenta in terms of the protojet momenta kâ ,b and the singular momentum k r using the inverse antenna mapping (5.1). At this point, each denominator factor dependent on a recombining momentum depends implicitly on the singular momentum. (The denominator factors independent of the recombining momenta are also independent of the singular momentum.) In any given Feynman diagram, the denominator factors dependent on the singular momentum can therefore be organized into three groups: one for each recombining momentum. Analogously to similar practices when considering loop diagrams, we will draw only the propagators corresponding to these denominators, "pinching" the others. In any such pinched diagram, each of the three groups forms a line, and the three lines meet at a single "central vertex". We can thus visualize a denominator structure in terms of a three-pronged object, that is a triskelion. Examples of such triskelia can be found in figs. 3 and 4.
In the following sections, we will investigate the partial-fraction decompositions by enumerating the identities for all distinct triskelia with a fixed number of denominators. As an aid in this endeavor, let us count the number of kinematically distinct triskelia with n de-nominators. Each of these denominators must belong to one of the k i , k j or k k denominator groups. We begin by counting the number of ways to partition the n denominators into three groups. By a standard "stars and bars" argument 3 , it is easy to show that for d groups this is n+d−1 d−1 . There are therefore n+2 We denote a triskelion by the notation (n j n i , n k ), where n j is the number of external momenta attached directly to the same central-vertex line as k j , and similarly n i,k count the number of external momenta attached to the same lines as k i,k respectively. This notation does not distinguish whether the external momenta are massive or massless, nor the number or type attached directly to the central vertex. We can reduce the number of triskelia to be considered by relying on the n i ↔ n k symmetry, and choosing n i ≥ n k . For n = 5, there are 1152 triskelia that we must consider.

X. Survey of Numerator-Free Real-Emission Integrand Reductions
In section VIII, we gave an example of a five-factor denominator where we needed the machinery of computational algebraic geometry in to demonstrate its reducibility to a sum of four-denominator terms. We classified the distinct types of denominators we must consider for general tree amplitudes for any number of external legs in the previous section, using the notion of triskelia. Here, we summarize the results for all the different triskelia that arise.
It is sufficient to examine the configurations for five-factor denominators. When examining terms in amplitudes with more than five factors in denominators -that is, when examining real-emission amplitudes with more than eight external legs -we can always choose a subset of the denominator factors to reduce first. In general, configurations with 3 The partitioning of n objects into d groups is the same counting problem as finding all monomials of degree n in d variables, hence the counts agree.
(a)  four or fewer denominators cannot necessarily be partial-fractioned in a D-dimensional way.
Each reduction step will eliminate factors in the denominator, effectively yielding a term in a lower-point amplitude but with external legs massive rather than massless. We must therefore consider all possible configurations of massive or massless external legs, and also different configurations of external legs attached to the 'central' vertex: no leg, a massless leg, or a massive leg. In fig. 3, we display graphically the set of five-denominator skeleton triskelia. As it turns out, the key aspects of the algebraic geometry are independent of Triskelion r-indep. mass-indep. max.λ power max. τ power (1 2, 2) ✓ ✓ 3 4   TABLE I. Reduction of terms with five-factor denominators, for all distinct triskelia. The second column indicates whether all coefficients are independent of k r ; the third column whether that is independent of the masses of the fixed external legs and the momentum flowing into the central vertex. The fourth and fifth columns show the maximum powers ofλ and τ in the partial-fraction coefficients.
these details. The results of our case-by-case analysis are shown in table I. For all triskelia, we obtain the desired reduction to a sum of four-factor denominators, independent of the masses of the external legs. The resulting numerators are free of k r in the form of s xr , but do depend implicitly on k r throughλ and τ . The latter dependence can be thought of as merely changing the functional form of the Jacobian (5.9). The maximum powers ofλ and τ in the coefficients of the resulting four-denominator terms are also shown in the table.

XI. Reduction of Integrand Numerators
In sections VII and VII A, we examined the reduction of numerator-free loop integrands in a conventional approach and also as recast in the language of computational algebraic geometry. Of course, in seeking reductions to a basis for virtual contributions, we must also consider loop integrands with numerators dependent on the loop momentum. In general, the loop momentum will be contracted into either external momenta, external polarization vectors, or other four-dimensional vectors. For any given integral with five or more external legs, we can express vectors that aren't external momenta in terms of a standard basis of four external momenta. The numerator is then a polynomial in invariants of the loop momentum contracted with external momenta. Its degree is determined by the power counting of the theory; for Yang-Mills theory, we must consider polynomials of up to degree n for an n-point amplitude.
The conventional approach observes that we can write any of the invariants, The first two terms cancel denominators, giving rise to lower-point integrals, while the last two terms yield integrals with a lower-degree numerator. Repeating this partial-fractioning ultimately reduces any integral either to a numerator-free integral, or to a four-(or lower-) point integral.
In order to rephrase this argument in the language of algebraic geometry, consider the integrand of a five-point one-loop integral with external momenta in four dimensions. The denominators again take the form given in eq (7.1), but there are only five instead of six of them.
We consider an integrand of the form, where Poly is a polynomial. We take as variables the loop momentum contracted with four specified external momenta {k 1,2,3,4 }, and re-express the polynomial in terms of those variables. We label this set of variables W ∶4 . Invariants of external momenta are parameters, that is constants with respect to the variables. The statement of reduction is then that the remainder of the polynomial with respect to the Gröbner basis of the denominators is constant. For this to be true for any polynomial, it is necessary and sufficient for, For this to hold, here the Gröbner basis must of course be nontrivial, and the underlying ideal must be zero-dimensional. Zero dimensionality implies a finite number of solutions; here there is a unique solution, so that no Lorentz invariants of the loop momentum can appear on the right-hand side. In the loop case, there is a unique solution in appropriate variables, but this is not essential to the reduction.
Equation (11.3) allows us to write a decomposition of the form, 4) with ω N 0 independent of , when Poly 1 is linear in . In the loop case, the coefficients ω N j also turn out to be independent of .
What about the real-emission case? Instead of the five denominators we considered in section VIII, we now consider terms with just four denominators {T i } 4 i=1 , using the same set of variables V given in eq (8.14). We require the analog of eq (11.3), namely that, v mod GröbnerBasis ({T 1 , T 2 , T 3 , T 4 , R τ , Rλ}; V ) = Poly(τ,λ) ∀v ∈ {sâ r , s rb , s r1 , s r2 } .

(11.5)
We must again require that the Gröbner basis be nontrivial, and that the underlying ideal be zero-dimensional. We again check this equation for all different configurations; here that is given by the triskelia with four denominators. These are shown in fig. 4. We find that this condition is obeyed by all, as shown in Table II. This allows us to write the decomposition, 6) where ω N,r 0 is a polynomial in τ andλ, with maximum powers given in Table II for the different triskelia. Unlike the loop case, however, we cannot always ensure that the ω N,r j are independent of k r (even after reducing the coefficients with respect to the Gröbner basis of their syzygies, following the same discussion as below eq (8.19)). The r-independence or dependence of the coefficients is shown in the last column of Table II. In eq (11.6), we have dropped terms proportional to R τ and Rλ as these functions will vanish on physical configurations.
Triskelion r-indep. mass-indep. max.λ power max. τ power r-indep. coefficients is independent of k r ; the third column whether that is independent of the masses of fixed external momenta and the momentum flowing into the central vertex. The fourth and fifth columns give the maximal powers ofλ and τ in ω N,r 0 . The sixth column indicates whether the coefficients ω N,r j are independent of k r .

XII. Discussion and Conclusions
The modern machinery of loop calculations relies on the existence of a basis of master integrals for each process. At one loop, this basis is well understood, and the integrals that arise are process-independent. The existence of the basis both at one loop and at higher loops depends in part on partial-fractioning identities for integrands of general Feynman integrals.
In this article, we have taken the first step in constructing an analogous basis for the realemission contributions. We studied the tree amplitudes out of which these contributions are computed at next-to-leading order. We showed the existence of partial-fractioning identities reducing any term in any tree-level Feynman diagram to denominators with four or fewer factors. Furthermore, each term will have either a numerator free of explicit dependence on the emitted particle, or a denominator with three or fewer factors. We also showed how to construct the required partial-fractioning identities explicitly.
Our construction relies on the inverse of the antenna mapping, expressing a color-ordered set of three momenta sufficient to describe all soft and collinear singularities in a slice of phase space, in terms of physical proto-jets and an emitted momentum. We employ techniques from computational algebraic geometry in order to find the partial-fractioning identities.
The work described here leaves to future work the question of delineating the form of the basis for the complete integrand of the real-emission contribution. This contribution is given by the square of the tree amplitude. In addition, the form of the master-integral basis at one loop also depends on Lorentz-invariance identities, related to the vanishing of integrals of certain nontrivial integrands. The form of the basis at higher loops relies on the vanishing of a broader class of integrals, generated by total derivatives. We leave the extension of these ideas to phase-space integrals to future work as well.

A. Jacobian for Change to Protojet Variables
There are two factors in the change of variables from the partonic momenta {k i , k j , k k } to the protojet momenta {kâ, k r , kb}: that contributed directly by the measure d D k i , and that contributed by the on-shell delta functions. Because the singular momentum does not change form (k j = k r ), its Jacobian is 1, and we can focus on the remaining variables. Writing the partonic variables in the form, k x =c x,x kx +c x,r k r , where x ∈ {i, k} andx =â,b, the partial derivatives themselves split up into two terms, ∂k µ x ∂k ν x =c x,x η µ ν + ∂c x,ŷ ∂sẑ r k μ y + ∂c x,r ∂sẑ r k µ r ∂sẑ r ∂k ν x =c x,x η µ ν + 2 ∂c x,ŷ ∂sx r k μ y + ∂c x,r ∂sx r k µ r k r ν .
The form of the second term relies on the fact thatc x,x andc x,r depend only on K 2 , sâ r , and s rb , with K 2 taken to be constant. The form of the second term means that it is a rank-1 matrix, so that the determinant is linear in it. Define, with I denoting the x and µ indices (taken tensorially), and J thex and ν ones. The The inverse of A is, The contraction with η implies that the second term in the expression for B in eq (A3) = w D 0 (sâ r , s rb )τ D (sâ r , s rb ) 1 + w −1 0 (sâ r , s rb )τ −1 (sâ r , s rb )sŷ r × c k,b ∂c i,ŷ ∂sâ r −c i,b ∂c k,ŷ ∂sâ r +c i,â ∂c k,ŷ ∂s rb −c k,â ∂c i,ŷ ∂s rb . (A6) Using the relations arising from momentum conservation (5.12), the first factor in the Jacobian simplifies to, w D 0 (sâ r , s rb )τ D (sâ r , s rb ) 1 + w −1 0 (sâ r , s rb )τ −1 (sâ r , s rb ) sŷ r ∂c i,ŷ ∂sâ r − ∂c i,ŷ ∂s rb .
We have checked the product of eqs. (A11) and (A14) numerically.

B. Finite-Field Momentum Configuration
For the calculations described in sections VII A, X, and XI, we used the following momentum configuration, in the finite field Z p , with p = 275604541.