Double-Virtual NNLO QCD Corrections for Five-Parton Scattering: The Quark Channels

We complete the computation of two-loop helicity amplitudes required to obtain next-to-next-to-leading order QCD corrections for three-jet production at hadron colliders, including all contributions beyond the leading-color approximation. The analytic expressions are reconstructed from finite-field samples obtained with the numerical unitarity method. We find that the reconstruction is significantly facilitated by exploiting the overlaps between rational coefficient functions of quark and gluon processes, and we display their compact generating sets in the appendix of the paper. We implement our results in a public code, and demonstrate its suitability for phenomenological applications.


I. INTRODUCTION
Scattering processes that produce multiple jets in the final states are abundant at the Large Hadron Collider.In the next decades, precise theory predictions suitable for analyzing the expected large event samples will be crucial for advancing our understanding of high-energy interactions.Historically, pure-QCD predictions aimed at multi-jet production have also played an important role in the discovery of new structures and methods in field-theory, which inspire neighboring fields.
The goal of this work are predictions for the fiveparton processes in QCD, which have been an active field of research for many years.Tree-level five-point amplitudes in QCD were derived a long time ago in 1985 [1].About a decade later the one-loop amplitudes were obtained [2][3][4].The increase in analytic and algebraic complexity of two-loop computations required significant theoretical and technical developments over the following 20 years.This led to the five-point two-loop amplitudes for the all-plus helicity configuration to be first computed numerically [5] and later in analytic form [6,7].By now, the five-parton two-loop amplitudes with all helicity configurations are known in the leading-color approximation [8][9][10][11][12][13][14][15].
These results opened the door for the first computation of next-to-next-to-leading order (NNLO) QCD predictions for three-jet production [16] (see also [17]), and the follow-up measurements of the strong coupling constant at high scales [16,18,19].The double-virtual corrections [15], contributing on average about 10% [16], have been included in the leading color approximation in these works.This highlights the potential importance of including subleading-color effects.
In the preceding publication [20] we have obtained the analytic expressions for the five-point gluon am-plitudes including all color contributions.Computing the mixed quark and gluon channels is the central goal of this work, which completes the two-loop five-parton scattering amplitudes.
The computation of the non-planar five-point amplitudes is a multi-layered challenge, due to its analytic and combinatorial complexity.We rely on a number of recent results and methods.We use the massless five-point Feynman integrals [21][22][23][24][25], rely on the geometric methods for obtaining integrationby-parts identities [26][27][28][29][30], and we employ analytic reconstruction methods [31][32][33][34][35][36][37][38][39][40][41].In our computation we closely follow our recent work [20].We apply the numerical unitarity method [9,28,[42][43][44] to compute numerical values for the scattering amplitudes, which we use to reconstruct analytic results.The computation relies on the spinor-helicity formalism, which in turn yields very compact expressions as displayed in the appendices.In order to reduce the computational load, we develop an efficient way to obtain a large portion of quark amplitudes from gluon amplitudes.In fact, we rescale gluon amplitudes in a way reminiscent of supersymmetry Ward identities [1,45,46].The remaining parts of the amplitude are obtained applying a variant of a recent analytic reconstruction technique [44].
The analytic results for amplitudes are provided in ancillary files, suitable for future theoretical and phenomenological studies.In particular, we provide a C++ library for fast and stable numerical evaluation of our analytic results.Together with the gluon channel, presented in the first part of this work [20], these results will provide crucial input for NNLO predictions for three-jet and for N 3 LO di-jet production in hadron collisions.
Note added: while this work was in preparation, partially overlapping results were reported [47].We thank the authors of ref. [47] for correspondence on numerical comparison between our results.

II. NOTATION AND CONVENTIONS A. Helicity amplitudes
We consider the channels of the five-parton scattering in QCD that involve external quarks.We associate indices 1 and 2 to the initial state partons and 3, 4 and 5 to final states.There are three channels with a pair of quark and anti-quark, ū−h1 Without loss of generality we denote the quarks with u for up-type quark.There are four quark channels with distinct quark flavors, chosen to be up-type and down-type quarks, All other channels involving distinct quark flavors are related to the above by charge conjugation and permutation of labels.The three channels with four identical quark flavors are obtained as linear combination of the distinct quark ones, Representative diagrams for the two-loop contributions are collected in table I and II.
Throughout this article we use two interchangeable notations to specify the particles' helicities.We either give signs ± to label positive/negative helicities or, equivalently, the (half-)integers h = ±1/2 and h = ±1 to provide a better distinction between the quark and gluon helicities, respectively.
Whenever we use the arrow notation (→), the first two particles are to be understood as crossed to be incoming, such that their quantum numbers and momenta are given in incoming convention, while the final states are understood in out-going convention.
If no arrow (→) is used we consider the all-outgoing convention.
We will also use a shorthand for spinor chains, and we write tr 5 as 1 Under Lorentz transformations, spinor contractions have a residual covariance under little-group transformations.In fact, a little-group transformation of the i th leg with helicity h i reads (λ i , λi ) → (z i λ i , λi /z i ).Accordingly, helicity amplitudes transform as A → z −2hi i A. We refer to the exponent of the z i as the little-group weight.In summary, helicity amplitudes are homogeneous expressions of spinor brackets not just w.r.t. the total degree but also w.r.t each little-group weight.
Finally, throughout the article we will denote the group of cyclic permutations of n elements {i 1 , . . ., i n } by Z n (i 1 , . . ., i n ) and the set of all permutations of n elements, the symmetric permutation group, by S n (i 1 , . . ., i n ).

C. Color space
The external gluons are in the adjoint representation of SU (N c ), with indices denoted as a, b, c or a k , which run over over N 2 c − 1 values.Quarks carry (anti-)fundamental color indices i, j ( ī, j) or i k ( īk ), which assume N c values.Both color indices will be collected in tuples ⃗ a = {a 1 , . . ., a n , i 1 , ī1 , . . ., īn }.We explicitly represent the parton amplitudes in the color space through the trace basis as in ref. [49].
We use the shorthand notation tr(i 1 , . . ., i n ) = tr(T ai 1 • • • T ai n ), and (T a ) ī i are the hermitian and traceless generators of the fundamental representation of SU (N c ).
The generators T a are normalized as, and fulfill the commutator relations, The color algebra used in this paper is obtained from applying the Fierz identity, which allows to evaluate the summation over adjoint indices.
Using the above notation, the three-gluon twoquark scattering amplitudes are decomposed in color structures as Here σ = {σ 1 , . . ., σ 5 } denotes permutations which acts on all external-particle labels as σ(i) = σ i .The sums run over all permutations that do not leave the respective color structures invariant, i.e. over 6, 3, and 2 elements in the three lines of eq. ( 18), respectively.
Similarly, the one-gluon four-quark amplitudes are given by, where the sums run over exchanging quarks pairs.The amplitudes A i admit an expansion in terms of the bare QCD coupling constant α 0 s = (g 0 s ) 2 /(4π), with L denoting the number of loops.

D. Renormalization
We employ the 't Hooft-Veltman scheme of dimensional regularization to regularize ultraviolet (UV) and infrared (IR) divergences of loop amplitudes, where the number of space-time dimensions is set to D = 4−2ϵ.For helicity amplitudes with external quarks we employ the prescription of ref. [11].To cancel the UV divergences, we renormalize the bare QCD coupling constant in the MS scheme.This is accomplished by the replacement in eq.(20), where S ϵ = (4π) ϵ e −ϵγ E , with γ E = −Γ ′ (1) the Euler-Mascheroni constant, and µ 0 , µ are regularization and renormalization scale respectively.The QCD β-function's expansion coefficients are We then expand the renormalized amplitudes through the renormalized coupling as in eq.(20).
The remaining divergences are of IR origin and can be predicted by the universal factorization [52][53][54][55]: where the finite remainder R is obtained by the application of the color-space operator Z, which is derived by exponentiation [54], of the anomalous dimension matrix Here the sum runs over pairs of external partons, and the color-space operators T i act on the color representation of the i th parton.For adjoint indices the action is given by (T a i ) bc = −if abc , and for fundamental indices as (T a i ) k j = ±(T a ) k j ; n q and n g are the number of quarks and gluons in the process respectively, and the functions γ cusp , γ q = γ q , γ g can be found in ref. [56, Appendix A]2 .
After absorbing all divergences of amplitudes into UV and IR renormalization through eqs.(23) and (25) we recover expansions of eqs.(18) to (22) at the level of finite remainders R. We then define partial finite remainders which will be the elementary building blocks considered in this work.

E. Generating set of finite remainders
Using parity, charge conjugation transformations and permuting momentum assignments to the external states, we find a generating set of finite remainders that we need to compute.Identities inherited from color factors allow one to further restrict the set of required functions, as discussed below section IV.
Focusing first on two-loop data and suppressing the labels specifying the N c and N f decomposition we consider a generating set of helicity assignments.To this end we generate all helicity assignments, and select one representative from the orbits of the charge conjugation and parity transformations.Furthermore, given that we are considering analytic amplitudes, we chose a convenient momentum assignment for each helicity amplitude, which we line up with the little group transformation properties in the single-minus and the MHV amplitudes.
For the single-minus helicity configuration we have For the MHV configurations we have five generating remainders, The analogous analysis yields the following generating set of finite four-quark remainders with the MHV helicity configuration, In section IV we discuss further identities between remainders associated to distinct terms in the N c ,N f decomposition.

F. NNLO hard function
The two-loop NNLO QCD corrections for partonic cross sections are obtained from squared helicityand color-summed partial remainders, which we call hard functions H, Here the summation is performed by mapping each partial remainder into one from the generating sets eqs.( 29) to (31).The hard function is expanded perturbatively up to O α 2 s as in eq. ( 20).We further expand H in powers of N f , while we keep the value of N c implicit,3  [2] .(33c)

G. IR scheme change
It is important to highlight that the finite remainders encompass all the physical details related to the underlying scattering process.Specifically, one can compute any observable by utilizing finite remainders (see e.g.[57]).In this way it is possible to cancel many undesirable side effects of dimensional regularization.
In this work we use the minimal subtractions scheme of IR renormalization, following the conventions of refs.[54,56].One might be interested in obtaining finite remainders defined in a different IR renormalization scheme, e.g.Catani's scheme [52,[54][55][56].In the following we show that it is possible to convert the finite remainders in our scheme to any other scheme by an additional finite renormalization, i.e. the knowledge of higher orders in ϵ of amplitudes is only initially required to derive finite remainders in arbitrary scheme.We take advantage of this fact in our computational framework and circumvent analytic reconstruction of amplitudes.
Suppose we are interested in a different scheme where the finite remainder is defined as where A is the same UV renormalized amplitude as in eq. ( 25).Here we remind the reader that A and R are vectors and Z is an operator in color space.We assume that both Z and Z have a perturbative expansion as in eq. ( 20), and We then consider the difference which is finite by definition.Therefore the operator δZ = Z − Z Z −1 must not have ϵ poles, except possibly the ones that cancel upon action on R. Provided the latter cancellation does not rely on the existence of a non-trivial null space, both δZ and R and can be truncated at O ϵ 0 .We can therefore express the remainder R through R order-by-order in α s by a finite renormalization δZ, whose perturbative expansion starts at O(α s ).
For the squared finite remainders we can write more explicitly through two loops ⃗ h,⃗ a .
It is then straightforward to perform helicity and color summation to derive hard functions in the new scheme.
We have explicitly calculated the operator δZ to convert the minimal subtractions scheme to the Catani scheme.We verified that the latter must be supplemented by both types of tripole color correlation terms added in the later revisions of [56, eq. ( 17)] for the poles to cancel at the level of partial remainders for five-parton scattering.In addition, we have verified numerically that the O(α 2 s ) contributions to the squared helicity-and color-summed matrix elements originating from the tripole correlation terms vanish in all five-parton channels.

III. NUMERICAL SAMPLING OF REMAINDERS
We will construct the partial remainders from analytic reconstruction, i.e. we compute the analytic form of the partial remainders from numerical evaluations in a finite field.Partial remainders can be expressed as a linear combination of transcendental integral functions h i and rational coefficient functions r i , The integral functions, referred to as pentagon function, are known [25].The computation of the analytic form of the function coefficients r i is one of the central results of this paper.
Here we summarize the input data required for our computation following ref.[20].For numerical evaluation of remainder functions in a finite field we use the program Caravel [43], which implements the multi-loop numerical unitarity method [9,28,42].In this approach amplitudes are reduced to a set of master integrals by matching numerical evaluations of generalized unitarity cuts to a parametrization of the loop integrands.For the five-parton process we use the recently obtained non-planar parametrization [20,44].Furthermore, for the quark processes we extended the set of planar unitarity cuts to nonplanar diagrams which are required for subleadingcolor partial amplitudes.We generated the cut diagrams with qgraf [58] and arranged them into a hierarchy of cuts with a private code.The unitarity cuts evaluated through color-ordered tree amplitudes are matched to the amplitude definitions in eqs.( 18) and (19), by employing the unitarity based color decomposition [59,60].The ϵ-dependence of cuts that originates from the state sums in loops are obtained by the dimensional reduction method developed in ref. [14,61,62].With these upgrades Caravel now computes the function coefficients r i of five-parton partial amplitudes up to two loops, given a kinematic point and a choice of polarization labels for the external gluons.
Next we will require two types of numerical samples of the remainder functions which we repeat from ref. [20] for convenience: 1. Random phase-space points: these are N randomly generated phase-space points which we label by the superscript n.We represent these points in terms of sets of spinor variables, They are subject to momentum conservation i λ n i λn i = 0. Below we will use the shorthand notation, to denote the spinor-helicity variables associated to a phase space point.
Here the reference spinor η is chosen randomly.The constants c i are obtained by solving the linear momentum-conservation condition.In addition we will use the anti-holomorphic slice which is obtained from eq. ( 40) by swapping λ i ↔ λi , η → η and renaming t → t.

IV. IDENTITIES BETWEEN PARTIAL AMPLITUDES
Partial amplitudes in the trace-basis representation are known to satisfy linear relations originating from symmetry properties and the adjoint representation gauge interactions of field theory.These relations, which we will refer to as color identities, can be exploited to reduce the number of partial amplitudes that need to be computed, and subsequently to improve the efficiency of the numerical evaluation of the hard functions (32).
Color identities were discussed for multi-loop gluon amplitudes [64] and recently in ref. [65].Much less is known about the scattering of quarks and gluons at two loops.Here we empirically identify all linear relations for the five-point amplitudes including quarks from numerical evaluations.We proceed as follows.
Linear relations can only hold between remainders with the same little-group weight, which we specify below by the labels ⃗ h = {h 1 , ..., h 5 }.We group all remainders with identical little-group weights into sets where the index i enumerates the S remainder functions.The superscript ⃗ λ specifies the momentum of the states (see eq. ( 38)).We employ random numerical samples (38), which associate a vector of numerical values to each remainder, We then search for vanishing linear combinations of these vectors, The set of nontrivial constant solutions c i yield the desired identities between partial remainders.
Technically, we simplify the search for identities by using finite-field arithmetic and exploiting the fact that the transcendental functions in the decomposition (37) form a basis [66].We first obtain identities between coefficients of selected transcendental functions.We then intersect them in order to find an identity valid for the entire remainder.
In principle we can distinguish two classes of identities: 1) helicity dependent ones, which hold in a single class of helicity configurations, i.e. in the single-minus or MHV configuration.2) Helicity independent relations, which hold for all helicity assignments.Color identities are of the second type.At two loop we find only helicity independent identities.
Let us note that algorithms to obtain color decompositions and relations are well understood for arbitrary multiplicity at one loop [50,[67][68][69][70].In particular, representations of two-quark and four-quark amplitudes in terms of so called primitive amplitudes are known [3,4,50,71], which imply the identities that we study here.

A. Two-quark channel
We consider first the remainders associated to the partial amplitudes of A(1 h1 u , 2 h2 ū , 3 h3 g , 4 h4 g , 5 h5 g ) eq. ( 18).To start with, we collect the set of all remainders with identical little-group weights.These are obtained from evaluating the remainders on permuted momenta, keeping helicity quantum numbers assigned to each momentum fixed.With this in mind, the full group of permutations is generated by the groups S 2 (1 h1 , 2 h2 ) and S 3 (3 h3 , 4 h4 , 5 h5 ), which do not mix gluons and quarks.In total, the permutation group contains 2 × 6 = 12 elements.This set of remainders can be further reduced, using symmetry properties of the color factors: R 1 : Charge conjugation symmetry considered for the amplitude A(1 h1 u , 2 h2 ū , 3 h3 g , 4 h4 g , 5 h5 g ) ( 18) forces the coefficients of the color factors (T a3,a4,a5 ) ī2 i1 and (T a5,a4,a3 ) ī1 i2 to match, e.g.R 1 (1, 2, 3, 4, 5) = −R 1 (2, 1, 5, 4, 3).This relation halves the number of independent momentum permutations.We chose representatives generated by S 3 (3 h3 , 4 h4 , 5 h5 ), namely R 2 : Using charge conjugation and the cyclicity of the trace tr(3, 4) = tr(4, 3) implies that partial amplitudes A 2 (and their remainders) are unchanged under S 2 (3 h3 , 4 h4 ), and S 2 (1 h1 , 2 h2 ).Its symmetry group thus has dimension 4, meaning that out of the total 12 permutations of momenta (and helicities) there are 3 inde-pendent permutations.We chose the representatives R 3 : Finally, charge conjugation symmetry and cyclicity of the trace tr (3,4,5) implies invariance of the partial amplitudes A 3 (and R 3 ) under the transformations Z 3 (3 h3 , 4 h4 , 5 h5 ) and S 2 (1 h1 , 2 h2 ).After modding the 12 total permutations by these 3 × 2 = 6 symmetry transformations 2 independent permutations remain, which we pick to be Here we suppressed again the superscripts specifying the partial remainder, namely their loop order and contribution in N c and N f in eq. ( 18).We have obtained a set of two-loop partial remainders with identical little-group weight, All remainders in this set of remainders can be expressed in terms of the generating set of section II E.
After following the steps discussed in the beginning of this section we find no nontrivial relations in the two-loop two-quark channel.
As a final point, we note that the identities that we have found do not allow us to express any of the partial remainders in terms of sums over permutations of the others, i.e. the generating set discussed in section II E cannot be further reduced for any n c , n f .This is in contrast to the well-known fact that in the five-gluon channel the most subleading in N c expansion partial amplitude can be eliminated [64].

V. ANALYTIC RECONSTRUCTION
As discussed in section III we have available numerical evaluations of the coefficients r i in the remainder function (58).We will now build upon ref. [20] to obtain compact analytic expressions for the function coefficients from such numerical samples.The starting point is the understanding that the coefficients admit the least common denominator representation, where the denominator factors D j are given by the letters in the symbol alphabet of pentagon functions [23,24,72] with integer exponents q ij [13].The goal of analytic reconstruction is to determine q ij and the numerator polynomials N i in eq. ( 55).First we determine the exponents q ij .To this end we follow the univariate-slice reconstruction [13] in spinor-helicity variables [37,44,63].In this approach the function coefficients are obtained as univariate rational functions r i (t) and r i ( t ) on a holomorphic and an anti-holomorphic slice (40), respectively.Given the rational functions r i (t) and r i ( t ), their denominators are matched to products of the letter polynomials D j (t) and D j ( t ).This uniquely fixes the exponents q ij in each of the functions r i (55).In particular considering holomorphic and anti-holomorphic slices independently, ensures that one identifies purely (anti-)holomorphic terms, such as [ij] and ⟨i j⟩.
The importance of the exponents q ij is two-fold.On the one hand, we determine which of the letters actually appears as denominator factors.For the five-parton finite remainders we observe that the set of denominator factors consists of the 35 elements, where the set runs over all independent permutations of the spinor strings/brackets.None of the coefficients in the remainder has a tr 5 singularity.
On the other hand, the exponents constrain the ansatz (55), since the mass dimension and littlegroup weight of the numerator polynomial N i is uniquely fixed by those of the helicity amplitude and the denominator.Consequently, we can construct a finite dimensional ansatz for the polynomial N i .We have thus reduced the computation to the problem of finding finitely many polynomial parameters.
Before we obtain the functions r i we wish to identify linear dependences, to identify a minimal set of functions that we need to compute.To this end we first sort the functions r i according to complexity, namely, the mass dimension of the respective numerators N i , which in turn is correlated with the polynomial's parameters.Next we identify linear dependence numerically via Gaussian elimination, and determine the indices of the basis coefficient functions in the set B, In this way we further reduce the data, required to specify the scattering process to a basis of rational coefficient functions r i and a constant rationalvalued matrix M ij [14], The next task is to determine the set of numerator polynomials {N i } i∈B for the basis functions of all partial remainders from the random numerical evaluations of eq. ( 38).Naively, determining N i requires as many evaluations as there are free parameters in the polynomials N i .A large number of evaluations can be limiting due to the evaluation time of partial remainders.Reducing the size of the required numerical sample is thus an important goal.Below we exploit two observations about the structure of the coefficient functions r i to significantly reduce the size of the required numerical samples: In section V A we recycle the compact coefficients of the two-loop five gluon amplitudes [20] and, in section V B we exploit linear basis changes to find new coefficients in partial fractioned form.

A. Rescaled coefficient functions
We construct a class of candidate coefficient functions from the known functions for the five-gluon amplitudes [20].To this end we rescale the gluon coefficient functions by spinor brackets to match the little-group weight of the quark amplitudes.This rescaling is inspired by analogous factors in supersymmetry Ward identities which link tree amplitudes of gluon and gluino states [1,45,46] (see also [73]).Such a rescaling is expected from on-shell recursion relations [74] for loop amplitudes [75][76][77][78] and collinear factorization [79][80][81] in general.Let us start from a gluon function from appendix C of ref. [20], e.g.
First, in order to align the little-group weights with the two-quark-three gluon amplitudes, we permute momentum labels Then, we can build functions for the process (u 1/2 , ū−1/2 , g 1 , g 1 , g −1 ) by multiplying the gluon function by any function which raises and lowers the little-group weights of legs 1 and 2 by one unit, respectively.For instance, functions of the form correctly map the little-group weights.The function that we obtain is We then test numerically whether the resulting function belongs to the vector space spanned by the coefficients of the (u 1/2 , ū−1/2 , g 1 , g 1 , g −1 ) MHV partial remainders, If it is in fact part of the span, we keep the function.
In this way we obtain a set of analytically known functions which allows to parametrize part of the coefficient functions r i .We denote the set of these functions by and index them by the label i in the set B. To simplify the set, we remove linearly dependent functions ri .
Let us note that we do not aim here to explore all possible rescaling factors.For simplicity, we only consider factors such as that of eq. ( 61) and generalizations thereof with numerator and denominator mass dimensions not exceeding two.
After applying this rescaling procedure we obtain a significant number of coefficient functions of the quark amplitudes from the gluon ones.We obtain more than 50% of the two-quark three-gluon MHV functions by rescaling the five-gluon functions.Since the basis of two-quark three-gluon single minus functions can be reconstructed from a small number of sample evaluations, we do not apply this strategy for them.We also obtain the majority (more than 90%) of the four-quark one-gluon basis functions by rescaling two-quark three-gluon and the five-gluon functions.

B. Filling the space of coefficient functions
So far we have obtained a portion of the coefficient functions from rescaling gluon coefficients, as discussed above in section V A. We will now reconstruct new rational functions which we add to the set {r i } i∈ B of eq. ( 64) until it spans the full function space {r i } i∈B of eq.(57).In order to simplify the discussion, we now assume that the sets B and B correspond to a specific helicity class, i.e. single minus or the MHV configuration.By construction, the functions ri are in the linear span of the r i .In terms of the vectors of function values, we have, In the reconstruction of the missing functions, we now leverage the observation that partial fractioned coefficients take a simple form, with reduced mass dimension of numerator and denominator.We select the r i with lowest mass dimension numerator, that does not lie in the span of {r i } i∈ B .We then solve the partial fractioned ansatz [44] The parameters qij in the ansatz are set empirically, i.e. we set a degree bound on the sum of powers j qij ≤ q.We then walk through all possible choices of denominators j D qij j in eq. ( 66) of the chosen degree obtained as subsets of the original denominator j D qij j in eq. ( 55) and fit the numerator polynomial Ñi and coefficients c ik using the numerical sample evaluations from eq. (38).Given a successful fit, we then consider the function, ri = Ñi and all associated functions with permuted momentum labels and matching little-group weights.We add all such functions to the set of new, analytically known functions.We denote the updated set again with {r j } j∈ B with an adjusted index set B. The procedure is repeated until the linear span of the reconstructed functions covers the one of the coefficient functions, By construction, the only overshoot of the span in the left-hand side compared to that in the right-hand side can be in the permutation closure of the generating functions.No generating spinor-helicity function can be dropped while keeping eq. ( 68) valid.Thus, the cost of this overshoot is minimal, while potentially being helpful to further simplify the basis.We observe that this approach is very effective.The degree bound of q = 12 suffices for the analytic reconstruction of all remaining coefficient functions from approximately 2k random numerical samples (38).
The method is efficient due to a number of implementational improvements: we cache numerical values for both the right-hand side and the summation of eq.(66).For each choice of qij , we only need to regenerate an ansatz for N i , insert numerical values, and perform a Gaussian elimination.For ansätze of this size, both construction via OR-tools and row reduction via linac take O(100 ms), meaning hundreds to thousands of guesses can be checked in the time it takes to collect additional numerical samples for the remainder functions.

VI. RESULTS
We have obtained the analytic expression for the two-loop two-quark three-gluon and four-quark onegluon helicity finite remainders.The results are given in ancillary files as explained in detail in section VI C.
One of the central results of this work is the basis of rational coefficient functions ri , which we given in the appendix of this paper.We label the fourquark one-gluon MHV functions simply as r and present them in appendix B. The two-quark three gluon functions are labelled based on the helicity of

Particle Helicities
Vector-space dimension Generating set size Table III: For each helicity configuration, this table shows the dimension of the vector space of rational functions, and the number of functions in the generating set that spans the space upon closure under the symmetries of the little-group weights.the last gluon, as r+ for the single minus configuration, and r− for the MHV one.We present them in presented in appendix C and appendix D, respectively.An account of the size of the function bases is displayed in table III.
In order to facilitate future comparison with our results we collect reference values for finite remainders in appendix A. Finally, in section VI B we present an efficient implementation of our results and discuss its performance and stability.

A. Validation
In order to validate the amplitude computation we have performed a number of checks.The evaluation of the amplitudes with the numerical unitarity method was carried out with the well-tested program Caravel [43].At each phase-space point used in analytic reconstruction we check cancellation of ϵ poles in finite remainders (28).To validate the analytic reconstruction, we check that the analytic results match further evaluations in Caravel using finite-fields with a different characteristic.
Furthermore, we perform a number of checks on the hard function (32), which verifies the assembly of partial remainders based on color and helicity sums.We verified that all expected symmetries of the hard functions (32) under permutations of particles' momenta are satisfied, as well the parity-conjugation symmetry.We have compared our numerical evaluations of the one-loop hard functions H (1) (33) against BlackHat [78] in all physical channels and found perfect agreement.We have verified that upon taking the leading-color limit of H (1) , H (2) , we reproduce (after performing the IR scheme change) the results resented in [15].Finally, we have compared our results after performing the IR scheme change (see section II G) with the numerical benchmarks presented in [47,Table II].We find agreement with a revised version of [47].

B. Numerical evaluation
We implement our analytic expressions for the generating set of helicity partial remainders in eqs.(18) and (19), as well for the hard functions in eqs.(32) and (33) for all physical channels in eqs.( 1), ( 4) and ( 8) in a C++ library [82].Together with the five-gluon channel from ref. [20] this is the complete set required for the calculation of NNLO QCD corrections for three jet production at hadron colliders.
To demonstrate the numerical performance of our implementation, we sample three representative channels over 100k points from the phase space of ref. [83] (see also [15]).For completeness, we also study the five-gluon channel calculated in [20].The evaluations are compared to target values computed in quadruple precision, and the distribution of the base 10 logarithm of the relative error (correct digits) is shown in fig. 1 (c.f.ref. [15]).We observe that despite the markedly increased complexity of subleading color amplitudes, the numerical performance of our results is excellent and comparable to the leading-color results reported in [15].The rescue system developed in [15] is effective in capturing unstable points.This is especially relevant for the five-gluon channel where we observe the second small peak in the distribution at about 16 digits formed by the phase-space points rescued through a quadruple precision evaluation.
We observe (in the panels of fig. 1) average evaluation times per phase-space point of few seconds in all production channels, already taking into account the time required to detect and rescue unstable points in quadruple precision.In contrast to the leading-color approximation, where most of the evaluation time is spent on pentagon functions, the evaluation in full color is dominated by the contraction of indices in eq. ( 58).This hints that further improvements are achievable if a more tailored basis of transcendental functions is used.The observed stability and evaluation times will enable seamless lifting of the leading-color approximation which has been employed in cross-section computations so far [16,17].

C. Ancillary files
We provide ancillary files for all independent partial finite remainders, including all crossing [84].We organize the ancillay files in the same manner as ref.[20].Overall, for the complete five-parton computation, the folder structure is as follows: For each of these folders, representing external states and the associated helicity configuration, we provide the bases h j and ri (58), respectively in the files 1. basis transcendental , 2. basis rational .
Further subfolders contain the matrices M ij of rational numbers (58).The folders are labelled as where ⃗ h, L, n c and n f refer to helicities, number of loops, number of N c powers and number of N f powers.Since L, n c and n f do not always identify a unique partial, for A 2 and A 3 we extend this notation to where the extra "integers" represent the split in the fundamental generators, i.e. 2 1 for A 2 and 3 0 for A 3 .The rational matrices are named as rational matrix {permutation} , for each permutation of the external legs, which may involve crossings.
We further provide assembly scripts within the ancillary files.

VII. CONCLUSIONS
We have presented the computation of the twoquark three-gluon and four-quark one-gluon amplitudes at two loops in QCD.We derive compact analytic expressions in the spinor-helicity formalism for the finite remainders, including all contributions beyond the leading-color approximation and all crossings.
We systematically investigate linear relations among partial remainders in the trace basis of the color generators, relying on numerical amplitude evaluation and linear algebra.We do not find any non-trivial identities among the two-loop two-quark three-gluon partial remainders, while we obtain six identities among two-loop four-quark one-gluon partial remainders.
With regards to the analytic reconstruction, we explore a new method to obtain quark amplitudes from gluon ones, inspired by supersymmetry Ward identities.This entails rescaling the gluon spinorhelicity basis functions presented in [20] by simple factors carrying the appropriate little-group weights, such as ratios of the respective tree amplitudes.
Finally, we provide the efficient C++ code [82] for the computation of color-and helicity-summed squared matrix elements, suitable for immediate phenomenological applications.
Together with our earlier results of the compact gluon amplitudes [20] this completes our computation of the two-loop five-parton amplitudes in full color.We envisage their application to both new phenomenological studies, as well as novel theoretical investigations into the perturbative structure of QCD.32) for representative physical channels contributing to three-jet production at hadron colliders.The dashed curves represent respective cumulative distributions.Here N f is set to 5, and the renormalization scale is set dynamically to (half of) the sum of the transverse momentum of the final-state partons.The phase space definition is taken from ref. [83].

Figure 1 :
Figure1: Distributions of the base 10 logarithm of the relative error (correct digits) of the NNLO hard functions defined in eq.(32) for representative physical channels contributing to three-jet production at hadron colliders.The dashed curves represent respective cumulative distributions.Here N f is set to 5, and the renormalization scale is set dynamically to (half of) the sum of the transverse momentum of the final-state partons.The phase space definition is taken from ref.[83].

Table I
Table II: Representative Feynman diagrams for four-quark one-gluon amplitudes, contributing at order N ] 2 ⟨45⟩