Leading-color two-loop QCD corrections for three-jet production at hadron colliders

We present the complete set of leading-color two-loop contributions required to obtain next-to-next-to-leading-order (NNLO) QCD corrections to three-jet production at hadron colliders. We obtain analytic expressions for a generating set of finite remainders, valid in the physical region for three-jet production. The analytic continuation of the known Euclidean-region results is determined from a small set of numerical evaluations of the amplitudes. We obtain analytic expressions that are suitable for phenomenological applications and we present a C++ library for their efficient and stable numerical evaluation.


Introduction
Multi-jet events are ubiquitous at high-energy hadron colliders such as the Large Hadron Collider (LHC).Their precise theoretical description including higher-order quantum corrections offers a robust path to test the Standard Model of particle physics.This has been demonstrated, for example, by comparing measurements of the strong coupling α s extracted from ratios between three-and two-jet production cross sections with theoretical predictions computed at next-to-leading order (NLO) in QCD [1][2][3][4].Given the amount of experimental data that has been collected since these studies were conducted and the amount of data forecast to be collected in the near future at the LHC, there is great interest in obtaining theoretical predictions that will allow these types of studies to be carried out at a new level of precision.In this paper, we provide a key missing ingredient required to reach this goal, namely the leading-color two-loop contributions to the three-jet production cross section at hadron colliders at next-to-next-to-leading order (NNLO) in QCD.These are the so-called double-virtual 2 → 3 contributions, which together with the real-virtual 2 → 4 and the double-real 2 → 5 contributions complete the O(α 2 s ) corrections to the production of three jets at hadron colliders.
The calculation of NNLO QCD corrections to two-jet production at hadron colliders has been a major success of the field in recent years.Results were obtained first at leading color [5][6][7] and later subleading-color terms were incorporated [8] (see also [9]).Three-jet production is a more complex process.Results are currently known at NLO in QCD [10,11] (see also [12,13]) and also including electroweak corrections [14,15].Work on the NNLO QCD corrections to three-jet production has so far been reported only for the leading-color two-loop scattering amplitudes, first through numerical evaluations [16][17][18][19][20] and then in analytic form [21][22][23][24][25]. Progress in including subleading-color contributions has been made for related supersymmetric theories [26][27][28] as well as for special external gluon states [29].In this work we present the leading-color double-virtual contributions to the scattering cross section, starting from the known analytic amplitudes [24] which we analytically continue from the (non-physical) Euclidean region to the physical region of phase space.Furthermore, we implement our expressions in a publicly available C++ library that allows one to evaluate these contributions efficiently and in a stable manner, making our results readily available for a future full calculation of the corresponding NNLO QCD corrections. 1ur work builds on a number of recent technical advances that have paved the way to the calculation of two-loop scattering amplitudes for 2 → 3 processes with phenomenological applications [32][33][34][35][36][37].We employ the known pure two-loop master integrals [21,26,38,39], which can be written in terms of a special set of multi-valued transcendental functions-the so-called pentagon functions [40,41].In ref. [41], the five-point integrals are expressed in terms of a basis of functions for all assignments of the incoming momenta required for the double-virtual corrections we compute, and they can be efficiently evaluated in the physical region corresponding to three-jet production.
We obtain the analytic form of a generating set of partial amplitudes in the physical region by analytically continuing the known amplitudes from the Euclidean region [24].Based on considerations on the analytic continuation of the pentagon functions, the continuation of each helicity amplitude is reduced to a linear-algebra problem that can be efficiently solved using evaluations of the amplitudes in a finite field [18,42,43].These evaluations are performed with the program Caravel [44].This is a new generic procedure that allows to upgrade results for amplitudes in the Euclidean region to results valid in other regions.Finally, we assemble all partial amplitudes into squared color-and helicity-summed finite remainders, which can be numerically evaluated over phase space with the provided C++ library [45].We perform several checks on our results and verify the efficiency and numerical stability of our implementation.
The paper is organized as follows.In section 2 we establish our conventions and define the objects we will be computing.In section 3 we describe our procedure to compute the partial amplitudes in the physical region and the different permutations that are required to obtain squared finite remainders.In section 4 we present our results in the form of ancillary files and showcase the efficiency and numerical stability of our code for the nu-merical evaluation of the double-virtual NNLO contributions to the production of three jets at hadron colliders, and discuss the checks that we have performed on our results.We summarize our results in section 5.

