Planar Two-Loop Five-Gluon Amplitudes from Numerical Unitarity

We present a calculation of the planar two-loop five-gluon amplitudes. The amplitudes are obtained in a variant of the generalized unitarity approach suitable for numerical computations, which we extend for use with finite field arithmetics. Employing a new method for the generation of unitarity-compatible integration-by-parts identities, all helicity amplitudes are reduced to a linear combination of master integrals for the first time. The approach allows us to compute exact values for the integral coefficients at rational phase-space points. All required master integrals are known analytically, and we obtain arbitrary-precision values for the amplitudes.


I. INTRODUCTION
The expected precision of upcoming cross-section measurements at the Large Hadron Collider at CERN currently drives the development of new computational methods in perturbative quantum-field theory, and in particular in QCD. At the same time, formal advances in the field of amplitudes are finding more and more traction and are helping to devise more efficient approaches to multi-loop computations. Here, we aim to contribute to these developments by refining and extending methods to compute phenomenologically relevant two-loop scattering amplitudes. With our results, we advance the state of the art in the calculation of loop amplitudes with many external particles. Besides their importance for phenomenology, the new results highlight the potential of the methods we employ.
In this paper, we present the computation of two-loop gluonic amplitudes which contribute to three-jet production at hadron colliders at next-to-next-to-leading order (NNLO).
Predictions for these processes can be used for constraining the strong coupling, as it can be extracted from precision measurements of three-to two-jet-production ratios (see e.g. ref. [1,2]). While the finite parts of our results can be integrated over phase space, significant additional developments will be required to obtain the full NNLO predictions.
Probably the most challenging step will be the consistent combination of all virtual and real-radiation contributions to NNLO accuracy. Nevertheless, we are optimistic that by providing key ingredients towards this goal, such as the two-loop scattering amplitudes, we will spur progress in the field to obtain NNLO predictions for three-jet production in the coming years.
Currently, only a limited set of amplitudes with five or more external particles are known to two-loop order in QCD. First benchmark results for two-loop five-gluon amplitudes were obtained for the all-plus-helicity configuration [3,4], whose computation relied on compact analytic expressions for the integrands. These were computed with the unitarity method [5][6][7][8], coupled with efficient integrand-reduction techniques [9][10][11][12][13]. The final amplitudes were first obtained by integrating the optimized integrands numerically and later by explicit reduction to master integrals [14]. The relevant master integrals have been independently evaluated in ref. [15]. More recently, the all-plus two-loop amplitudes were obtained from a recursive approach [16] that was also extended to six and seven external gluons [17,18].
During the final stages of this work, benchmark results for planar five-gluon amplitudes with generic helicity assignments were presented [19]. The amplitudes were integrated with a combination of numeric and analytic techniques, after partial reduction using unitarity cuts.
In this paper, we apply the numerical unitarity approach [9,[20][21][22], which was recently extended to two-loops [23][24][25], to compute the full set of planar five-gluon two-loop amplitudes. This computation marks the first time that a full set of two-loop five-scale amplitudes is reduced to a linear combination of integral coefficients multiplied by corresponding master integrals. The set of master integrals is minimal, as all possible integral relations have been imposed. The potential of this approach has recently been demonstrated in the context of planar four-gluon two-loop amplitudes [25], and in order to handle the complexity of the five-point amplitudes we refine it in the present work in two ways.
Firstly, we use an improved approach to obtain analytic expressions for integration-byparts (IBP) identities which are important in our construction of the amplitude. The unitarity method [5,6,8] matches unitarity cuts to an ansatz of the loop amplitudes written in terms of master integrals. When considered at the integrand level, one-loop unitarity approaches [9,21] relate cuts to an ansatz for the amplitude's integrand written in terms of surface terms (spurious numerators that integrate to zero) and master integrands (the integrands corresponding to the master integrals). The numerical unitarity method at the two-loop level [23][24][25] is a natural generalization of the one-loop integrand level approach: the unitarity cuts are matched to a complete parameterization of the integrand in terms of specially tailored IBP relations [26][27][28] (the surface terms) and master integrands. As such, the method directly achieves a reduction to master integrals at the integrand level.
Compared to more canonical approaches, the challenging inversion of the IBP systems required to obtain reduction tables for tensor integrals [27] is avoided. While we leave a more detailed discussion of the way we generate the required IBP relations to future work [29], we describe the central formal improvements for their derivation.
Secondly, we adapt numerical unitarity to exact numerical computations based on the finite field approach [30] and its application to unitarity cuts [31]. The finite field arithmetics allow us to compute exact values for the integral coefficients at rational phase-space points.
Combined with the analytic expressions for two-loop master integrals with five massless external legs [15], we obtain stable results whose precision only depends on the number of digits to which the multiple polylogarithms appearing in the master integrals are evaluated.
The extension to finite fields will be crucial for a future use of functional reconstruction techniques to determine analytic forms of the amplitudes from the exact numerical results.
The article is organized as follows. In Section II we summarize the numerical unitarity method, as well as the two refinements already highlighted. In Section III we give explicit details of our implementation. In Section IV we present our results for the four independent leading-color helicity amplitudes required for all five-gluon scattering processes at a given phase-space point. In Section V we give our conclusions and outlook. Finally, in Appendix A we discuss the universal infrared structure of the amplitudes which we use to validate our results.

II. NUMERICAL UNITARITY APPROACH
The main goal of this paper is the calculation of two-loop corrections to five-gluon scattering amplitudes in the leading-color approximation for any helicity configuration in pure Yang-Mills theory. We write the perturbative expansion of a bare five-gluon helicity amplitude as where α 0 = g 2 0 /(4π) is the bare QCD coupling. The set σ denotes all non-cyclic permutations of five indices, and we make it explicit that the color structures beyond tree-level give only a factor of N C at each order in the expansion in the leading-color approximation. The planar amplitudes A (j) are functions of the momenta p σ(i) , helicities h σ(i) and color indices a σ(i) . In the leading-color approximation there is a single color structure and it is thus sufficient to specify an helicity assignment for the ordered set of legs, In this section we describe our approach to the calculation of the A (i) . Although in this paper we focus on the leading-color contributions, the approach we use is completely generic and applicable beyond this approximation.

A. Overview
We apply a variant of the unitarity method [5][6][7][8] suitable for automated numerical computations of multi-loop amplitudes [23][24][25]. Our approach starts from the standard decomposition of the amplitude in terms of master integrals, where the index Γ labels the different propagator structures and i denotes an additional degeneracy for cases where multiple master integrals with identical propagator structures appear, as can happen beyond one loop. The set ∆ denotes the collection of diagrams which specify the possible propagator structures of the amplitude, and only depends on the kinematics of the amplitude (i.e. on the momenta and masses of the particles involved).
In figs. 2 and 3, we give an explicit example for a set of propagator structures ∆ and the corresponding set of master integrals together with their multiplicity, in the case of twoloop five-gluon leading-color amplitudes. The aim of the computation is to determine the coefficient functions c Γ,i which contain the process specific information.
Next, we promote the decomposition of the amplitude of eq. (II.3) to the integrand level.
The integrand, which we denote A( l ) with l representing the loop momenta 1 , . . . , L in a L-loop calculation, is decomposed as [23] where the ρ j are the inverse propagators and P Γ the set of propagators associated with the topology Γ. Compared with (II.3), we have extended the sum on i to also go over so-called surface terms in the set S Γ , which vanish upon integration but are required to parameterize the integrand. We work in dimensional regularization with D space-time dimensions, and thus allow D-dependent surface terms. All in all, the decomposition (II.4) is rather universal, depending only on the set of diagrams ∆ and the allowed power counting of the theory under consideration. As in the decomposition of eq. (II.3), the spectrum-specific information is carried by the coefficients c Γ,i . Obtaining a complete parameterization of the integrand of the amplitude is a non-trivial process, and in the next subsection we explain how this was achieved. For now, it is sufficient to assume that such a decomposition exists.
The ansatz (II.4) holds as a function of the loop momenta l . Provided one can evaluate A (k) ( l ) at specific values of l , it can be used to construct a system of equations for the coefficients c Γ,i . In a generalized unitarity calculation, the coefficients c Γ,i are determined from a set of cut equations obtained by setting l to specific on-shell configurations. Indeed, the leading terms in the various on-shell limits of the ansatz (II.4) factorize, yielding the cut . (II.5) Here, the set T Γ labels all tree amplitudes corresponding to the vertices in the diagram Γ, and the state sum runs over all internal states allowed by the theory. The Γ l correspond to a configuration of the loop momenta where the propagators in the set P Γ associated with the topology Γ are on-shell. Although this limit probes all topologies Γ for which P Γ ⊆ P Γ , through a top-down approach, i.e. starting from the topologies with the most propagators, we can make sure that the only unknowns in eq. (II.5) are the coefficients c Γ,i belonging to Γ itself. By sampling over enough values of Γ l , we can build a system of equations big enough to determine all the coefficients. We note that beyond one loop there are topologies that correspond to subleading terms in the on-shell limit Γ l , in which case no factorization of the ansatz (II.4) is known. This is a minor obstacle that can be easily overcome [24].
For the process we are concerned with in this paper, two-loop gluonic amplitudes in the leading-color approximation, the state sum in eq. (II.5) goes over the (D s − 2) gluon helicity states. Furthermore, as the surface terms are functions of the space-time dimension parameter D, the coefficients c Γ,i obtained by solving the cut equations are functions of both D s and D. The value chosen for D s defines different regularization schemes. We have implemented both the 't Hooft-Veltman (HV) scheme [32] where D s = D = 4 − 2 and the four-dimensional helicity (FDH) scheme [33,34] where D s = 4.

B. Efficient Construction of Unitarity-Compatible IBP Identities
We now return to the discussion of how the decomposition in eq. (II.4) is achieved. Here we present the generic framework of our method but leave more specific details for a future dedicated publication [29]. The surface terms are constructed from a complete set of socalled unitarity-compatible integration-by-parts (IBP) identities [23,[35][36][37]. These relations are built such that they do not involve any new propagator structures besides those in ∆, the full set of propagator structures we started with. As is the case in more canonical approaches to multi-loop amplitude calculations, our approach requires a complete set of IBP identities for each of the propagator structures [25]. In the context of the ansatz in eq. (II.4), completeness means that surface terms and master integrands span the full function space prescribed by the power counting of the theory. The system of equations in eq. (II.5) is then invertible and by solving it we achieve a full reduction of the amplitude to master integrals.
In principle, the desired decomposition can be achieved by IBP reducing all allowed insertions of irreducible scalar products (ISPs) 1 on a given propagator structure. The reduction identities may then be used as the set of surface terms in eq. (II.4). In practice, the currently available reduction programs [38][39][40][41][42][43] struggle with the high tensor-rank appearing in two-loop five-gluon amplitudes as the size of the IBP system that needs to be solved becomes prohibitive. Instead, we follow an alternative path to obtain unitarity-compatible relations which has recently received much attention [23,25,[35][36][37][44][45][46][47]. By only generating unitarity-compatible IBP relations from the start, we ensure that the size of the system remains manageable. In our approach, the solution of the IBP system is done by solving the cut equations (II.5), which allows to efficiently organize the calculation.
We start from the construction of the 'transverse IBP identities' from refs. [23,25], but further improve the construction of the remaining identities by solving syzygy module equations. These are closely related to the defining equation for the so-called IBP-generating where 1 ≤ k ≤ L is the loop momentum label, ν is the Lorentz index and there is no summation over j. We require f j to be polynomial in the dot products built from loop momenta and external momenta. The vectors {u ν k } are required to be polynomial vectors in the loop momenta. The above equation ensures that the IBP relation contains no integrals with raised propagator powers and is therefore a suitable unitaritycompatible IBP relation. In the context of eq. (II.4), the left-hand-side of (II.7) gives a valid surface term.
Our method uses inverse propagator variables to trivialize the separation of terms that vanish on-shell from ISPs that survive on the cut, as in [23,25]. However, we further reduce the polynomial degree of the syzygy equations by simply rewriting the defining equation for IBP-generating vectors (II.6) in terms of the components of the loop and the external momenta. If starting from a formulation in terms of so-called Baikov polynomials [37], this method of solving eq. (II.6) may be viewed as a variant of the intersection method [46]. We write down the following ansatz for the IBP-generating vector where we have an implicit sum over the label a for independent loop momenta and label c for independent external momenta. Consistent with the fact that the original vectors {u ν k } are polynomial we require the vector components u loop ka and u ext kc to be a polynomial in the dot products built from loop momenta and external momenta. Since the latter can be expressed in terms of ISPs and inverse propagators, the variables u loop ka , u ext kc and f j are polynomials in ISPs and inverse propagators. Eq. (II.6) then becomes . . .
where there is implicit summation over a, c, k, and ν, and the inverse propagator labels j(i) run over all propagators in the set P Γ , for the unknown polynomials u loop ka , u ext kc and f j . The equations are defined over the freely generated polynomial ring given by ISPs and inverse propagators. Syzygy equations can be solved using algorithms in computational algebraic geometry implemented in e.g. the Singular computer algebra system [48]. We will give more details of how we solved the syzygy equations using Singular in the next section.
Once a generating set of vector components u loop ka and u ext kc has been obtained, a sufficient set of IBP relations is obtained by multiplying the generators with irreducible numerators.
At this stage, we explicitly impose off-shell power-counting conditions (i.e. we consider both on-shell and propagator variables), sometimes after forming linear combinations of IBP relations. In a second step the independence of IBP identities has to be established. Given the manifest dependence on propagator variables and the vectors compatibility with unitarity cuts, it is natural to switch between on-shell and off-shell as well as numerical and analytic perspectives to simplify computational steps [23,37]. In particular, the validation of the linear independence of the off-shell relations [23,25] is most easily performed by setting propagator variables to their on-shell values, as redundancy of the relations is then manifest.
Finally, we note it is often convenient to count the number of master integrals on-shell [23,25,37,44].

C. Numerical Unitarity over an Arbitrary Number Field
Numerical computations have better scaling properties than analytic ones for processes with a large number of scales. However, the stability and efficiency of analytic results for amplitudes often offer major advantages over numerical evaluations, specially when used in conjunction with Monte Carlo programs. The boundary between the two approaches is becoming blurred in the field of scattering amplitude calculations with the recent introduction of functional reconstruction techniques [25,30,31]. In this approach, numerical evaluations are used to completely reconstruct analytic expressions from numerical samples.
Whilst already successful in the reconstruction of two-loop amplitudes with two scales [25], it is clear that multi-scale problems such as the five-gluon amplitudes are significantly more complicated. Not only must more efficient techniques be used for the functional reconstruction [31], but the loss of precision associated with floating-point arithmetics in multi-scale computations poses a considerable difficulty for the reconstruction procedure. To address this issue, in this section we reformulate numerical unitarity for gluonic amplitudes so that it only employs operations defined on a field (addition, subtraction, multiplication and di-vision). In this way, we open the door to the use of exact arithmetics, such as those of rational numbers or finite fields. This eliminates any question of numerical stability, whilst still leveraging the speed of modern computer hardware.
The extension of numerical unitarity to employ only field operations is non-trivial. Indeed, operations such as taking square roots or, more generally, solving generic polynomial equations are not allowed, and these feature heavily in many formulations of generalized unitarity. Here, we build on the ideas proposed in ref. [31] and use them in the context of multi-loop numerical unitarity. More concretely, there are two components of a typical numerical unitarity formulation that must be re-examined: (1) the generation of a set of four-momenta satisfying momentum conservation and on-shell conditions for the external kinematics, and (2) the parameterization of the internal loop momenta, in particular for the solution of the quadratic on-shell equations and the construction of specific sets of surface terms. In this section we will describe how these obstacles can indeed be overcome. Furthermore, this will be achieved without needing to introduce any notion of complex numbers.
In the following we work over an arbitrary field which we denote by F, and which for all practical purposes can be considered to be the rational numbers Q.

External Kinematics
The problem of generating a set of F-valued external momenta which satisfy on-shell conditions as well as momentum conservation can be solved in a number of ways. In this work, we choose to parameterize the external kinematics based upon so-called "momentum twistors" [49] as proposed in refs. [12,31], where the reader can find further details. In  .
(II. 12) implementation. We also note that a wisely chosen parameterization such as that of [31] means that the invariants and Gram-determinants are compact functions of the parameters.
This simplifies the functional dependence of the amplitude on these parameters, making them particularly suitable for a future reconstruction of the analytic expressions.

On-shell Momenta
As discussed around eq. (II.5), in order to compute integral coefficients in a unitarity approach we need to generate loop momenta which satisfy a set of conditions which set some propagators on-shell. For a two-loop calculation, this set of topology specific, kinematicallydependent quadratic conditions define an algebraic variety in the 6-dimensional space in which we embed the loop momenta. A direct approach for finding a rational parameterization of this variety, i.e., a parameterization in terms of a set of F-valued parameters that only uses the operations that are defined on a field, is non-trivial. Instead we take inspiration from the fact that the integrand is a rational function of irreducible scalar products (ISPs), and the µ ij variables which we shall define shortly. We will see that in a set of adapted coordinates it is trivial to generate loop momenta such that the ISPs and µ ij are F-valued. Therefore, the integrand evaluated on such a loop-momentum configuration will also be F-valued. We then represent the loop momenta in a phase-space dependent way, circumventing the rational parameterization required when using the standard 6-dimensional representation. We now give more details about this procedure.
We begin by parameterizing the loop momenta l (l = 1, 2, 3 see fig. 1) as [23,25]  where (G l ) ij is the inverse of the Gram matrix, The index set B p l labels the external momenta which leave the strand l. These momenta are completed with other independent external momenta p i , with i ∈ B t l , so as to span the whole scattering plane. This parameterization follows the conventions of ref. [25], with the caveat that the vectors spanning B ct are no longer normalized.
The inverse coordinate transformation is often useful and is given by The on-shell variety is then defined by setting the propagator variables ρ li to zero. In D-dimensions the variables α li form an independent complete set of coordinates on the variety, corresponding to polynomials in momentum variables. They are the irreducible scalar products we already mentioned previously.
By considering the α li as a set of independent coordinates and taking them as F-valued, the constraint equations (II.15) imply that the µ ij are also F-valued. We must now construct an explicit representation of the loop momenta. This is required, for instance, for the calculation of the tree amplitudes in eq. (II.5). In the 4-dimensional slice this is trivial to achieve with a standard cartesian basis. However, if we were to also do this for the (D − 4)-dimensional space we would generically be required to take square roots.
We now present our solution to this problem in the context of a two-loop calculation, but with µ ij = − µ i · µ j , here using the Euclidean scalar product. We stress that the basis vectors n 1 andñ 2 used to represent the (D − 4)-dimensional space are no longer of unit norm. For each on-shell point, which corresponds to a different value of the µ ij , this affects how we calculate the scalar product between two vectors w a and w b . Explicitly, Such a scalar product then allows F-valued representations of loop momenta which satisfy the on-shell conditions for an arbitrary topology.
This momentum representation is almost all that is required to use a Berends-Giele recursion [50] to calculate products of tree amplitudes. The remaining difficulty is in performing the sum over helicities at each cut line, as this requires explicit representations of the D s -dimensional polarization vectors. To remedy this, we choose to avoid constructing the polarization states by trading the helicity sum for the insertion of a light-cone projection operator (see e.g. [51]), where η is an arbitrary F-valued light-like reference vector satisfying η · l = 0 and η ·ñ = 0.
Note that in all but one cut line per loop, the Ward identity allows us to drop the second term of equation (II.24). For the remaining cut line of each loop, we re-express the projector as a sum of a direct product of vectors over F.
We note that as the irreducible scalar products of eq. (II.19) are expressed in terms of a set of non-normalized vectors, this affects the representation of so-called 'traceless completion' surface terms [25], i.e. surface terms associated with variables in B ct ∪B . This normalization now explicitly appears in the parameterization of the loop momenta, see eq. (II.13). For example, consider a traceless completion surface term associated to the transverse vector Only by including the factor of n 2 i is the numerator insertion (II.25) a surface term. We finish this section by noting that this procedure applies both for planar and nonplanar cases, and is easily generalized to higher (and lower) loops. In summary, using the steps described in Sections II C 1 and II C 2, all contributions are manifestly F-valued, and as an added benefit we never needed to introduce complex numbers.

III. IMPLEMENTATION FOR PLANAR FIVE-GLUON AMPLITUDES
We have implemented the techniques described in section II in a C++ code. In this section, we first discuss the implementation of the cut equations. Then we describe how these are solved in practice to compute the master integral coefficients and, finally, how we obtain the amplitude at a given kinematic point.

A. Construction of Cut Equations
We first construct the full set of propagator structures that appear in the problem (see figure 2). This is achieved by generating all cut diagrams for the full color process using QGRAF [52] and then color decomposing them according to [53] in Mathematica. By taking the leading-color limit of this decomposition and extracting the coefficient of a given trace, we build the set of propagator structures relevant for the color ordered amplitude, which we then organize hierarchically. This is then passed on to a C++ code. The next step in constructing the cut equations is the parameterization of the integrand of the amplitude in terms of master integrands and surface terms. In order to solve the syzygy equations described in section II B for each topology in fig. 2 we use the package Singular. For these computations we find that the 'slimgb' algorithm [54]  We note that a fast evaluation of the surface terms is required for an efficient implementation. To this end we express the surface terms as IBP-generating vector components multiplied by derivatives of the irreducible scalar products, as well as the total divergences of the vectors multiplied by irreducible scalar products. Furthermore, to improve both evaluation time and compilation time of the IBP vector components, we found it useful to employ the facilities provided by FORM [55] for the simplification of multivariate polynomials [56,57].
Given the parameterization of the integrand with surface terms, any linearly-independent tensor insertion can be used as a master integral. As such, it is trivial for us to change basis of master integrals, by filling the master space of each topology with any set of integrals independent under the IBP relations. This freedom was especially useful during testing, where we could make use of finite integrals to control the structure. In order to numerically calculate the appropriate products of tree amplitudes necessary to reconstruct the integrand in our algorithm, we implement a Berends-Giele recursion [50] tailored to directly compute multi-loop cuts using D-dimensional momenta and states living in D s dimensions. This is both a high-performance and flexible choice, as changing the field content requires only implementing new Feynman rules. The setup is particularly useful for our numerical computations in dimensional regularization, as it is straightforward to evaluate the products of tree amplitudes at different values of D s in order to reconstruct the functional dependence on the parameter, in a way that automatically recognizes if a given two-loop cut has a linear or quadratic dependence on D s . Caching for multiple D s values is built-in. External momenta are taken to live in 4 dimensions, as do the associated gluonic polarization states. The D-dimensional loop momenta are represented in 6-dimensions, as described in detail in section II C.

B. Amplitude Evaluation
Having constructed the cut equations (II.5) as described in the previous section, we now solve them for the integral coefficients at fixed values of the kinematics. Both equation (II.5) and its hierarchically subtracted analogue describe linear systems in the ansatz coefficients.
For fixed values of D and D s , we evaluate the equations numerically over enough randomly chosen on-shell momenta configurations to form a linear system that constrains the coefficients. We then solve this system for the coefficients using standard PLU factorization and back substitution.
In order to reconstruct the D and D s dependence of the coefficients we first sample D s over three distinct values [25] in a generalization of [21], allowing us to easily implement both the FDH and HV flavors of dimensional regularization. As D s is restricted to be greater than or equal to the embedding space of the loop momenta (6 for a two-loop calculation) we pick D s = 6, 7, 8. The values of D are chosen randomly. The D-dependence is obtained by repeating the computation for sufficient (a priori unknown) values of D from which the rational dependence is reconstructed [30,31] using Thiele's formula [31,58]. After evaluation at a single phase-space point, the exact denominator as well as rank of the numerator function can be stored, allowing the use of simple polynomial inversion techniques and less sampling.
In practice, implementing our approach over the rational numbers Q will result in a slow calculation. However, the final master integral coefficients are strongly constrained by physical properties and it is expected that their resulting form will be compact. We thus follow the approach outlined in [30,31], using finite fields F p . We implement the algorithm over the finite fields provided by Givaro [59]. We use various cardinalities of order 2 30 , implementing Barrett reduction [60,61] in order to improve the speed of finite field multiplication. For a given kinematic point in Q, we perform the computation in a sufficient number of finite fields to apply a rational reconstruction algorithm, also provided by Givaro.
In order to obtain an expansion for the amplitudes, we combine the coefficients with expansions of the appropriate master integrals. We list the topologies with master integrals in figure 3. For the five-point master integrals we choose the basis of [15] in order to make use of the publicly available implementation distributed with the paper. For lower point integrals, we implemented the analytic expressions provided in [62]. In the case of factorizable topologies we choose the scalar integral as a master and calculate them independently. In order to evaluate the necessary multiple polylogarithms we use the implementation found in GiNaC [63], which can be tuned to the desired precision. All integrals have been independently numerically validated with Fiesta 4 [64]. We note that the choice [62] of an integral with a squared central propagator for the second slashed-box master is important in order to avoid spurious weight-five contributions in intermediate stages.  [15,62], the precision of the results we present is only limited by the number of digits we ask from GiNaC [63] when evaluating the master integrals, which can be arbitrarily increased.
(1 + , 2 + , 3 + , 4 + , 5 + ) -5.000000000 -3.8931790255 5.9810885816 (1 − , 2 + , 3 + , 4 + , 5 + ) -5.000000000 -16.322002103 -10.383813287  To validate our results, we first reproduce the universal ultraviolet/infrared pole structure of the amplitudes [65], which we summarize in appendix A in the context of two-loop fivegluon amplitudes in the leading-color approximation. Computing the prediction for the pole structure requires the five-point one-loop amplitudes up to order for all four independent helicity configurations. For this, we used our own implementation of numerical unitarity at one loop, which we checked up to order 0 against results from BlackHat [22]. Within our check, the precision of the numerical results we obtain for the poles is limited by the fact that we use Fiesta 4 [64] for the one-loop pentagon integral. Using our code, we confirm the published results for the all-plus helicity amplitudes [3,14,16]. We also validate the results of [19] which appeared during the final stages of the preparation of this article.
Finally, we note that although we only present results for the four independent helicity configurations, we have also verified the pole structure for the other helicity configurations obtained by parity conjugation or permutation of external legs. Since our calculation is based on a numerical setup, these amount to independent calculations that give us an internal consistency check of our implementation.

V. CONCLUSION
The main result of this paper is a numerical implementation of two-loop five-gluon amplitudes in the leading-color approximation, for any helicity configuration. The calculation is done in the multi-loop numerical unitarity approach, which we extended for use with finite-field arithmetics. Using a new approach to the generation of IBP relations, we obtain a complete decomposition of the amplitude in terms of master integrals. Because of the extension to finite-field arithmetics, the corresponding coefficients computed from unitarity cuts are exact. To illustrate our calculation, we presented the results at a given kinematic point, for which we used the available master integrals [15,62]. Our implementation was validated by checking the universal pole structure [65], and confirms all results available in the literature [3,14,16,19].
The high polynomial degree and the many scales associated with multi-loop computations with several external particles make it challenging to set up accurate numerical approaches.
While analytic results are stable, the complexity of intermediate computational steps often makes them difficult to obtain for challenging amplitudes of that kind. Here, we chose an alternative approach which removes the strict separation of numerical and analytic perspectives and keeps the best features of each. We employed exact arithmetics as previously explored in the context of generalized unitarity in ref. [31]. Our approach is universal and is obtained after refinements of the method introduced in [24,25]. The exact approach is valuable for a number of reasons: the validation of results is simplified because the results are exact, which allows us to evaluate the amplitudes in singular limits, and, in addition, it paves the way for the functional reconstruction of their analytic form. This is left for future work.
The potential to generalize to other processes is significant, based on our experience with the present computation. The approach is general and requires little spectrum-dependent work. In particular, the parameterization of the integrand of the amplitude in terms of master integrals and surface terms, which allows to achieve a full reduction to master integrals, only depends on the kinematics of the process and the power-counting of the theory. The current setup handles five kinematic scales and the generalization to other massless particles is straightforward. Based on the present computation, we also believe that general five-scale processes are attainable and that adding further mass scales will be achievable in the near future. Finally, extensions to non-planar amplitudes appear well within reach.
Given the importance of analytic results for inspiring new methods and for identifying hidden symmetry principles, it will be exciting to learn about the analytic form of the integral coefficients. With the presented methods, we are well set up to explore this direction.

Renormalization of Leading-Color Two-Loop Gluonic Amplitudes
Renormalization of the amplitude is performed in the MS scheme. It is implemented by replacing the bare coupling by the renormalized one, denoted α S , in eq. (II.1). The bare and renormalized couplings are related through with S = (4π) e − γ E , where γ E = −Γ (1) is 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 = µ 2 0 = 1. For purely gluonic amplitudes, the coefficients of the QCD β function are The renormalized amplitude is written as with the renormalized A In the leading-color approximation, the color structure of loop corrections is the same as that of the leading-order contribution up to a factor of N C that was included in the perturbative expansion parameter, see eq. (A.3). For a n-gluon amplitude, the operator I (1) is then where s i,i+1 = (p i + p i+1 ) 2 with the indices defined cyclically. The operator I (2) is (A.7) The poles of the bare amplitudes can be recovered from those of the renormalized amplitude by using eq. (A.4). For amplitudes that are finite at one-loop (such as the all-plus and single-minus helicity configurations) the pole structure of the unrenormalized amplitude is particularly simple because the 1/ renormalization term in eq. (A.4) matches the 1/ term in the I (1) operator [14]. For example, in the all-plus five-gluon case we have