Planar Two-Loop Five-Parton Amplitudes from Numerical Unitarity

We compute a complete set of independent leading-color two-loop five-parton amplitudes in QCD. These constitute a fundamental ingredient for the next-to-next-to-leading order QCD corrections to three-jet production at hadron colliders. We show how to consistently consider helicity amplitudes with external fermions in dimensional regularization, allowing the application of a numerical variant of the unitarity approach. Amplitudes are computed by exploiting a decomposition of the integrand into master and surface terms that is independent of the parton type. Master integral coefficients are numerically computed in either finite-field or floating-point arithmetic and combined with known analytic master integrals. We recompute two-loop leading-color four-parton amplitudes as a check of our implementation. Results are presented for all independent four- and five-parton processes including contributions with massless closed fermion loops.


Introduction
The progress in our understanding of the analytic properties of loop amplitudes has recently led to the computation of the first two-loop five-point amplitudes in QCD [1][2][3][4][5][6]. These computations focused on the leading-color contributions to the five-gluon process. In this paper we take a further step and compute the scattering amplitudes of all five-parton processes in the leading-color limit, including corrections with massless closed fermion loops. Two-loop five-parton amplitudes without closed quark loops were recently presented in ref. [7], and related work on the complete reduction of two-loop five-parton amplitudes appeared in refs. [8,9]. Our results are an important step towards the automation of the calculation of two-loop partonic amplitudes, which are in turn an important ingredient towards obtaining theory predictions to three-jet production at hadron colliders at next-tonext-to-leading order (NNLO) in QCD.
The main result of this paper is a numerical method for the computation of two-loop multi-parton amplitudes, including massless quark states in addition to gluons. We apply a numerical variant of the unitarity method [10][11][12][13], which was extensively used for oneloop computations [14][15][16][17] and recently generalized to two loops [18][19][20]. In this paper we extend its implementation to two-loop processes involving fermions. Four-parton twoloop corrections interfered with their tree-level amplitudes [21][22][23][24] as well as the associated helicity amplitudes [25][26][27][28][29] were computed with analytic methods some time ago, and we recompute these results as a check of our implementation. The two-loop numerical unitarity method we employ avoids the challenging algebra of analytic multi-scale computations and is at the same time sufficient for the numerical phase-space integration required in crosssection computations. To showcase its potential we provide numerical benchmark values for five-parton amplitudes. To this end we compute integral coefficients with exact or floating point arithmetic, and combine them with the numerical evaluation of the two-loop master integrals [30,31]. We leave an analysis of the integration over the physical phase space to future work.
A number of developments is necessary for handling fermions. Fermion amplitudes have been computed within the numerical unitarity method at one loop [32][33][34][35][36], and in this paper we propose a generalization for two-loop amplitudes. We first discuss the treatment of fermion states in dimensional regularization. At higher loop orders, the subtleties in this procedure become increasingly relevant (see [37] for a recent review). In particular, we discuss in detail the definition of dimensionally-regulated helicity amplitudes with pairs of external quarks. We present a prescription for how to define a helicity amplitude at two loops which can be used to compute interference terms in the 't Hooft-Veltman scheme. This prescription can be implemented numerically, extending well known analytic methods [25][26][27][28][29] which are not directly applicable in a numerical calculation since they rely on abstract algebraic manipulations of the γ-matrix algebra.
A second technical advance concerns the implementation of the numerical unitarity method in exact arithmetic [6] for amplitudes with fermions, based on the use of finitefields techniques for amplitude computations [38]. The main obstacle to overcome relates to the fact that generic polynomial equations do not have solutions in a generic algebraic field. This is in tension with the fact that in a unitarity-based approach one needs to generate loop momenta satisfying the set of quadratic on-shell conditions. We describe a way to handle this difficulty while maintaining the power of exact finite field computations when considering both gluons and fermions. We stress that the ability to perform exact calculations on rational phase-space points is an additional feature of our computational framework, and the numerical unitarity method has also been implemented in more standard floating-point arithmetic.
The article is organized as follows. In section 2 we define leading-color dimensionallyregulated helicity amplitudes with external fermions. In section 3 we present explicit details of our implementation. In section 4, we first present results for the leading-color four-parton helicity amplitudes, and then present our numerical results for a complete set of independent five-parton amplitudes. We also describe the checks we performed. We give our conclusions and outlook in section 5. Finally, we present some useful γ-matrix identities in appendix A, as well as some details on the infrared structure of the amplitudes in appendix B.

Dimensionally Regulated Helicity Amplitudes
Dimensionally regulated scattering amplitudes are functions of the continuous dimension parameter D = 4 − 2 , which regulates both ultraviolet and infrared singularities of loop integrals. Within this framework, the state spaces of the external particles are also naturally understood to be formally infinite dimensional. In practice, however, we are interested in obtaining predictions for physical external states that are strictly four dimensional. For instance, in this paper we will be computing helicity amplitudes. In order to compute a dimensionally-regulated amplitude with a given set of four-dimensional external states, we must find how to represent them in the D-dimensional space. That is, we must find a consistent embedding of the physical four-dimensional state in D dimensions. This is trivially achieved for gluon helicity states: these are vector particles and any four-dimensional polarization state can be embedded in a generic D-dimensional space by filling the remaining components of the vector with zeros. For fermion states, however, the embedding is less trivial as the nature of the D-dimensional Clifford algebra means that there is no single associated state in D-dimensions. As such, one might wonder if it is possible at all to unambiguously define four-dimensional helicity amplitudes with external fermions. In this section we describe how we address this problem, inspired by the approach of refs. [26,28], and precisely define the objects that we will be computing in subsequent sections.