Helicity Amplitudes
We consider all the channels contributing to the production of three partons at hadron colliders.We label the initial-state particles with indices 1 and 2, and the final-state particles with indices 3, 4 and 5.There is a single channel with five gluons, three channels with a pair of quark and anti-quark, and four channels with two distinct pairs of quark and anti-quark In eqs.(2.1) to (2.3), we identify each particle by its type (g for gluon, q or Q for distinct quarks, and q and Q for the corresponding anti-quarks), a superscript h i for the helicity state and a subscript p i for the four momentum.All other channels involving quarks are related to the ones above by charge conjugation and permutation of incoming or outgoing labels.While in eq. ( 2.3), we assumed the two quark pairs to be of distinct flavours, we should also consider the case with two quark lines of identical flavour.There are three such channels which, schematically, can be obtained from In order to obtain them we thus also compute the (redundant) channels We remark that, when labelling channels, such as in eqs.(2.1) to (2.5), the momenta p 1 and p 2 are taken to be incoming and p 3 , p 4 and p 5 are taken to be outgoing.In all other contexts we use the symmetric "all-outgoing" convention of refs.[20,24], such that momentum conservation is p 1 + p 2 + p 3 + p 4 + p 5 = 0. We will consider each of the channels above in the leading-color approximation.Let us denote the set of helicities by a vector h and the set of momenta by a vector p.The amplitudes M(h, p) describing each of the channels above, which are vectors in color space, can be expanded in powers of the bare strong coupling α 0 s = (g 0 s )2 /(4π), where we drop the dependence on h and p for compactness.The order of α 0 s is aligned with the number of loops of the Feynman diagrams contributing to the associated coefficient, so we refer to M (1) and M (2) as the one-loop and two-loop amplitudes respectively, and to M (0) as the tree-level amplitude.All calculations will be done in the 't Hooft-Veltman (HV) scheme of dimensional regularization with D = 4 − 2ϵ space-time dimensions and we mostly keep the dependence on ϵ implicit.For our definition of helicity amplitudes with external quarks, we refer the reader to the detailed discussion in ref. [20].
For each channel, the coefficients M (i) can be decomposed in terms of different color structures, and at leading color this decomposition is independent of the loop order.We write where N c is the number of colors, C (σ) denotes a vector in color space and the sum is over a channel-dependent set of permutations Σ.We defer a more detailed discussion of the sets Σ and the vectors C (σ) to section 3.2.For now, we simply define the action of a permutation σ = {σ 1 , σ 2 , σ 3 , σ 4 , σ 5 } ∈ S 5 on momenta and helicity labels as σ(p i ) := p σ i and σ(h i ) := h σ i .The Φ (σ(h), σ(p)) are normalization factors, and the A (i) (σ(h), σ(p)) are partial amplitudes.
For each term in the sum, the coefficient of each C(σ) has a fixed ordering of the external momenta determined by the permutation σ.The Φ (σ(h), σ(p)) are defined so that A (0) (σ(h), σ(p)) ≡ 1 whenever the tree-level amplitude is non-vanishing, and for the other cases we choose 2 qq→Q Qg ∼ + + + . . . .where ε(•, •, •, •) is the fully anti-symmetric Levi-Civita tensor, whose sign is fixed such that Im(tr 5 ) > 0 corresponds to the choice of the basis helicity amplitudes in ref. [24].The action of σ on quantities that depend on momenta is inherited from the action on momentum labels.For instance, the action on the Mandelstam invariants s ij = (p i +p j ) 2 is given by σ(s ij ) = s σ i σ j , and on the parity-odd invariant in eq.(2.10) by σ(tr 5 ) = sgn(σ) tr 5 .
We work in the leading-color approximation where the number of colors N c is large and the ratio N f /N c is kept fixed, where N f denotes the number of massless quark flavors.In this limit, each of the A (i) (h, p) can be further decomposed in powers of N f .We write (2.11) In fig. 1 we depict characteristic diagrams that contribute to A (2) [j] for representative channels with no external quark lines, a single external quark line or two external quark lines.
The renormalized amplitudes M R are obtained from the bare amplitudes M by considering an expansion in powers of the renormalized coupling α s (µ).In the MS-scheme, the later is related to the bare coupling α 0 s through where S ϵ = (4π) ϵ e −ϵγ E , and µ 0 and µ are the dimensional regularization and renormalization scales respectively, which from now on we assume to be equal.We will also suppress the µ dependence in the coupling and write α s ≡ α s (µ).The β i are the coefficients of the QCD β-function, where we have set ), and only kept terms that contribute at leading color.The renormalized partial amplitudes are related to their bare counterparts as (2.14)

Finite Remainders
Amplitudes contribute to physical observables only through so-called finite remainders (see e.g.[46]), which we denote R (h, p).We define a remainder for each partial amplitude A (h, p).They are obtained from the renormalized amplitudes by removing the remaining singularities of infrared origin that can be determined from lower-loop amplitudes and well-known universal factors [47][48][49].The finite remainders can be expanded in powers of α s , with the R (i) defined as The operators I (i) are channel specific, and we collect them in appendix A for all the channels in eqs.(2.1) to (2.3).In these expressions, we extend the expansion of R to also include terms of order ϵ 0 .This subtracts non-trivial contributions from the finite term of R that are related to the lower-loop amplitudes.It is clear that the finite remainders R (k) contain the genuine new information at order k.We stress that even though our calculations are performed in the HV scheme of dimensional regularization, the remainders computed in this scheme agree with the ones computed in the conventional dimensional regularization (CDR) scheme, see e.g.[50].Finally, the finite remainders can also be decomposed into powers of N f , similarly to eq. (2.11), (2.17) Since the functions R (i) [j] are dimensionless by construction, their scale dependence can be restored at any point by rescaling the momenta as p → p/µ (we set the factorization scale µ F to be equal to the renormalization scale, µ F = µ).Unless it happens to be important for the point we wish to make, we will keep this dependence implicit.The color-dressed remainder R M (h, p) is trivially obtained from the decomposition of the color-dressed amplitude M (i) (h, p) into partial amplitudes, see eq. (2.7): (2.18) Whenever the distinction is important, we will call the R (i) (σ(h), σ(p)) partial remainders, in analogy with the concept of partial amplitudes.We also recall that, even though this is kept implicit, the color-dressed remainders are vectors in color space (in general, the C (σ) have open (anti-)fundamental and adjoint color indices, see section 3.2 for their explicit form).

Squared Finite Remainders
When computing physical observables, we are interested in the squared remainder, summed over color states and, in the case of unpolarized observables, over helicity states.We denote this object as H.It admits an expansion in powers of α s as In terms of the quantities in eq.(2.18), H is defined as where the normalization N H is fixed such that H (0) = 1, (2.21) To keep the notation light, in this equation and the following we do not write explicitly the sum over contracted colour indices which is understood to be performed in the squaring procedure.At order α s we have and at order α 2 s M (h, p) , (2.23) Figure 2: Schematic representation for the two kinds of contributions to H (2) : the square of one-loop amplitudes, and the interference between tree-level and two-loop amplitudes.
with each of the two contributions depicted in fig. 2. We note that, when summing over contracted color indices, we have C † (σ ′ ) C(σ) ∼ δ σσ ′ in the leading-color approximation we adopt.This greatly simplifies the construction of the squared remainders.
To be clear about which terms we keep in the leading-color approximation, we explicitly write the decomposition of our remainders in powers of N f /N c .For H (1) , we have For H (2) , For most channels the subleading-color corrections are suppressed by a factor of 1/N 2 c .However, for the channels in eq.(2.4) with two identical quark lines the suppression is by a single power of N c .

Physical Region Helicity Amplitudes
Our goal in this paper is to obtain the double-virtual contributions for three-jet production at hadron colliders at NNLO.The central object we need to compute is the squared finite remainder defined in section 2.3.In this section, we discuss how we construct the different contributions starting from the Euclidean-region partial amplitudes obtained in ref. [24].We recall the comment below eq.(2.5) regarding the change of conventions for denoting each channel between this paper and ref. [24].