Embedding of Fermionic States in Dimensional Regularization
There are several consistent regularization schemes that can be chosen, see e.g. ref. [37] for a recent review, and our discussion applies to both the conventional dimensional regularization (CDR) and the 't Hooft-Veltman (HV) schemes. The two schemes differ in the way vector particles (gluons in our case) are treated, and we follow the description given in the reference above. 1 In CDR, all vector fields are vectors in a space of dimension D s = 4 − 2 .
In HV, one distinguishes between regular and singular vector fields. The former do not lead to any singularities and are considered to be strictly four-dimensional objects. The latter are a source of singularities and are D s -dimensional vectors. For our purposes, this means that gluons whose momentum we integrate over are D s dimensional, and external gluons whose momentum we do not integrate over are four dimensional. The two schemes are consistent, in that their contributions to NNLO computations can be related by known transition rules [39], but as we shall see below the HV scheme introduces some simplification in the calculation.
We consider fermions in D s dimensions, as is for instance necessary when a CDR gluon or a singular HV gluon is emitted from a quark line. If the fermion line closes upon itself, as in e.g. the N f corrections to gluon amplitudes (i.e., corrections with a closed massless-quark loop), we only need the defining property of the Clifford algebra where we explicitly write the dimension D s as a subscript of the γ-matrices and the metric, and use a metric with mostly-minus Minkowski signature, g [Ds] = diag{1, −1, ..., −1}. Here 1 [Ds] is the identity operator in the representation space of the Clifford algebra. In the presence of external fermions, however, we must also describe the corresponding states and an explicit representation of the D s -dimensional Clifford algebra is required. Since we are ultimately interested in specifying four-dimensional external states, it is furthermore convenient to construct the representation in a factorized way starting from four dimensions (see e.g. refs. [40,41]). We thus consider a Clifford algebra in D s dimensions as the tensor product of a four-dimensional and a (D s − 4)-dimensional one: whereγ [4] ≡ i(γ 0 [4] γ 1 [4] γ 2 [4] γ 3 [4] ), such that (γ [4] ) 2 = 1 [4] is the identity operator in the fourdimensional algebra. The indices a, b denote the spinor indices in the four-dimensional algebra and κ, λ the ones of the (D s − 4)-dimensional one. The γ µ [Ds −4] form themselves a (D s − 4)-dimensional Clifford algebra with signature g [Ds−4] = diag{−1, ..., −1}. In amplitude calculations we naturally encounter products of γ matrices, and in this paper we will mainly focus on chains of γ [Ds−4] matrices. We thus define a convenient basis for these chains, constructed by anti-symmetrizing over their Lorentz indices and given by (see e.g. [40]) where S n denotes the set of all permutations of n integers and sgn(σ) the signature of the permutation σ ∈ S n .
The spinor states associated with the D s -dimensional Clifford algebra live in a D tdimensional space. 2 For four-dimensional momenta they can be constructed from a tensorproduct representation as where we have introduced an index s = {h, i} to denote the polarization states in terms of spinors of the four-and (D s − 4)-dimensional subspaces. Without loss of generality we can require that (η i ) κ and (η i ) κ be dual to each other, where the {f, f ν 1 , f ν 1 ν 2 , · · · } are the constant coefficients in the decomposition. In loop calculations, this naively introduces reference vectors and tensors that can yield linear dependence on the components of the loop-momenta beyond four dimensions, see e.g. [34,36], in contrast with what happens, for instance, in the case of amplitudes with only gluons. We only mention this here as an observation since, as we shall see in the remaining of this paper, this will not be an issue with our definition of helicity amplitudes. The tensor product representation of the Clifford algebra is particularly useful to separate four-and (D s − 4)-dimensional spinor indices in γ-matrix chains. Indeed, a product of γ matrices where some Lorentz indices are within four dimensions (denoted µ i ), and some are beyond four dimensions (denotedμ i ) is split into two blocks, a four-dimensional and a (D s − 4)-dimensional one. For instance, we have (2.7) Consider now contracting the above product of γ matrices with a four-dimensional fermion state, such as the u andū spinors: The result is a tensor with open indices in the (D s − 4)-dimensional space. We recall that these are in one-to-one correspondence with a (D s − 4)-dimensional state, and the above expression (2.8) is thus equivalent to a contraction with on-shell helicity states in D s dimensions. Our ultimate goal is the calculation of amplitudes relevant for cross-section computations, and we must then understand which tensor structures beyond four dimensions are necessary. This will be done in the next subsections. We finish this section with two comments. First, we note that in the HV scheme this tensor decomposition results in simpler expressions than in CDR. Consider for instance the tree-level qq → QQ amplitude In the HV scheme, the gluon between the two quark lines is four dimensional, i.e., µ ≤ 3, while in the CDR scheme, the gluon is D s dimensional. From eq. (2.2) we thus get where the M (0) i are coefficients that are determined from products of four-dimensional γmatrices contracted with four-dimensional spinors. In the HV scheme the amplitude is determined by a single coefficient, while in CDR two are needed. In the remainder of this paper we thus choose to work in the HV scheme. Nevertheless, our discussion generalizes to the CDR scheme in a straightforward way.
The second comment we wish to make is that, although in this section we consider D s = 4 − 2 , which means the Clifford algebra defined in eq. (2.2) is infinite dimensional, in numerical calculations one might need to construct an explicit representation of the Clifford algebra and thus take D s to be an even integer (larger than 4). The construction of eq. (2.2) still holds and, in fact, it can be iterated: any even D s can be reached by constructing a tensor product of the (D s − 2) algebra with a 2-dimensional algebra, even if the (D s − 2) algebra was already constructed as a tensor product of two algebras.