Partial Amplitudes in the Physical Region
The first step towards the construction of the squared finite remainder in the physical region is to obtain the amplitudes computed in ref. [24] in such region.Here we consider only the basis set of partial helicity amplitudes of ref. [24], and we will return to their permutations required to assemble eq.(2.23) in section 3.2.The amplitudes of ref. [24] were computed in the (non-physical) Euclidean region where The latter condition is trivially satisfied by real-valued momenta.In this section, for clarity, we will denote a phase-space point in the physical region as p φ and a phase-space point in the Euclidean region as p E .For region-agnostic statements we simply denote a phase-space point as p.
In order to explain how we convert the results of ref. [24] to the kinematic region in eq.(3.2), let us briefly summarize how they were obtained.The discussion below holds for each power of N f , that is for each A (i) [j] , but for simplicity of the expressions we will suppress this dependence in intermediate steps.The first step is a computation of the decomposition of the bare partial amplitudes at one and two loops into a set of master integrals: where we make the dependence on the dimensional regulator ϵ explicit.This decomposition was obtained using the two-loop numerical unitarity approach [18,20,51,52].At this stage, the master integrals were replaced by their expansion in terms of Euclidean-region pentagon functions. 3Denoting the vector space spanned by relevant pentagon-function monomials in this region as ⃗ f (p E ) = {f l (p E )} l∈B , with B denoting the set of multi-indices l, we write In this decomposition, we assume that we have a pure basis of master integrals [53] at one and two loops, in which case the coefficients d k,l are rational numbers.Furthermore, we assume all integrals are normalized so that their Laurent expansion around ϵ = 0 has no poles, and stress that the set of pentagon functions is common to one and two loops.By inserting the decomposition (3.4) into eq.(3.3) for both one-and two-loop amplitudes, and then using these expressions in the definitions for the finite remainders given in eq.(2.16), we obtain a decomposition of the two-loop remainders into pentagon functions: where we have reintroduced the superscript [j] to stress that the decomposition holds for each power of N f .Obtaining the analytic form of the coefficients r (h, p) was the main result of ref. [24].Let us finish this summary by commenting on the structure of the coefficients r (i)[j] k (h, p).They can be written as where we briefly suppress the i and j indices, we recall ⃗ s = {s 12 , s 23 , s 34 , s 45 , s 15 }, and r ± k (h, ⃗ s) are rational functions of the Mandelstam variables.The analytic continuation of the r k from Euclidean to physical region is therefore trivial.For a given i and j, the r k (h, p) are not all linearly independent of each other (over Q) for different values of k, and we can define a basis of rational functions ⃗ v(h, p) = {v 1 (h, p), . ..} such that all r k (h, p) live in the Q-linear span of ⃗ v(h, p).We can then write where is a matrix of rational numbers.To understand our approach to the analytic continuation of eq.(3.7), let us begin by considering what would happen if we had chosen a different basis of pentagon functions where N is a matrix of rational numbers.It is clear that to write eq.(3.7) in this new basis of functions we do not need to determine the basis ⃗ v(h, p E ) of rational functions.The remainders would have exactly the same form, but with the matrix M E replaced by M which is just a different matrix of rational numbers.
Let us now discuss what happens with the pentagon functions under analytic continuation, which we will do through an example.Consider the simplest non-trivial pentagon function, log(−s 12 ).The choice to write it in this form follows from the fact that in the Euclidean region s 12 < 0 and the master integrals are real valued.We note that, in this region, we could also have written log(|s 12 |).Let us now continue this function to the physical region in eq.(3.2),where s 12 > 0. We find where the sign is determined by the usual Feynman i ε-prescription of the propagators.We see that the analytic continuation into the physical region introduces a new pentagon "function", the transcendental constant i π.Crucially, however, the coefficient of this function is generated by the analytic continuation of log(−s 12 ).This implies that the coefficient of i π is trivially related to the coefficient of log(−s 12 ) when analytically continuing eq.(3.4) or eq.(3.5), or more precisely that it is in the Q-linear span of ⃗ v(h, p φ ) (we recall that p φ is used to denote a point in the physical region).We assume that this is a general feature of the analytic continuation of pentagon functions: while the dimension of the vector space in eq.(3.5) might increase, the coefficients of the new pentagon functions are not independent from the coefficients of the pentagon functions in the Euclidean region. 4Even though we do not prove this statement, as we shall see below we can explicitly check that it holds for the remainders in eq.(3.7).
To analytically continue the remainders computed in ref. [24], we use the set of pentagon functions defined in ref. [41], and we denote the vector space of their monomials as ⃗ g(p φ ).These are valid precisely in the kinematic region of eq.(3.2), and are related to the set ⃗ f (p E ) used in ref. [24] by a change of basis and by analytic continuation.Given the points raised above about these two operations, this means that remainders in the physical region can be written in terms of the same bases of rational functions that were found for the Euclidean region.Hence we can write the partial remainder in the physical region as where for completeness we made the dependence on µ explicit.We stress that this representation is valid only when p 1 and p 2 are in the initial state and p 3 , p 4 and p 5 are in the final state, because that is the region of validity of the pentagon functions ⃗ g(p φ ).We discuss the generalization to other configurations in section 3.2.The M (h) can be determined from a small number of numerical evaluations of the amplitudes.To do this, we take equation eq.(3.9) as an ansatz for the form of the remainder when expressed in terms of the physical region pentagon functions ⃗ g(p φ ), with the entries of the M (i)[j] φ (h) as unknowns.We then constrain the unknowns by performing the analogous procedure of eqs.(3.4) and (3.5) at a sufficient number of phase-space points. 5Linear algebra then allows us to extract the values of In practice, we use the implementation of two-loop numerical unitarity in Caravel [44] to generate all the required numerical data.For efficiency reasons, the procedure is performed over finite fields and we construct the Q-valued M (i)[j] φ (h) via the Chinese remainder theorem and standard rational reconstruction techniques.We remark that we used at most three finite fields of cardinality O(2 32 ).Following these steps, we obtain the analytic continuation of the partial amplitudes computed in ref. [24] to the physical region defined in eq.(3.2), corresponding to the channels in eqs.(2.1) to (2.5).Crucially, this linear algebra exercise would not allow one to compute the matrices M (i)[j] φ (h) in the ansatz (3.9) if the assumption on the analytic continuation of the pentagon functions did not hold.

From Partial Amplitudes to Amplitudes
For physical applications, it is not sufficient to know a single partial amplitude for each channel.Indeed, we must be able to evaluate all the terms in the sum in eq.(2.7) for each helicity configuration.We start by summarizing this information for each of the channels listed in eqs.(2.1) to (2.5).We remind the reader of the comment on momenta conventions under eq.(2.

5).
Five gluon: There is a single channel, see eq. (2.1).The set Σ = S 5 /Z 5 contains 24 elements, and the color factor for σ 0 = {1, 2, 3, 4, 5} is For each of the 24 partial amplitudes, there are 32 different helicity configurations.Accounting for relations under charge conjugation and parity, there are 192 partial helicity amplitudes to compute.We can further reduce this set by considering the S 2 × S 3 permutations of initial state and final state labels which, crucially, still correspond to the physical region defined at the start of section 3.1.This leaves us with a generating set of 32 amplitudes.
Single quark line: There are three channels, see eq. (2.2).For the channel qh 1 p 1 q h 2 p 2 → g h 3 p 3 g h 4 p 4 g h 5 p 5 , the set Σ = {1, 2} ⊗ S 3 (3, 4, 5) contains 6 elements, and the color factor for For each of the three channels there are 6 partial amplitudes, each with 16 different helicity configurations.Accounting for relations under charge conjugation and parity, there are 144 partial helicity amplitudes to compute.Considering the S 2 ×S 3 permutation of initial state and final state labels further reduces the cardinality of the generating set to 40.
Two quark lines: There are four channels, see eq. ( 2.3), but we also include the two channels in eq. ( 2.5) in order to consider two identical quark lines.For the channel qh p 4 g h 5 p 5 , the set Σ = {{1, 2, 3, 4, 5}, {3, 4, 1, 2, 5}} contains 2 elements, and the color factor for σ 0 = {1, 2, 3, 4, 5} is For each of the six channels there are 2 partial amplitudes, each with 8 different helicity configurations.Accounting for relations under charge conjugation and parity, there are 48 partial helicity amplitudes to compute.Considering the S 2 × S 3 permutation of initial state and final state labels further reduces the cardinality of the generating set to 24.
Having classified a generating set of partial amplitudes for each channel which respect the region constraints of the pentagon functions, we can now apply an analogous procedure to that of section 3.1 to each of these partial amplitudes.Let us consider the remainder of a partial amplitude where the helicities and momenta are related by a permutation σ to those in eq.(3.9).Clearly, the master integrals in the decomposition (3.3) of the permuted amplitude are related to those of the original amplitude by the permutation σ.We remind the reader that the pentagon functions of ref. [41] form a basis of the space of transcendental functions with p 1 and p 2 initial state and p 3 , p 4 and p 5 final state.Hence we can express any permutation of the master integrals with this configuration of the momenta as where B φ is the set of labels distinguishing the elements of ⃗ g(p).It therefore follows that all permutations of the partial remainders contributing to the remainder in eq. ( 2.18) can be expressed in the same basis of pentagon functions.That is, a partial finite remainder evaluated with a permutation of the input helicities σ(h) and on a permutation of the input point σ(p) can be expressed in the form where, as in eq.(3.9), we made the dependence on µ explicit.We recall that we can easily construct the bases ⃗ v (i) [j] (σ(p), σ(h)), as the action of σ on rational functions is trivial to evaluate.The matrices (σ(h), σ) can then be determined by the same ansatz procedure described for the analytic continuation of the results of ref. [24].Altogether, this allows us to obtain all the expressions required to evaluate the squared finite remainders defined in section 2.3 for the channels in eqs.(2.1) to (2.5).