Tensor Decomposition of Helicity Amplitudes
We consider a helicity amplitude M , expanded in perturbation theory, with the k-th order term written as M (k) . We saw previously that these are tensors in the (D s − 4)-dimensional spinor space, see eq. (2.10) for an explicit example. Here we introduce a basis for the associated tensor space in the spinor indices beyond four dimensions, whose elements are denoted as v n . In general, the basis depends on the physical process described by M and on the order k in the perturbative expansion. We will suppress this dependence for simplicity of the notation and write where the M (k) n are computed from γ µ [4] matrices and external states in four dimensions, and the tensor structure of the amplitude in the spinor indices beyond four dimensions is fully contained in the v n . In the following, we explicitly construct the basis {v n } for two families of amplitudes: those with a pair qq of external quarks and any number of external gluons, and those with two pairs qq and QQ of external quarks (of either different or identical flavor) and any number of external gluons.
The different tensors v n are constructed by contracting the Lorentz indices of chains of γ [Ds] matrices with other Lorentz vectors in the amplitude after all loop integrations have been performed. The remaining objects that carry Lorentz indices are four-dimensional external momenta, four-dimensional polarization vectors and chains of γ [Ds] matrices. Any Lorentz index in a γ [Ds] -matrix chain that is contracted with a four-dimensional object becomes four-dimensional, contributing only a trivial tensor structure in the (D s − 4)dimensional space. For instance, if ε µ represents a four-dimensional polarization vector of an external gluon, (2.12) Similarly when two Lorentz indices are contracted inside the same chain of γ [Ds] matrices, the tensor structure beyond four dimensions is trivial, as follows from: Non-trivial tensors v n are obtained by contracting Lorentz indices of two chains of γ  matrices. The basis introduced in eq. (2.3) for these chains is particularly useful for computing these contractions. Let us consider an amplitude with a pair qq of external quarks and any number of external gluons. There is a single chain of γ  matrices and, as there are no other objects with (D s − 4) indices, it follows from the discussion above that for this case there is a single term in the sum of eq. (2.11): (2.14) We define the dual tensor w 0 such that w 0 · w 0 = 1, with more details given in appendix A. Let us now consider an amplitude with two quark pairs of different flavors, qq and QQ, and any number of gluons. We can now contract Lorentz indices between two different chains of γ matrices, and the basis {v n } is then larger in this case. Using the basis for the γ-matrix chains introduced in eq. (2.3), we construct the associated basis {v n }: where we have made explicit the indices in the (D s − 4)-dimensional space. The basis {v n } is infinite dimensional for D s = 4 − 2 (because there are infinitely many independent terms of the form of eq. (2.3)), but at each order in the perturbative expansion only a finite number of basis elements contribute, as follows from inspecting the corresponding Feynman diagrams. We thus have In the HV scheme, the decomposition is independent of the number of external gluons. In particular, the value of n k can be determined from the amplitude with no external gluons, by examining the Feynman diagrams with the most singular gluons. These are ladder-type four-point diagrams with the gluons in the rungs. We find for instance that n 0 = 0, n 1 = 3 and n 2 = 5 for tree-level, one-and two-loop amplitudes, respectively. Our decomposition is similar to the one presented in ref. [28], but differs in the choice of basis tensors in eq. (2.11).
In practical calculations, one is interested in computing specific coefficients M (k) n in the decomposition of eq. (2.11). We construct the basis {v n } such that this operation is trivial, The calculation of the coefficients c n requires some technical operations on γ matrices that we present in appendix A. We then construct the dual basis {v n }, with elements Using the dual basis, we directly get Finally let us consider an amplitude with two identical quark pairs, which can be constructed by anti-symmetrizing the distinct-flavor amplitude M (k) over the two flavors [28,29]. It is then easy to see that the decomposition of eq. (2.11) requires an enlarged basis compared to the distinct-quark case of eq. (2.15). We thus define the tensors {ṽ n } as 20) and the decomposition of eq. (2.11) is over the sets {v n } and {ṽ n }. The basis tensors satisfy where the set {ṽ n } is constructed to be dual to {ṽ n } in the same way as in eq. (2.18). We finish this subsection with a comment on the case where D s is a finite integer D 0 s . All the discussion above holds, but one must be careful with a small detail. The basis of the Clifford algebra in eq. (2.3) now contains only a finite number of terms, and the basis of tensors {v n } is consequently restricted by the dimension D 0 s . If one wants to compute the coefficient of a given tensor v i , one must thus choose D 0 s large enough such that v i ∈ {v n }. Nevertheless, one can check that a calculation done in D 0 s dimensions agrees with the D s = D 0 s limit of the same calculation done in generic D s .