Results and Validation
In this section we describe our main results, which are of two types.The first is a collection of ancillary files containing analytic expressions for a generating set of (partial) remainders for each of the channels in eqs.(2.1) to (2.5).The second is a C++ library that allows for an efficient numerical evaluation of the finite remainders defined in eq.(3.9) for all required permutations, as well as the squared finite remainders defined in eqs.(2.22) and (2.23).Finally, we describe the checks we have performed on our results.

Analytic Results
In the ancillary files we give all necessary ingredients to assemble analytic expressions for a generating set of two-loop partial remainders R (2)[j] φ (σ(h), σ(p)), for each of the channels in eqs.(2.1) to (2.5).From this generating set, all the permutations required in eq.(2.18) can be constructed by applying (a combination of) charge conjugation, parity transformations, 6  and the S 2 × S 3 permutations of initial and final state labels. 7 The expressions we present in the ancillary files allow to explicitly construct the generating set of remainders in the decomposition of eq.(3.14).The basis of transcendental functions ⃗ g(p) is universal across all channels and powers of N f , and is given in the file anc/pentagonFunctionBasis.m, written in the notation of the Mathematica package provided with ref. [41].The basis of rational functions ⃗ v (2)[j] (σ(h), σ(p)) is channel and power-of-N f specific but trivially dependent on σ (more precisely, the action of σ on ⃗ v (2)[j] (σ(h), σ(p)) simply amounts to the permutation of the labels in the Mandelstam invariants appearing in the rational functions).For each channel and power of N f , we simply reuse the bases in the ancillary files of ref. [24] (we note again the comment below eq.(2.5) regarding the conventions to denote each channel).For completeness, we include them in the file anc/rationalBases.m in the format {channel, helicity, Nf} -> {list of rational functions in (⃗ s, tr 5 )}.
We note that one can alternatively use the form of the same bases of rational functions determined in ref. [25], which are more compactly expressed in the spinor-helicity formalism.Finally, the matrices M (2)[j] φ (σ(h), σ) for each member of the generating sets of remainders are given explicitly in anc/twoLoopProjectionMatrices.m.Each is expressed in the form {channel, helicity, NfPower, permutation} -> matrix. 6Parity transformation of the analytic expressions are performed by changing the signs of all parity-odd pentagon functions and tr5. 7One may wish to express the S2 × S3 permutations of the partial remainders from the generating set in terms of the pentagon functions evaluated on a non-permuted phase-space point.This can be achieved by following the procedure described in section 3.2.
where matrix is in Mathematica's SparseArray format.As described in section 3, these matrices constitute the main analytic result of this paper.Indeed, they were the only missing information required to obtain analytic expressions for the five-parton amplitudes in the physical region corresponding to three-jet production at hadron colliders given in eq.(3.2).

Numerical Evaluation
Having phenomenological applications in mind, we implemented a C++ library that allows for the efficient evaluation of finite remainders and squared finite remainders.In this section we comment on this numerical implementation.
Let us start by describing the quantities which can be evaluated with our C++ library [45].The code performs the numerical evaluation of one-and two-loop (partial) remainders for each channel, each helicity configuration and each power of N f .That is, we evaluate each of the R (i) [j] (σ(h), σ(p)) for i = 1, 2 and present each power of (N f /N c ) j separately.We stress that we have implemented all permutations, and not only the generating set for which we provide analytic expressions.For physical applications, one might be interested in the squared finite remainders defined in eqs.(2.22) and (2.23).Our code also outputs the numerical values for these quantities, once again presenting each power of N f /N c separately, that is the H (i) [j] in eqs.(2.24) and (2.25).
We employ PentagonFunctions++ [41] for the numerical evaluation of the pentagon functions.For the efficient evaluation of the rational functions, we generate optimized code with FORM [54,55].In any of the helicity remainders R (2)[j] (σ(h), σ(p)) that we compute numerically, the dominant contribution to the evaluation time is that of the pentagon functions.When computing H (1) and H (2) , we have thus organized our numerical code such that the pentagon functions are evaluated only once per phase-space point, independently of the number of permutations contributing to the sums in eqs.(2.22) and (2.23).We note nevertheless that, compared to e.g. the case of triphoton production [32], the evaluation of H (2) for five-gluon channels receives contributions from a large number of partial amplitudes, and therefore the time spent in the evaluation of the other components in eq.(3.14) becomes noticeable.We find that the typical double-precision evaluation time for any H (2)  is O(1s) on a single CPU, while the evaluation time for H (1) is negligible at O(10ms).
To be suitable for phenomenological applications the numerical evaluations should not only be fast but also stable across phase space.To guarantee the latter, we develop a precision-rescue strategy that detects and corrects unstable numerical evaluations.The main source of numerical instabilities is the presence of large cancellations between different contributions in eq.(3.14), when a phase-space point gets close to the zero sets of the rational functions' denominators.As we have already argued in refs.[23,24,32], these are determined by a set of 25 quantities {W i (⃗ s)}, with i = 1, . . ., 25, which are a subset of the so-called alphabet associated with the five-point two-loop master integrals (see e.g.ref. [56] for the definition of these quantities).This subset is rational in the Mandelstam variables ⃗ s defined in eq.(2.9) and all elements have the same dimensions.We find that potentially unstable points can be efficiently identified by computing a quantity that characterizes the spread of scales for a given phase-space point ⃗ s,8 We introduce a threshold κ thr , and if κ(⃗ s) > κ thr the point is considered stable.Otherwise, the point is potentially unstable and we perform a second evaluation on an infinitesimally perturbed point ⃗ s δ defined as where ε double ≃ 10 −16 is the machine epsilon of double-precision floating point numbers. 9e estimate the accuracy of potentially unstable points as We introduce another threshold value ∆ thr , and if ∆(⃗ s) < ∆ thr the point is considered stable.If ∆(⃗ s) > ∆ thr , the point ⃗ s is considered unstable and evaluated in quadruple precision.The performance of this rescue strategy and the values we choose for κ thr and ∆ thr will be discussed below.
To demonstrate the performance of our numerical C++ implementation, we employ Sherpa 2.2 [57] to sample phase-space points from a Monte-Carlo integration grid optimized on the Born matrix elements for three-jet production.We sample a sufficient number of points to draw 110k events for each of the channels gg → ggg, qg → qgg and qq → Q Qg, which are representative of all channels in eqs.(2.1) to (2.3).We use the same phase-space definition that was used in the 7 TeV ATLAS analysis [58], which we summarize here for completeness.We require three hard jets defined with the anti-k t algorithm, with the radius parameter of R = 0.4, as implemented in FastJet [59].We require the jets to be within the rapidity range |y| < 3 and that they all have transverse momenta larger than 50 GeV.In addition, the leading and subleading jets are required to have transverse momenta larger than 150 GeV and 100 GeV respectively.We then evaluate H (2) at the dynamical scale and sum the contributions of different powers of N f /N c , setting N c = 3 and N f = 5.The accuracy of double-precision evaluations is determined by comparing them to an evaluation in quadruple precision as