Two-loop Helicity Amplitudes for NNLO Phenomenology
We have established that helicity amplitudes in dimensional regularization are tensors in the (D s −4)-dimensional space and introduced a basis of that space on which we can decompose the amplitude. We should in principle compute all coefficients in the decomposition. However, it turns out that in a given phenomenological application not all coefficients may be relevant. We discuss below the two cases involving external quarks pertinent to the subject of this paper, the amplitudes with only external gluons being trivial in this regard.
Two-loop qqg . . . g amplitude: For the case of an amplitude with a pair qq of external quarks and any number of external gluons, there is a single coefficient to determine, see eq. (2.14). At order k in perturbation theory we call this object A (k) . It is computed using i.e. by tracing over the (D s − 4)-dimensional indices of the fermion line. In this paper we are mostly interested in k = 2.
Two-loop qqQQg . . . g amplitude: For a two-loop amplitude with two quark pairs of different flavors, qq and QQ, and any number of gluons there are in principle six coefficients to determine. However, in an NNLO computation (that is not loop-induced) the two-loop amplitude is interfered with the tree amplitude, which has a single tensor structure in the HV scheme. The contribution we must compute is of the form where we have used the orthogonality of the tensors v n and the fact that c 0 (D s ) = 1, see eq. (2.17). For NNLO corrections, it is thus sufficient to compute the coefficients M which amounts to computing the (D s − 4)-dimensional trace of M (2) on each fermion line. We define the amplitude A (k) (q,q, Q,Q, g, . . . , g) for any order k in an analogous way. This approach is similar to the one of ref. [28] and is in agreement with the prescription of ref. [36]. On a first look, it might however look inconsistent with the way qqQQ helicity amplitudes are defined in ref. [29]. Written in the formalism we have introduced in this section, the authors computeṽ and, given the relations of eq. (2.21), this would not necessarily give the same A (2) (q,q, Q,Q) defined in eq. (2.24). For phenomenological applications, however, one can show that only the so-called finite remainder is relevant [42], and we now show that the choices of eqs. (2.24) and (2.25) give the same result for this quantity. 3 We first recall that the infrared poles of a renormalized QCD amplitude M R have a universal structure, and we can write an amplitude in terms of its universal pole structure and a finite remainder which we will denote F [43][44][45][46]. More explicitly, for a two-loop amplitude we have where I (1) and I (2) are operators in color space. We refer the reader to appendix B for explicit expressions for these operators in the leading-color approximation of the amplitudes considered in this paper. Since 27) and the remainder computed from eq. (2.24) thus agrees with the one computed from eq. (2.25).
Finally, we now show that in the case of two pairs of identical quarks we can also use the definition of eq. (2.24) for NNLO phenomenology. The relevant contribution is the interference of the tree-level amplitude with the remainder, i.e.
where we denote with tildes the flavor exchanged objects. Here, we have used the orthogonality of the v n andṽ n up to O( ) to simplify the expression. Importantly, the right hand side of eq. (2.28) now only contains terms that can be computed through the definition of eq. (2.24).

Leading-Color Amplitudes
In this paper we compute a complete set of independent four-and five-parton helicity amplitudes in the leading-color approximation. More concretely, we keep the leading terms in the formal limit of a large number of colors N c , and scale the number of massless flavors N f whilst keeping the ratio N f /N c fixed. Each amplitude can be decomposed in terms of color structures whose coefficients are related by symmetry, and in this section we define our notation for the color decomposition of the amplitudes. We denote the fundamental generators of the SU (N c ) group by (T a ) i , where the adjoint index a runs over N c 2 −1 values and the (anti-) fundamental indices i andī run over N c values. We use the normalization Tr(T a T b ) = δ ab .
In this work, we will compute amplitudes where the external partons have well defined (either positive or negative) helicities, following the conventions of ref. [47]. We first discuss the four-point amplitudes. We will consider amplitudes for the scattering of four gluons, one quark pair and two gluons, and two distinct quark pairs. In the leading-color approximation we write where S n denotes all permutations of n indices and S n /Z n denotes all non-cyclic permutations of n indices. We write the particle type explicitly as a subscript, and all remaining properties of each particle (momentum, helicity, etc.) are implicit in the associated number. In the case of amplitudes involving quarks, we recall that the amplitudes A have been defined in eqs. (2.22) and (2.24). For the five-point case, we will consider the amplitudes for the scattering of five gluons, one quark pair and three gluons, and two distinct quark pairs and one gluon. In the leading-color approximation we write with similar notation as in the four-point case. For both the four-and five-point cases, the amplitude with two identical quark pairs can be obtained by anti-symmetrizing over the distinct flavors q and Q as discussed in the previous subsection. The kinematic coefficients of equations (2.29)-(2.34), denoted by the various A, are known as the leading-color partial amplitudes. They can be perturbatively expanded up to the two-loop order as where α 0 = g 2 0 /(4π) is the bare QCD coupling and A (k) denotes a k-loop partial amplitude. The partial amplitudes can be further organized in terms of the number of closed fermion loops, ranging from none up to the loop order, which each contribute one power of N f . We write (2.36) Figure 1. Representative Feynman diagrams for leading-color A (2) (g, g, g, g, g) amplitudes, contributing at order N 0 f , N 1 f and N 2 f .

Figure 2.
Representative Feynman diagrams for leading-color A (2) (q,q, g, g, g) amplitudes, contributing at order N 0 In the leading-color approximation, the structure of these amplitudes simplifies, receiving contributions only from planar diagrams. Representative diagrams for each of the fiveparton amplitudes we consider are given in figs. 1, 2 and 3.

Calculation of Planar Multi-Parton Amplitudes
In order to compute two-loop four-and five-parton amplitudes, we apply a variant of the unitarity method [10][11][12][13] suitable for automated numerical computations of multi-loop amplitudes [18][19][20]. The aim of the computation is to determine the coefficient functions c Γ,i and combine them with the master integrals I Γ,i in the standard decomposition of the amplitude: Here ∆ is the set of all diagrams that specify different propagator structures Γ in the amplitude. The index i runs over the set M Γ of master integrals associated with each propagator structure. In order to determine the coefficient functions c Γ,i , we promote eq. (3.1) to the integrand level. The integrand is denoted A( l ), where l represents the loop momenta, and we decompose it as [18] where P Γ is the set of propagators in the diagram Γ, and the ρ j denote inverse propagators. We extended the sum in eq. (3.1) to also run over surface terms contained in the set S Γ . These surface terms vanish upon integration but they are necessary to parametrize the integrand. The surface terms are constructed from a complete set of so-called unitaritycompatible integration-by-parts identities [18,[48][49][50]. For all the processes considered in this article we use the master/surface-term parametrization given in ref. [6], which only depends on the kinematics of the processes. While in amplitudes with fermions additional Lorentz-symmetry breaking terms may appear prior to integration (see e.g. [34,36] and the discussion below eq. (2.6)), they do not in our definition of helicity amplitudes in eq. (2.19).
The cancellation of these terms will be discussed in section 3.1.
In the numerical unitarity method, the coefficients c Γ,i in the ansatz (3.2) can be determined by building systems of linear equations through sampling of on-shell values of the loop momenta l . In the on-shell limit the leading contributions of eq. (3.2) factorize, where we label the set of the tree amplitudes associated to the vertices in the diagram Γ by T Γ , and the sum on the left-hand side represents the sum over all internal states on the internal edges of the diagram Γ. In eq. (3.3) the loop momenta Γ l is such that all propagators in P Γ are on-shell, and so in these limits we also probe diagrams Γ such that P Γ ⊆ P Γ (a relation that we denote as Γ ≥ Γ). Beyond one loop there exist diagrams in ∆ with doubled propagators. The numerators of such diagrams correspond to leading and subleading terms in their on-shell limits, and for the latter no factorization of the integrand into tree amplitudes is known. Nevertheless, as shown in ref. [19], one can systematically organize the set of cut equations (3.3) in such a way that all master-integral coefficients necessary to obtain the full amplitude can be computed.
In the following, we discuss the details of the procedure when applied to processes with fermionic degrees of freedom. First, in section 3.1, we discuss an approach that allows the use of finite fields in the presence of fermions. Next, in section 3.2, we discuss the implementation of the products of tree amplitudes with fermions. Finally, in section 3.3, we describe how these components come together to compute the integrated amplitude.

Finite Fields and Spinors
The extension of unitarity approaches to employ only operations defined in an algebraic field was proposed in ref. [38]. A finite-field based calculation allows to compute exact values for the integral coefficients c Γ,i of eq. (3.1) in a numerical framework. This idea was applied recently in [5,6] for pure gluon-scattering amplitudes, and here we discuss our implementation for amplitude computations with fermions.
From here on we denote by F an arbitrary number field. In practice, we will be interested in F being the field of rational numbers Q or the finite field Z p of all integers modulo a prime number p. In general, polynomial equations do not have solutions in F. This is at odds with the fact that in a unitarity-based approach one needs to generate loop momenta which satisfy a set of quadratic conditions corresponding to setting propagators to zero. In ref. [6], this was resolved by making sure that all scalar products between the momenta in the problem were F-valued. In the presence of fermions, the situation becomes more complicated due to the extension of the Clifford algebra beyond four dimensions. More specifically, terms such as µ γ [Ds]µ exhibit the (D − 4)-dimensional components of the loop momenta, which are in general not F-valued for on-shell momenta (more concretely, if we work on the field of rational numbers these components are in general irrational), leading to terms in the sub-currents of the Berends-Giele recursion that are not F-valued. To address this issue, we start with a parametrization of the on-shell spaces as in ref. [6] but always use normalized basis vectors. We write the two-loop momenta as where we denote their (D − 4)-dimensional components as µ 1 and µ 2 . Next, we choose an orthonormal basis n i of the (D − 4)-dimensional space with n 1 in the direction of µ 1 and write µ 1 = r 1 n 1 , µ 2 = µ 12 µ 11 r 1 n 1 + r 2 n 2 where r 1 = √ µ 11 , r 2 = µ 22 − µ 2 12 /µ 11 , (3.5) with µ ij = µ i · µ j . In a theory containing only vector particles we only ever need the values r 2 i , which are F-valued both on-and off-shell [6]. In contrast, in a theory with fermions, components of Berends-Giele currents will take the generic form a 00 + a 10 r 1 + a 01 r 2 + a 11 r 1 r 2 , which is not F-valued. In order to nevertheless be able to work in the field F, we consider the algebra V over the field F, with V the vector space spanned by the basis {r 0 = 1, r 1 , r 2 , r 1 r 2 } and equipped with the standard addition and multiplication. All components of the Berends-Giele are elements in the algebra, and can thus be written as a linear combination of the r i with F-valued coefficients. More concretely, this means we only need to determine the a ij in eq. (3.6) which are F-valued by construction. An important observation is that, although the coefficients a 10 , a 01 and a 11 in eq. (3.6) are non-zero in intermediate stages of the calculations, they vanish for the integrands of helicity amplitudes as defined in eq. (2.11). This cancellation of the r i terms holds in the HV scheme and is due to the projection onto the invariant tensors v n of eq. (2.11), which yields polynomials in the Lorentz invariants µ ij at the integrand level. To see this point more explicitly, consider the integrand M k ( l ) of an amplitude with an arbitrary number of quark lines where the subscript k encodes the dependence on the (D s − 4) spinor indices. 4 We can write the integrand in the form with the tensors f ρ 1 ···ρn,σ 1 ···σm with an invariant tensor are carried by metric tensors only. Hence, after contraction with invariant tensors, the integrand only depends on µ ij . In contrast, evaluating amplitudes that introduce a reference axis in the (D s − 4)-dimensional space would lead in general to a dependence on the components of the µ i and thus on the r i terms. This is the case for instance when considering gluonpolarization components in the (D s − 4) dimensions (as required in the CDR scheme) or generic values of the (D s − 4)-dimensional spinors η i , as written explicitly in eq. (2.6).
We finish with a comment that is not related to the use of finite fields but follows from the discussion above. Since our representation of fermion amplitudes is manifestly Lorentz invariant in (D s −4) dimensions prior to loop integration, the integrands of fermion amplitudes can be decomposed in terms of the same set of master integrands and surface terms as those used for amplitudes with gluons only [6,20].

Tree Amplitudes
In order to numerically calculate the necessary products of tree amplitudes used in the cut equations (3.3), we implement a Berends-Giele recursion [51]. The presence of the fermionic degrees of freedom means that we require concrete representations of the Clifford algebra in D s dimensions, where D s is even. It can be shown that integrands of the HV amplitudes defined in eqs. (2.22) and (2.24) depend at most quadratically on the parameter D s . As we must also take D s ≥ 6, we implement the recursion for three values of D s , specifically 6, 8 and 10. Explicit constructions can be found (for example) in [40,41] or obtained using the factorized definition in (2.2). Importantly, to obtain manifestly real representations of the Clifford algebra, we continue components of momenta to imaginary values keeping kinematic invariants real valued. For gluon amplitudes, the analytic continuation can be equivalently interpreted as changing the metric signature to g [Ds] = diag{+1, −1, +1, ..., −1}. As far as spinor representations are concerned, the two perspective are not equivalent as the latter would also alter the inner product of the spinors. Effectively we work in the alternating signature while maintaining the conjugation operation for spinors as defined in Minkowski signature.
In order to implement the prescription of eqs. (2.22) and (2.24) for computing amplitudes with external fermions, we first construct four-dimensional states with a specific helicity from Weyl spinors using the conventions of ref. [47]. To handle the (D s − 4)dimensional Clifford algebra we work with a canonical basis for the associated spinors η i κ = δ i κ and fixη κ i to be its dual as in eq. (2.5). Through eq. (2.4) we then construct the full set of D s -dimensional states associated with a given four-dimensional state. The projections in eqs. (2.22) and (2.24) then amount to the evaluation of (normalized) traces over the (D s − 4)-dimensional indices.
We close with two technical remarks. First, within our implementation all internal Lorentz indices are taken to be D s -dimensional despite the HV prescription that this should only be the case for one-particle-irreducible diagrams. This is allowed because the difference between this prescription and the HV prescription does not contribute to the helicity amplitudes as defined in eqs. (2.22) and (2.24). Second, with an appropriate normalization of the spinor states and their conjugates, the components of the spinors in internal state-sums also take the form of eq. (3.6), so no special treatment is needed in finite-field computations.

Amplitude Evaluation
We start by constructing the set ∆ of all propagator structures which are associated to a given amplitude in the decomposition of eq. (3.2). For this task we produce all cut diagrams in the full-color process employing QGRAF [52], followed by a color decomposition performed in Mathematica according to ref. [53]. In the latter step, tree-level decompositions for processes involving several fermion lines are necessary and we perform them following ref. [54]. We then take the leading-color limit and extract a hierarchically-organized set of propagator structures associated to the color-ordered amplitudes in eqs. (2.29)-(2.34). This decomposition is then processed by a C++ code. The master/surface-term decomposition which we employ is the same as the one used in refs. [6,20]. It was constructed using the computational algebraic geometry package SINGULAR [55] to solve syzygy equations that allow to obtain a set of unitarity-compatible surface terms.
Solving the multiple systems of linear equations associated to all cut equations (3.3) is achieved through PLU factorization and back substitution. To reconstruct the dependence of the master-integral coefficients on the dimensional regulators D and D s we sample over enough values to resolve their rational or polynomial dependence, respectively. In a generalization of ref. [32], the quadratic D s dependence is reconstructed from the evaluations at D s = 6, 8 and 10. Explicit D dependence on the integral coefficients is induced by the D-dependent surface terms. We sample multiple values of D randomly to extract the rational dependence of all master coefficients by using Thiele's formula [38,56].
The above approach to obtaining the coefficients in the decomposition of the amplitude in eq. (3.1) is implemented in a numerical framework which allows two independent computations. The first involves evaluation over the finite fields provided by Givaro [57]. We use cardinalities of order 2 30 and, to improve on the multiplication speed, we imple-ment Barrett reduction [58,59]. For a given rational kinematic point, we perform the computation in a sufficient (phase-space-point dependent) number of finite fields to apply a rational-reconstruction algorithm after using the Chinese Remainder theorem. The second mode of operation carries out the evaluations in high-precision floating-point arithmetic. In this case we do neither employ the technology to control algebraic terms described in section 3.1 nor do we use the refined computational setup based on real-valued operations.
The coefficients are then combined with master integrals as in eq. (3.1) to give the integrated amplitudes. For four-parton amplitudes we use the same implementation of the integrals as the one used in ref. [20] and for five-parton amplitudes we use the same as in ref. [6]. In the former case, we used our own calculation of a set of master integrals. In the latter case, we used the integrals of ref. [30] for the five-point master integrals, the integrals of ref. [60] for the lower point integrals, and our own calculation of the one-loopfactorizable integrals. In all cases, the polylogarithms in the -expansion of the master integrals are evaluated with GiNaC [61], which can be tuned to the desired precision.

Numerical Results for Helicity Amplitudes
In this section we present numerical values for leading-color two-loop multi-parton helicity amplitudes. We first present our results for four-parton amplitudes. These are known in analytic form [25,[27][28][29] and we use them as a validation of our approach. Then we present our new computation of five-parton helicity amplitudes. We include all N f corrections corresponding to closed massless-quark loops.
For each different choice of external partons we consider, we will show tables of numerical results for a full set of independent helicity assignments corresponding to a single partial amplitude in the color decompositions of eqs. (2.29)-(2.34). Furthermore, we present results only for distinct flavor configurations: as discussed in section 2, see eq. (2.28), results for finite remainders of amplitudes with identical quarks can be obtained by antisymmetrizing on the flavor assignments. In appendix B we give all the ingredients required for computing these remainders, in particular results for one-loop amplitudes expanded through order 2 .

Four-parton Amplitudes
We evaluate the four-gluon, two-quark two-gluon and four-quark amplitudes at the phasespace point 5 , A (2) 37639897 Table 1. The bare two-loop four-parton helicity amplitudes evaluated at the phase space point in eq. (4.1). We set the normalization factor A (norm) to A (1)[N 0 f ] ( = 0) for the amplitudes with vanishing trees, and to A (0) otherwise.
In table 1 we show numerical results for the bare two-loop four-parton helicity amplitudes. In order to expose the pole structure of the amplitudes (see appendix B) we normalize them to the corresponding tree-level amplitude if it is nonvanishing, or to the corresponding A (1)[N 0 f ] ( = 0) amplitude otherwise. The results have been obtained with exact values for the integral coefficients and with the master integrals evaluated to a precision that allows to show 10 significant digits.
We have validated our results by carrying out a set of checks. We verified that they satisfy the expected infrared pole structure [43]. We summarize the relevant formulae for this check in appendix B. Furthermore, we have carried out a systematic validation of the 0 contributions of all of our results against their known analytic expressions from refs. [25,[27][28][29]. For the four-gluon amplitudes we have compared directly the 0 pieces of our results with the analytic expressions of [25]. For the two-quark two-gluon and fourquark amplitudes we have used the one-loop results given in appendix B.3 to compute the corresponding finite remainders F (2) as defined in eq. (2.26). After accounting for the different choices of normalization for the H [n] ( ) operators (see appendix B) made in refs. [27] and [28], we have found perfect agreement.

Five-parton Amplitudes
We present results for the five-parton amplitudes evaluated at the phase-space point  We set the regularization scale µ to 1 and the normalization of the results is fixed by the expansion in eqs. (2.35) and (2.36). All results have been computed in the HV scheme. Given that our integral coefficients are computed as exact rational numbers, the final precision of our results is determined by how many digits we require from GiNaC [61] in the evaluation of the polylogarithms in the master integrals we use [30,60]. In table 2 we present results with 10 significant digits. All results in table 2 have been checked to satisfy the pole structure of two-loop amplitudes [43]. The one-loop amplitudes required for these checks have been obtained from our own setup, and cross-checked up to order 0 with BlackHat [17]. We present their numerical values in appendix B. The N 0 f piece of the all-plus five-gluon amplitude have been checked to reproduce the analytic result of [3], and for the other helicity configurations we have validated the results of [5] with our implementation. We also find agreement with the numerical results of the N 0 f terms of the two-quark three-gluon and four-quark one-gluon two-loop amplitudes which have been presented in the revised version of ref. [7]. Finally, we also cross-checked the pole structures of other helicity configurations, not explicitly shown. As our setup is a numerical one, this amounts to internal consistency checks of our computational framework.

Conclusion
We have presented the calculation of the planar two-loop four-and five-parton helicity amplitudes, extending the numerical variant of the two-loop unitarity method already used in [6,20] to amplitudes with fermions. Our results include all corrections associated with closed massless fermions loops. Numerical results for some of the amplitudes we have computed have been presented recently [7]. Given our results, the complete set of twoloop amplitudes required for a NNLO QCD calculation of three-jet production at hadron colliders in the leading-color approximation are now available.
We first described a formalism for computing multi-loop helicity amplitudes with external fermions in dimensional regularization, consistent with the approaches of refs. [26,28,36]. This was achieved by embedding the four-dimensional external fermionic states in D s dimensions and preserving the invariance of the amplitude under Lorentz transformations in the (D s − 4)-dimensional space. Within this formalism, we precisely stated our definition of helicity amplitudes and devised a numerical method to compute parton scattering amplitudes in the HV scheme. After interference with the Born amplitudes, changing to other regularization schemes can be achieved by known transition rules [39].
Our computational approach relies on a parametrization of the two-loop four-and five-point massless integrand in terms of master integrands and surface terms. With our definition of helicity amplitudes, we can reuse the same parametrization already used for four-and five-point gluon amplitudes independently of the type of partons. We extended the finite-field implementation of ref. [6] to fermion amplitudes, allowing us to compute exact master-integral coefficients for all partonic subprocesses at rational phase-space points. The computations were also performed in an alternative setup using floating-point arithmetic and we find agreement between the two variants of our numerical method.
We present reference values for helicity amplitudes. These are obtained by combining the master-integral coefficients we compute with the corresponding master integrals, in particular using the recently obtained analytic expressions for five-point integrals [30,31]. We have validated our results in a number of ways: we reproduce the results for two-loop four-parton helicity amplitudes computed from their known analytic expressions [25,27,28], we find the correct infrared structure of each amplitude and we validate the finite pieces of recently published five-parton results [1][2][3][4][5][6][7] as detailed in section 4.2.
The techniques developed in this paper show the potential for the automation of twoloop multi-particle amplitude calculations in the Standard Model. Our numerical approach is relatively insensitive to the addition of scales. Having already implemented vector and spinor fields, we are now ready to explore processes of phenomenological relevance that include jets, (massive) gauge bosons and leptons in the final state. While techniques for computing two-loop master integrals progress and new methods appear for handling infrared divergent terms in real-real and real-virtual contributions, we expect to provide a program that can deliver one-and two-loop matrix elements necessary for computing precise QCD predictions for the LHC.

Acknowledgments
We thank J. Dormans and M. Zeng for many useful discussions. We thank Z. Bern and L.J. Dixon for providing analytic expressions from references [25,27], which we employed to validate our results for the two-loop four-gluon and two-quark two-gluon amplitudes. We thank S. Badger, C. Brønnum-Hansen, H. B. Hartanto and T. Peraro for correspondence concerning the numerical results presented in ref. [7]. We thank C. Duhr for the use of

A Operations on γ Matrices
In the following we derive the values of the contraction of the tensors w 0 and v n which are used in eqs. (2.14) and (2.16). In this appendix we take d to be an even integer denoting the dimension of the space for which the Clifford algebra is defined and we denote the dimension of the γ-matrix representation by Since amplitude computations are homogeneous in the factor d t , it can be factored out and replaced by a suitable value in order to suit the four-dimensional limit. In this appendix, we keep the parameter d t in analytic form in order to maintain a consistent finite-dimensional algebra and for clarity of the equations. In the main text we use formulas with the replacement d t → 1 imposed, which is the value consistent with a calculation in dimensional regularization [41]. We start with the trivial case of w 0 which appears in eq. (2.14). It is easy to find that For the tensors v n of eq. (2.16) we must first consider traces of γ-matrix chains of the form .
Unitary representations for the Clifford algebra can always be found as explained for example in ref. [40]. We will require the following traces of antisymmetric γ-matrix chains, where the summation runs over all permutations S n of n elements. The traces are computed in fixed integer dimensions where the dimensions of the γ-matrix representation is taken to be d t dimensional. For contracted Lorentz indices we will also use that The sum counts the number of antisymmetric tensors of rank n in d dimensions, which is the number of ways to choose an ordered subset of n elements from a fixed set of d elements.
With these preparatory equations we can compute the inner products of the v n tensors of eq. (2.15) which yield the normalisation factors c n of eq. (2.17): (A.7) In the above formulas the summation over the indices ν i is trivially performed. In the next step, we isolate a contribution of the form that we computed in eq. (A.6), which gives the same result for each permutation σ but multiplied by sgn(σ). The final results follows trivially. As expected, for each n the result has zeros in the dimensions d for which there are insufficient distinct labels µ i and ν i available to form antisymmetric index configurations of n indices. We recall that in the main text we set d t = 1 and d = D s − 4.
Finally we collect the results for the contractions of the tensorṽ m and v n required for amplitudes with two quark lines of identical flavor, see eq. (2.20). These contractions lead to a single trace instead of a product of traces as was the case in eq. (A.7). We refer to the above intuitive argument: tensor contractions including v n orṽ n vanish whenever the dimensionality d is insufficient to accommodate the respective antisymmetric index arrangements in the Lorentz indices µ i and ν i . In particular this implies that contractions including the tensors v n are proportional to d for n = 0. We find that Here p mn (d) is a polynomial-valued matrix which we will not require explicitly for the present paper, and in the last equality we made explicit the fact that in this paper we are interested in the case d = D s − 4 = O( ).

B Divergence Structure of Two-Loop Five-Parton Amplitudes
We use the HV dimensional regularization scheme to handle both ultraviolet and infrared divergences. UV divergences are removed through renormalization and the remaining infrared poles can be computed from the corresponding lower-order amplitudes [43][44][45][46]. In this appendix we detail this procedure. Reproducing the pole structure of the amplitudes we have computed is an important check of our results.

B.1 Renormalization
We perform renormalization of the QCD coupling in the MS scheme. It is implemented by replacing the bare coupling by the renormalized one, denoted α s , in eq. (2.35). The bare and renormalized couplings are related through where S = (4π) e − γ E , with γ E = −Γ (1) the Euler-Mascheroni constant. µ 2 0 is the scale introduced in dimensional regularization to keep the coupling dimensionless in the QCD Lagrangian, and µ 2 is the renormalization scale. In the following, we set µ 2 0 = µ 2 = 1. The leading-color coefficients of the QCD β-function are The perturbative expansion of the renormalized amplitude is where λ is the power of g 0 in the tree amplitude, with α 0 = g 2 0 /(4π) and similarly for α s . For four-parton amplitudes λ = 2, and for five-parton amplitudes λ = 3. The renormalized amplitudes A (i) R are related to the bare amplitudes A (i) as follows:

B.2 Infrared Behavior
The poles of renormalized amplitudes are of infrared origin and can be predicted from the previous orders in the perturbative expansion [43][44][45][46]: [n] ( )A [n] and I (2) [n] depending on the number and the type of the scattering particles. This dependence is denoted by the subscript [n]. For amplitudes in the leadingcolor approximation and for which all quark lines have distinct flavor, the operators I (1) [n] and I (2) [n] are diagonal in color space and can be written in a very compact form. The operator I (1) [n] is given by with the indices defined cyclically. The index a i denotes a type of particle with momentum p i , i.e., in the context of our paper, a i ∈ {g, q,q, Q,Q}. We introduced the auxiliary symbols γ a,b , symmetric under the exchange of indices, γ a,b = γ b,a , and defined according to: γ g,q = γ g,q = γ g,Q = γ g,Q = γ g,g + γ q,Q 2 , γ q,q = γ Q,Q = 0 . (B.7) The operator I and H [n] ( ) is a diagonal operator at leading color that depends on the number of external quarks and gluons in the process, The poles of the bare amplitudes, as presented for example in tables 1 and 2, can be recovered from those of the renormalized amplitude by using eqs. (B.4).

B.3 Numerical Results for One-loop Amplitudes
To predict the expected pole structure of the amplitudes computed in section 4 it is necessary to compute corresponding one-loop results up to high enough order in . For completeness, we present one-loop results in tables 3 and 4 which we have obtained with our own implementation of one-loop numerical unitarity. The expansion has been performed up to O( 2 ) in order to allow the evaluation of finite remainders as in eq. (2.26). This was used to reproduce the analytic results for finite remainders of the qqgg and qqQQ amplitudes of refs. [27,28]. The results are normalized to remove overall phase ambiguities in the amplitudes, choosing the tree-level amplitude if it does not vanish, or the leading term of the one-loop amplitude otherwise. We present numerical values with 10 significant digits. (1 + q , 2 − q , 3 + g , 4 + g ) 0 0 −0.1818181818 −1.074210422 −3.518712119  Table 3. The bare one-loop four-parton helicity amplitudes evaluated at the phase space point in eq. (4.1). We set the normalization factor A (norm) to A (1)[N 0 f ] ( = 0) for the amplitudes with vanishing trees, and to A (0) otherwise. ( Table 4. The bare one-loop five-parton helicity amplitudes evaluated at the phase space point in eq. (4.2). We set the normalization factor A (norm) to A (1)[N 0 f ] ( = 0) for the amplitudes with vanishing trees, and to A (0) otherwise.