Validation
The main analytic results of our paper are the matrices described in section 4.1.These are obtained by building on the analytic results of ref. [24], which have undergone very stringent checks.Compared to ref. [24], our current results are written in terms of a different basis of transcendental functions [41] which have also undergone several direct checks (e.g., comparison with the functions of ref. [40]) and indirect checks [32,33,35,36].Furthermore, we have verified that we reproduce the results independently obtained in ref. [19] for the numerical evaluation of five-gluon amplitudes in the physical region defined in eq.(3.2).Finally, we compared our results for the function H gg→ggg in the purely gluonic channel against an independent computation [60] and found agreement.
We have also performed several consistency checks on our numerical results.Since we obtain all the helicity partial remainders required for the functions H (1) and H (2) from sums over the permutations of the generating sets of remainders introduced in section 3.2, it is important to check that these operations are correctly implemented in the C++ library we have developed.We have made three types of consistency checks to validate this procedure.
Symmetries.We have checked that the remainders and squared remainders satisfy all expected charge, parity conjugation and permutation symmetries.
Validation of H (1) .We have compared the one-loop squared remainders H (1) of all the channels in eqs.(2.1) to (2.5) against the evaluations by the BlackHat library [61] and found full agreement. 10The color-dressed remainders R M in eq.(2.18) are constructed from generating sets of partial remainders, as described in section 3.2, and then are directly fed into the definition of H (1) , see eq. (2.22).The procedure for assembling the squared remainders is independent of the number of loops, hence this check validates the corresponding assembly of the contributions for the calculation of H (2) .Scale dependence of remainders.Given that we evaluate the finite remainders directly, the correct pole structure of the underlying scattering amplitudes is satisfied by construction.We note nevertheless that finite remainders retain a memory of the poles of the amplitudes and of their explicit (scheme-dependent) definition.This information is encoded in the scale dependence of the finite remainders.Quite generally, we find that ∂ R (1) where we denote by λ the power of g 0 s in the tree-level amplitude (for five-parton amplitudes λ = 3, see eq. (2.6)) and the β i are the coefficients of the QCD β function explicitly given in eq.(2.13).H [n] is given explicitly in eqs.(A.6) and (A.7) and controls the collinear singularities of scattering amplitudes.We have verified that the scale dependence of the remainders obtained from our numerical implementation is described by eq.(4.6).Furthermore, we have independently checked that, starting from this scale dependence, we reproduce the scale dependence of the functions H (1) and H (2) , which again validates the correct assembly of these functions.

Conclusion
In this paper we present the NNLO double-virtual contributions to three-jet production at hadron colliders at leading color.While all required amplitudes were already available in the literature in analytic form, they were only available in the non-physical Euclidean region, and as such could not be directly used for phenomenological applications.Aside from analytically continuing those expressions to the physical region, we also provide a C++ library that allows to efficiently evaluate these contributions and obtain stable results across phase space.
The analytic continuation of the results of ref. [24] to the physical region corresponding to three-jet production is performed by rewriting them in terms of the pentagon functions defined in ref. [41].By noticing that this rewriting amounts to computing a matrix of rational numbers, we obtained the physical-region expressions from a few numerical evaluations of the amplitudes (57 in the most complicated case).These evaluations were performed with the implementation of the two-loop numerical unitarity approach in Caravel [44].We provide analytic expressions for a generating set of partial remainders.For all channels in eqs.(2.1) to (2.5), these allow one to evaluate the color-dressed remainders defined in eq.(2.18) as well as the squared remainders H (1) and H (2) defined in eqs.(2.24) and (2.25) in the physical region.
Our driving motivation was to provide results that are ready to be used for phenomenological applications.To this end, we developed a C++ library which efficiently evaluates individual partial remainders as well as the one-and two-loop squared remainders.The numerical evaluations can be performed in double or quadruple precision.We find that double precision is sufficient for the vast majority of phase-space points, and have implemented a rescue system which detects unstable points and reevaluates them in quadruple precision.With the rescue system enabled, we find very good overall numerical stability and our evaluation time for the two-loop squared remainders is O(1s).This shows that our results and the numerical library we provide are ready to be used in phenomenological studies. Channel H (1) [1]  H (2)[0] H (2) [1]  H (2) [2]   gg

B Reference Evaluations of Squared Finite Remainders
In order to facilitate the comparison with our results, we present in table 2 reference values for the evaluation of squared finite remainders at each power of N f , as defined in eqs.The results are obtained with our C++ library.

Figure 1 :
Figure1: Characteristic Feynman diagrams which contribute to A (2)[j] for representative channels with no external quark lines, a single quark line and two external quark lines, for j = 0 (left column), j = 1 (middle column) and j = 2 (right column).