Soft triple-real radiation for Higgs production at N3LO

We present the two first terms in the threshold expansion of Higgs production partonic cross-sections at hadron colliders for processes with three partons in the final state. These are contributions to the inclusive Higgs cross-section in gluon fusion at N3LO. We have developed a new technique for the expansion of the squared matrix-elements around the soft limit and for the reduction of the required phase-space integrals to only ten single-scale master integrals. We compute the master integrals building upon modern techniques for the integration of multidimensional integrals in dimensional regularization. Our results constitute an important step towards a systematic computation of the Higgs boson cross-section as an expansion around the threshold limit.


Introduction
The first years of experiments at the Large Hadron Collider (LHC) have probed decisively the physics of electroweak symmetry breaking, demonstrating a great understanding of the theory at this energy scale. The success of particle theory must be attributed not only to the Standard Model (SM), but also to remarkable progress in perturbative calculations. Results in this field of research are outstanding and have resulted in very accurate comparisons with data. Indeed, partonic cross-sections for a large variety of processes are now routinely computed at next-to-leading order (NLO) in the strong coupling constant α s 1 .
In addition, modern Monte-Carlo programs have reached a high level of sophistication with the development of methods to match and merge parton showers with fixed-order calculations 2 . Finally, basic processes at the LHC with a small final-state multiplicity can now be computed fully differentially at next-to-next-to-leading order (NNLO) 3 .
A classic example where perturbative calculations have been crucial for the interpretation of LHC data is the production of a Higgs boson. Perturbative corrections for Higgs boson signal processes are generically large, especially so for production via gluon fusion [49][50][51][52][53][54][55]. For this process, the total and fully differential cross-sections are known through NNLO with a typical uncertainty of about 10% due to variations of the factorization and renormalization scales [46,56]. Indeed, without the knowledge of the NNLO corrections, predictions at LO or NLO would be assigned theory uncertainties which are larger than the already achieved experimental uncertainties. With future LHC data, a comparison with theory at a level of precision of a few percent will be revealing the fine details of the mechanism of electroweak symmetry breaking and potentially uncovering new fundamental laws of physics.
Besides being of importance to Higgs physics, NNLO predictions have promoted the Drell-Yan process to a standard candle for studies at hadron colliders, where observables in this process are computed with a typical precision of better than 2% as indicated by scale variations [57,58]. Measurements of the mass of the W boson and the weak mixing angle from Drell-Yan data are one of the most stringent constraints on the Standard Model. Furthermore, Drell-Yan production data is an essential input for the extraction of parton distribution functions from hadron collider processes. Last but not least, the clean detector signatures of Drell-Yan events, combined with the excellent theoretical predictions of their rates, allow for a precise determination of the luminosity. Indeed, luminosity determination from Drell-Yan is relatively insensitive to high pile-up conditions, a fact which will be even more important during the high-luminosity phase of the LHC. We believe that it is possible to develop efficient methods for computing Higgs boson and Drell-Yan hadroproduction cross-sections at the next perturbative order, "N 3 LO". The benefits of such a program are potentially very important.
In Higgs production via gluon fusion the uncertainty due to scale variations at N 3 LO is anticipated to be half of what is found at NNLO [59]. Scale variation estimates can-not be derived from first principles, but entail a certain degree of subjectiveness. With the knowledge of more orders in the perturbative expansion, the remaining perturbative uncertainty can be estimated more reliably, not only from scale variations but also from the progression of the series. Improvements of the Higgs boson production cross-section uncertainty will propagate into the determination of the couplings of the Higgs boson.
Although existing NNLO calculations for the Drell-Yan process indicate that the crosssection is already known very precisely, an explicit N 3 LO computation will educate us further by challenging or verifying the traditional prescriptions for uncertainty estimates. Furthermore, high-luminosity runs at the LHC require the triggering of leptons at a higher transverse momentum, resulting in a deterioration of the NNLO scale uncertainty [45]. An excellent theoretical precision can then be recovered by including the N 3 LO perturbative corrections.
So far, there has been no N 3 LO calculation for hadron collider processes performed in the literature. In this paper, we attempt a first step towards computing inclusive hadroproduction cross-sections at this perturbative order. We focus on one of the most complicated N 3 LO contributions to the inclusive Higgs boson cross-section in the gluon fusion channel, namely, the contributions coming from partonic cross-sections for tree-level radiative processes with three additional partons emitted in the final state. A direct integration of the corresponding matrix-elements over phase space is however tedious. To simplify the problem, we perform a threshold expansion and devise a method to compute the coefficients of the expansion analytically. One of the main results of our paper is the computation of the first two terms of the threshold expansion for triple real corrections to inclusive Higgs production.
Our expansion method builds upon reverse unitarity [55,58,[60][61][62], a technique developed for the calculation of the inclusive Higgs cross-section and the rapidity distribution of electroweak gauge bosons at NNLO. The main idea is that phase-space integrals of matrixelements are dual to loop integrals in their algebraic properties: recurrence identities in the powers of propagators and the number of dimensions as well as differential equations satisfied by loop integrals also apply to phase-space integrals. This is achieved by associating on-shell conditions or other phase-space constraints in the form of delta functions with differences of otherwise identical propagators with opposite infinitesimal imaginary parts, in the spirit of Cutkosky's rules. In this paper, we observe and exploit that the duality rules of reverse unitarity can be expanded in kinematic parameters.
After performing a threshold expansion of loop integrals dual to the cross-sections for Higgs plus three partons processes, we apply automated reduction algorithms to reduce the coefficients of the expansion to master integrals. We have performed the reduction using existing public programs and programs developed specifically for the purpose of this computation.
The reduction yields ten master integrals for the first two coefficients in the threshold expansion. The master integrals themselves are not specific to Higgs production and will appear in the threshold expansion of other hadroproduction processes, such as Drell-Yan. For processes with a single colorless final state our set of master integrals is complete and our expressions are universal for the partonic cross-sections with the same initial state at leading order in the threshold expansion.
For the computation of the soft master integrals we have employed a combination of methods. The master integrals are evaluated by using a parametrization of the phase space in terms of the energies and angles of the momenta of the soft partons. When needed, we perform the phase-space integrals by introducing Mellin-Barnes representations. We evaluate some of the Mellin-Barnes master integrals in a closed form as hypergeometric functions or directly as a Laurent-expansion in the dimensional regulator ǫ with infinite nested sums as coefficients, which naturally evaluate to multiple zeta values in all cases. Some master integrals cannot be evaluated easily from their Mellin-Barnes representations. We have developed a procedure to turn Mellin-Barnes integrals into integrals over positive real parameters, which are easy to expand in ǫ whenever the integral is finite. The resulting parametric integrals are then evaluated by integrating out the integration variables one-byone in terms of multiple polylogarithms. The last integration then again produces multiple zeta values in a natural way.
We have restricted our computation of the partonic cross-sections to the first two terms of the threshold expansion. There are two obvious ways to extend this work in the future. We could continue with the computation of subleading terms in the threshold expansion. This requires a more intensive but not prohibitive computer algebra and the evaluation of soft master integrals from topologies which contribute only to higher orders in the soft expansion. The soft expansion is rather fast converging for Higgs production but rather slow for Drell-Yan, as has been noticed at NNLO. A second possibility is to perform the phase-space integrations for arbitrary kinematics, reducing them to a different set of master integrals. The soft master integrals which we present in this article can serve as boundary conditions for solving the differential equations satisfied by the master integrals in arbitrary kinematics.
This article is organized as follows: In Section 2 we present our general strategy based on reverse-unitarity to perform the threshold expansion for real-emission cross sections in terms of soft integrals, and we illustrate the method on some simple examples in Section 3. In Sections 4 and 5 we study some properties of soft integrals in general, and we show that there is an easy and canonical way to derive dimensional recurrence relations and Mellin-Barnes integral representations for generic soft integrals. The soft master integrals contributing to the first two terms in the threshold expansion of the leading-order partonic cross sections for p p → H + 3 partons are discussed in Section 6, and the corresponding results for the cross sections are presented in Section 7. Section 8 contains technical details about the computation of the soft master integrals, and in Section 9 we draw our conclusions and give an outlook for future work. The paper contains several appendices discussing phase parametrisations and a generic formula for the phase-space volume for H + n partons, as well as a new method to derive a parametric integral representation from a Mellin-Barnes integral and a description of an algorithmic way to perform the analytic integration of certain classes of parametric integrals.

Reverse unitarity, threshold expansion and soft integrals
We consider the production of a Higgs boson in association with j = 3 . . . N massless partons in the final state from two massless partons i = 1, 2 in the initial state, 1 + 2 → H + 3 + . . . + N. (2.1) The inclusive cross-section for this process in dimensional regularization is given by a phase-space integral over the momenta q j of the final-state partons, σ = dΦ N −1 (q H , q 3 , . . . , q N ; M 2 ; s; D) |A| 2 ({q j }, q 1 , q 2 ; D).
We work in D = 4 − 2ǫ dimensions and denote by q 1 , q 2 the momenta of the initial-state partons and we also use the shorthand notation q 12 = q 1 + q 2 , q 345 = q 3 + q 4 + q 5 , etc.
In the following we will often drop the functional dependence on the dimension D for clarity. The mass of the Higgs boson is denoted by M , and we denote the (squared) center-of-mass energy by s = (q 1 + q 2 ) 2 . |A| 2 represents the squared matrix-element multiplied with the appropriate flux and symmetry factors. The D-dimensional phase space measure is given by dΦ N −1 (q H , q 3 , . . . , q N ; M 2 ; s; D) with δ + (q 2 − m 2 ) = δ(q 2 − m 2 )Θ(q 0 ). Integrating out the momentum of the Higgs boson, we can rewrite eq. (2.2) as ..N − q 12 ] 2 − M 2 |A| 2 ({q j }, q 1 , q 2 ). (2.4) In this article, we restrict ourselves to the case of real-radiation matrix-elements without virtual corrections. We introduce the variables z = M 2 s andz = 1 − z . (2.5) We now rescale the momenta of all the partons, which captures the scaling of the partonic momenta in the final-state. We emphasize that this is not, as yet, an approximation, but rather a convenient change of integration variables which captures the correct asymptotic behavior at threshold asz → 0. In the following we assume s = 1, and we find Note that the full s-dependence can easily be recovered from dimensional analysis. The squared matrix-element |A| 2 consists of a rapidly growing number of terms with N , yielding a correspondingly large number of phase-space integrals. The method of reverse unitarity, developed in refs. [55,58,[60][61][62], allows the reduction of phase-space integrals to a basis of fewer master integrals by establishing a duality of phase-space and loop integrals, where the latter are amenable to algebraic methods [63,64] based on integration by parts [65,66].
Following the reverse unitarity approach, on-shell and other phase-space constraints are dual to propagators: "Cut" propagators can be differentiated in a similar way to ordinary propagators with respect to their momenta, leading to identical integration-by-parts (IBP) identities for phase-space integrals as for their dual loop integrals. Solving the system of IBP identities for phase-space integrals proceeds in the same way as for loop integrals, with the exception that for cut-propagators we can use the simplifying identity: According to the reverse unitarity method, we find a dual forward scattering loopamplitude with N − 1 cut-propagators for the real radiation contribution of eq. (2.7), namely, In this article, we take one further step and expand cut-propagators and the squared matrix-elements around z = 1, (2. 13) In this approximation, the cross section can be expanded into a power series inz, 14) The coefficients of the power series are given by where dΦ S N −1 denotes the "soft" phase space measure (2.16) The integrals which emerge after thez expansion depend trivially on one dimensionful parameter p 2 12 = s. If we put s = 1, the integrals are number integrals whose only functional dependence is through the space-time dimension D = 4 − 2ǫ. We will refer to such integrals as soft (phase space) integrals, and they are the main subject of this paper. We note that, apart from the cut Higgs boson propagator, the integrands of soft phase space integrals are homogeneous functions under a simultaneous rescaling of the final-state momenta. In addition, a soft integral can be reduced to a set of "soft" master integrals using IBP identities by exploiting the duality to loop integrals via reverse-unitarity. We will illustrate this property in the next section where we check our method on several examples.

Validation of the method and examples
In this section, we study the validity of the method described in the previous section at NLO and NNLO -two perturbative orders that are well studied in the literature and so we can compare our results readily with known results. In particular, we show that our method reproduces the correct results for the leading behavior of NLO and NNLO real emission amplitudes in the soft limit, as well as for the subleading terms in the expansion of the phase space volume up to N 3 LO and for a non-trivial double real emission master integral at NNLO.
At NLO, all phase space integrals that contribute to the real emission amplitude in general kinematics can be reduced to the phase space volume for H + 1 parton, As there is only one master integral which is a monomial inz, our method trivially gives the correct answer at NLO. At NNLO all double real emission phase space integral can be reduced in general kinematics to a linear combination of 18 master integrals [55]. The leading contribution of all master integrals in the soft limit to all orders in ǫ was computed in ref. [67], and it was observed that in this limit 17 master integrals are proportional to the soft limit of the phase space volume for H + 2 partons, where we defined More precisely, it was shown in ref. [67] that if X S i (z; ǫ) denotes the leading term in the soft limit of the double real emission master integrals, then we can write 4 where S i (z; ǫ) are monomials inz and rational functions of ǫ. Using the method described in the previous section, we can easily explain the structure of eq. (3.4). Indeed, we observe that in the soft limit all the double real emission phase space integrals can be reduced to only two master integrals. In particular, the IBP identities in the soft limit allow us to express all but one of the X S i in terms of the phase space volume, and the coefficients appearing in the reduction are precisely the functions S i . In other words, in the soft limit all double real emission phase space integrals can be reduced to linear combinations of the following two soft master integrals Our method thus provides the correct leading soft behavior of the double real emission contribution at NNLO. We emphasize that all the diagrams in this paper represent soft phase space integrals, i.e., all the diagrams represent integrals with respect to the soft phase space measure of eq. (2.16). In addition, the invariants appearing in the integrands of the soft integrals are defined with respect to the rescaled momenta defined in eq. (2.6), Our method does not only allow us to compute the leading soft behavior, but we can consistently expand around the soft limitz = 0. In the following we show that we can correctly reproduce the first few terms in the soft expansion of double and triple emission phase space volumes, as well as for the NNLO master integral X 18 of refs. [55,67,68].
Let us start with the phase space volume for H + 2 partons in the limit where the two partons are soft. On the one hand, from eq. (3.2) we immediately see that Φ 3 admits the expansion On the other hand, using eq. (2.13) we obtain the diagrammatic expansion where the dashed lines indicate numerator factors and dots represent additional powers of the propagators or the numerators. The diagrams appearing in eq. (3.9) are in one-toone correspondence with the terms in the expansion (3.8). Indeed, IBP reduction of the integrals in eq. (3.9) reveals 11) in perfect agreement with eq. (3.8). We checked explicitly that our method reproduces correctly the first ten terms of the soft expansion of the phase space volume for H + 2 partons.
As a second example we derive the subleading terms in the soft expansion of the double real emission master integral X 18 . Unlike for the phase space volume, no result is known for X 18 valid to all orders in ǫ in general kinematics, but the integral was evaluated explicitly up to O(ǫ) in terms of harmonic polylogarithms [69] in ref. [55,67,68]. We can thus compare the result of our method order by order in ǫ to the expansion of the harmonic polylogarithms around z = 1. Using eq. (2.13) we obtain IBP reduction of the diagrams appearing in the subleading terms gives (3.14) We checked that using these identities we can reproduce correctly the first five terms in the soft expansion of X 18 . The aim of this paper is to compute the leading terms in the soft expansion of the triple real emission amplitude for inclusive Higgs production. In order to test our method at N 3 LO, we verified that we can reproduce the correct soft expansion of the phase space volume for H + 3 partons. The phase space volume for H + 3 partons in general kinematics can be written in the form (see Appendix B) where we defined Using our method, we obtain the following diagrammatic expansion All the diagrams in the expansion can be reduced to the soft phase space volume, as expected, To summarize, our method provides a systematic way to perform the threshold expansion of phase space integrals for the production of a heavy colorless state. Every term in the expansion corresponds to a soft integral, as defined in Section 2, which can be reduced to a small set of soft master integral using IBP reduction. In the next two sections we study some additional properties of soft integrals in general, before applying our method to compute the threshold expansion of the triple real emission contribution to inclusive Higgs production.

Relations between phase space integrals in different dimensions
It is well known that loop integrals in different space-time dimensions are related by socalled dimensional shift identities [70]. After reduction to master integrals, the dimensional shift identities reduce to recurrence relations in the space-time dimension D for the master integrals themselves [71][72][73][74][75]. Reverse unitarity rules show that the same dimensional shift identities hold for the dual phase-space integrals (see also ref. [76]).
In this section we present an easy way to derive the dimensional shift identities for phase space integrals. We stress that the results of this section are generic and apply to phase space integrals in general kinematics. To start, let us consider a phase space integral in D dimensions, where we explicitly indicate the dependence on the space-time dimension D. The integrand f can be written as a product, where the P l are polynomials in the rescaled kinematic invariants s ij , raised to some power ν l . In Appendix A we show that the phase space measure dΦ N −1 (D) for parton+parton → H + (N − 2) partons can be parametrized solely in terms of kinematic invariants, where we have defined the normalisation factor , , (4.4) and the Gram determinant G N Obviously, the D-dependent constants factor out of the integral, and so the actual integral depends on D only through the exponent of the Gram determinant. It is then easy to see that the phase space measures in shifted dimensions are related by We can thus express a phase space integral in D + 2 dimensions as For a given set of polynomials {P l }, the form of the integrand f depends only on the exponents {ν i }. One therefore introduces the operators I + i and I − i acting on f as  We thus obtain the following compact formula relating phase space integrals in different dimensions, Every term in the polynomial can be evaluated according to the action of the I − i operators, yielding a superposition of modified integrals in D dimensions. By applying this method to a master integral, we can express the master integral in D + 2 dimensions as a linear combination of integrals in D dimensions. Using IBP identities, we can reduce the integrals in D dimensions to master integrals and thus we find a relation between the master integral in D + 2 dimensions and D dimensions. This dimensional recurrence relation can formally be written as with coefficients c ij (D) that are determined from the IBP reduction.

From soft integrals to angular integrals and Mellin-Barnes integrals
In this section we show that there is a canonical way to derive a Mellin-Barnes (MB) representation for the soft integrals that appear in the threshold expansion of phase space integrals. The soft integrals we need to consider have the form where f ({p j }; p 1 , p 2 ) is a ratio of products of multi-particle invariants that is homogeneous under a simultaneous rescaling of the final-state momenta, i.e., for some a. Note that this is true only if the (cut) Higgs boson propagator it not raised to an additional power. In that case, however, we can always reduce the integrand to a homogeneous function using IBP identities. Next, we note that there is a subclass of soft integrals that have an additional property: they are homogeneous with respect to individual rescalings of the final-state momenta, i.e., for some a j . This subclass of soft integrals is precisely the one where the integrand consists of products of power of two-particle invariants, where the index k runs over all the two-particle invariants appearing in f . Every soft integral can be converted into an integral of this type, to the price of introducing additional MB integrations. Indeed, if we write every multi-particle invariant as a sum of two-particle invariants, then we can convert sums into product by using the usual formula where the contour separates the poles at z = n from those in z = λ − n, n ∈ N. Without loss of generality we can thus assume that our soft integral is homogeneous with respect to individual rescalings of the final state momenta. If we concentrate on soft integral that satisfy eq. (5.3), it is natural to choose a parametrization of the soft phase space that makes the homogeneity explicit. One possible parametrization with this property is the so-called 'energies and angles' parametrization, where the final-state momenta are parametrized as 5 The E i parametrize the energies of the final-state partons and β i is the D-velocity in the direction p i . In this parametrization the phase space measure for each final-state parton takes the form where dΩ is the measure on the unit sphere parametrizing the solid angle of particle i. Furthermore, the on-shell condition for the Higgs boson, coming from the cut propagator in eq. (2.16), can be rewritten as Thus, the soft phase space measure can be parametrized as Using this parametrization, we see that every soft integral with an integrand of the form (5.4) can be written as (5.10) Both the measure and the integrand can be written in a factorized form, and so we can integrate out the energies in terms of a generalized Beta function, Hence, the only non-trivial integration is a multiple angular integration over the solid angles of the final-state partons. Angular integrals can be written in the general form In ref. [77] it was shown that such integrals fall into a class of generalized hypergeometric functions known as H functions, and an MB representation for the most general angular integral of this type was derived. We have thus a general recipe to derive MB representations for generic soft integrals. Although the previous technique allows us to derive a multifold MB representation for every soft integral we need to consider, it can sometimes be useful to insert, if existent, explicit closed expressions for the angular integrals. Indeed, for small values of m the integrals (5.12) are very simple and can be evaluated in closed form. In the following we briefly review some results for angular integrals which will be useful in our case.
The case m = 0 corresponds to the volume of the solid angle Note that this integral is sufficient to compute the soft phase space volume, which corresponds to putting m = 0 in eq. (5.10): we simply obtain a factor Ω D−1 for each final-state parton. Thus, we obtain , (5.14) in agreement with the results of Appendix B.
As we are only interested in massless momenta, β 2 j = 0, Lorentz invariance implies that the angular integral with one propagator must evaluate to a constant. Indeed, we have For angular integrals with two massless propagators one obtains [78] Ω To our knowledge, there are no closed formulas for angular integrals with three or more massless propagators. We will however need later on the angular integral with three massless propagators, which admits the MB representation [77], with z ij = z i + z j , and where the contours separate the poles coming from Γ functions of the form Γ(. . . − z ij ) from those coming from Γ(. . . + z ij ).
6. Triple real emission phase space integrals in the soft limit

Triple real soft master integrals for Higgs production
After having studied some properties of soft integral in the previous sections, we will use the technology developed in the previous sections to compute the threshold expansion of the leading-order cross sections for H plus five partons. More details about the construction of the amplitude in this limit will be given in Section 7. Here it suffices to say that we have computed the squared amplitude and we have checked that in the limit where we only keep the first two terms in the threshold expansion, all the phase space integrals can be reduced to linear combinations of the following ten soft master integrals, We have normalized all the integrals to the soft phase space volume for H + 3g defined in eq. (3.16). In the remainder of this section we give the dimensional recurrence relations satisfied by the master integrals and present the analytic results for each master integral as a Laurent expansion in the dimensional regulator ǫ. Technical details about how to compute the master integrals analytically will be given in Section 8.

Dimensional recurrence relations
Using the technique described in Section 4, we can derive dimensional recurrence relations for all the master integrals defined in the previous section. The knowledge of these recurrence relations provides us with a strong check on our results. In addition, it turns out that the master integral F 9 (D) is easier to compute in D = 6 − 2ǫ dimensions, where it is finite, and the dimensional recurrence relations allow us to relate the six-dimensional and four-dimensional results in an easy way.
The recurrence relation for the soft phase space volume is trivial to obtain from the recurrence relation for the Γ function, As we have defined all our master integrals relative to the phase space volume Φ S 4 , we can simplify their recurrence relations by factoring out the above result. We therefore define the ratio , (6.12) where N was defined in eq. (4.4). We give the results for the remaining master integrals relative to R. The dimensional recurrence relations for the non-trivial master integrals are 14) F 6 (D) , 20)

Analytic results for the soft master integrals
In this section we present the analytical results for the master integrals contributing to hadronic Higgs production in the leading and next-to-leading soft approximation. As the explicit evaluation of the master integrals is rather long and technical, we delay all details about the computation to Section 8 and only summarize the results at this point. The first master integral, the soft phase space volume, was already given in eq. (3.16) and will not be repeated here. All the remaining master integrals have been evaluated as a Laurent series in the dimensional regulator up to terms involving zeta values of weight at most six.
We have checked that our results agree numerically with the MB integral representation for soft integrals derived in Section 5. In addition, the results satisfy the dimensional recurrence relations for the master integrals given in the previous section (integrals in the shifted dimension have been evaluated numerically using the MB representation). Finally, we make an intriguing observation in our results: if we express all the zeta values up to weight six in the basis {ζ 2 , ζ 3 , ζ 4 , ζ 2 ζ 3 , ζ 5 , ζ 2 3 , ζ 6 }, the coefficients in front of the values are integers in all cases. We note that this statement is only true in the specific basis of zeta values that we chose. The results for the master integrals are listed in the rest of this section. (6.28) 7. The threshold approximation to the partonic 2 → H + 3 parton crosssection

General setup
The production of a Higgs boson in the collision of two hadrons h 1 , h 2 is dominated by QCD processes. The hadronic cross-section is related to the partonic cross-section by the general factorisation formula Here f h i (x) are the parton-distribution functions for the parton i inside of the hadron h, σ i+j→H is the partonic cross-section and S is the square of the total centre-of-mass energy of the hadronic system. The centre-of-mass energy squared of the partonic system is consequently given by s = x i x j S. The partonic cross-section is expanded perturbatively in the strong coupling constant α S . The leading QCD contribution arises in the SM via a top-quark loop at O(α 2 S ). We consider the Higgs boson to be relatively light compared to the top quark. This justifies to work in the limit of an infinite top quark mass and consider N f light quarks. We describe the interaction of the Higgs boson with gluons by introducing the effective Lagrangian Here G a µν denotes the gluon field strength tensor and H is the Higgs field. The Wilson coefficient c H can be found, e.g. in ref. [79][80][81]. The next-to-next-to-leading (NNLO) order correction to the inclusive Higgs cross-section has been computed in the past by employing this effective theory. In this work we present a part of the next term in in the perturbative series (N 3 LO) in the effective theory.
At every order in perturbative QCD, the cross-section receives contributions from various real and virtual radiation processes. We consider only tree-level processes with three real-emission partons in the final state. A direct integration over the phase-space of the corresponding matrix elements is challenging. As we have discussed in the introduction, we pave the way towards the full computation by performing an expansion of the phasespace integrals around the kinematic limit where the Higgs boson is produced at threshold, s ∼ M 2 . In this limit, all partons emitted in the final state are soft and their momenta vanish asz → 0. We expand the cross section in the small parameterz defined in (2.5). In the following, we shall present the leading and subleading terms in the threshold expansion.
Although we considered in our calculation the Higgs boson H as a final state, we would like to stress the universality of our result for the leading term in the soft expansion for any other colorless final state produced by gluons in the initial state [82]. The subleading term in the soft expansion is no longer universal.

Calculation
To obtain the real-emission cross-section we generate Feynman diagrams using QGRAF [83] and compute squared matrix-elements using programs based on GiNaC [84] or FORM [85] and our own C++ code for color and spin algebra. We perform our calculation in Feynman gauge with D gluon polarizations in order to maintain a simple structure for the denominators of the squared amplitudes. To recover the result for physical gluon polarizations, we add matrix-elements with Faddeev-Popov ghosts as external states. We compare our results with a set of numeric cross-sections for different phase-space points obtained with MadGraph [86] and find perfect agreement.
In a next step we use reverse-unitarity and interpret the phase-space integrals as threeloop integrals as described in Section 2 and we expand them into a Laurent series inz. The leading and next-to-leading terms in this series are then referred to as the soft and next-to-soft limit of the real emission cross-section. We then derive integration-by-parts (IBP) identities for the scalar soft phase-space integrals. We reduce all integrals to a set of 10 master integrals by employing the Laporta algorithm [63]. We implemented this algorithm in a C++ code that was developed by us specifically for this project, as well as the program AIR [87]. The calculation of the remaining 10 integrals is discussed in Section 8.

Leading soft contribution
In this section we present the leading soft contributions to the i j → H + X real-emission cross-sections, where X where represents the three final-state partons. To obtain a MSrenormalised quantity we redefine the strong coupling constant where γ E = Γ ′ (1) is the Euler-Mascheroni constant and µ R is the renormalisation scale. In our calculation we consider N f massless quark flavours. We include a flux factor of 1 2s and average over incoming colours and spins, leading to a factor of 2Nc for (anti-) quarks. We abbreviate C F = N c 2 −1 2Nc and C A = N c . In the case that the phase-space integration contains n identical particles we include a symmetry factor of 1 n! . For cross-sections with 4 (anti-) quarks we distinguish between the case of only one or two quark flavours. The two-flavour case is indicated by explicitly labeling the (anti-) quarks q i orq i . If the (anti-) quarks do not carry any label they are considered all to be of the same flavour. We present our results first in terms of the soft master integrals of Section 6. The resulting expressions are valid to all orders in the dimensional regulator ǫ. Next, we insert the results for the soft master integrals and write the cross-section as a Laurent series in ǫ up to terms containing zeta values of weight six, keeping a factor corresponding to the volume of the soft phase-space Φ S 4 (ǫ), eq. (3.16), unexpanded.
The cross-section containing only gluons is given by Note that the only colour coefficient of the tree-level matrix element for a Higgs plus five gluons is C 4 A C F (see, e.g., ref. [88]). Obviously this fact is left unchanged by the soft limit.
In the case that two final-state partons are a quark anti-quark pair we sum over all possible quark flavours and obtain a factor N f . The corresponding cross-section is then given by In the soft limit only the processes with gluons in the initial state contribute. All other contributions containing at least one (anti-) quark in the initial state vanish in this limit. This can be understood in terms of soft factorisation: there is no born-level process with a massless initial state fermion in Higgs production.

Next-to-leading soft contribution
In this section we present our results for the next-to-leading term in the soft expansion of the real-emission cross sections. Only partonic subprocesses that have at least one gluon in the initial state give a non-zero contribution at this order. The results are

Analytic computation of the master integrals
In this section we present some details on how to evaluate analytically all the master integrals F i defined in Section 6. The analytic result for the soft phase space volume can easily be obtained from the expression for the phase space volume in general kinematics and will not be discussed here (see Appendix B for the derivation).
In general, our strategy is to follow the steps outlined in Section 5: we use the 'energies and angles' parametrization to obtain a representation for the integral where the energies and angles appear in a factorized form. This may require the introduction of MB integrations in order to factorize sums of invariants in a denominator. We then integrate out the energies and angles to obtain a multifold MB representation for each master integral. Whenever we are able to do so, we evaluate the remaining MB integrals to all orders in ǫ in terms of hypergeometric functions that can easily be expanded into a Laurent series ǫ using the HypExp package [89]. In those cases where we did no manage to perform the MB integral in closed form for finite values of ǫ, we only compute the Laurent expansion of the integral around ǫ = 0, e.g., by resolving singularities in ǫ and summing up harmonic sums or by converting the MB integral to a parametric integral which can be computed more easily. Details on how to perform these steps for the different master integrals will be given in the rest of this section.

The master integral F 3
In this section we compute the master integral F 3 defined by The integrand only involves two-particle invariants, so following the discussion in Section 5 we can immediately insert the energies and angles parametrization and integrate out the energies in terms of Γ functions. The remaining angular integrations can easily be carried out using the formulas given in Section 5. Note that in the present case all the angular integrals can be performed in closed form, so there is no need to introduce MB representations for the angular integrals. We obtain The remaining integral can be brought into a more standard form by the change of variable cos θ 3 = 2y − 1, The integral over y is now easily performed using the recursive definition of the hypergeometric function We immediately get The 3 F 2 function can be expanded to the desired order in ǫ using the HypExp package [89], and we arrive immediately at the Laurent series of eq. (6.23).

The master integral F 2
The integrand of the master integral F 2 involves a sum of two-particle invariants in the denominator. We replace the sum by a product to the price of introducing a MB integration via eq. (5.5), The phase space integral is now in the form (5.10), and so we can introduce the energies and angles parametrization and integrate out all the energy and the angular variables. This results in the following two-fold integral representation for F 2 , which is of mixed MB and Euler-type, The Euler integral over y could immediately be performed in terms of a 3 F 2 function, but after that we still need to integrate over the MB parameter z 1 . We therefore prefer not to perform the integration over y, but we rather insert the MB representation for the hypergeometric function in the integrand .
The integral over y evaluates to a Beta function, and we are left with the following twodimensional MB integral The integral over z 1 is easily performed using Barnes' first lemma, and the remaining onefold MB integral can immediately be recognized as a 3 F 2 function (see eq. (8.8)). We finally obtain the following result for the master integral F 2 , in agreement with eq. (6.22),

The master integral F 7
The integrand of F 7 only contains two-particle invariants, We can therefore immediately integrate out the energy and the angular variables. We obtain In order to perform the integral over y, we introduce an MB representation for each 2 F 1 function in the integrand and perform the y-integration. This leaves us with the following two-fold MB representation for F 7 , To proceed, we notice that one of the two integrations evaluates to a 2 F 1 , which can be reduced to Γ functions using Gauss' identity, (8.14) Inserting this result into the two-fold MB integral, we immediately see that the remaining MB integral evaluates to a 3 F 2 function, and we get

The master integral F 4
The integral F 4 is defined by The integrand only contains two-particle invariants and we can immediately integrate out the energy and angular variables in the usual way. We obtain a one-fold MB representation, Closing the integration contour to the right and summing up residues at z 1 = n and z 1 = −1 − ǫ + n, n ∈ N × , we immediately see that F 4 can be expressed as a combination of hypergeometric functions,

The master integral F 6
The integral F 6 contains two sums in the denominator, which we can replace by products to the price of introducing two MB integrations, We then proceed in the by now familiar way and integrate out the energies and the angles, and we arrive at the following two-fold MB representation for F 6 , Unlike in the previous cases, we were not able to reduce this integral for generic ǫ to simple hypergeometric functions. We therefore only compute the Laurent expansion of the integral. We proceed in the standard way: we apply the packages MB [90], MBresolve [91] and barnesroutines [92] to resolve singularities in ǫ and to expand the resulting integrals under the integration sign and apply Barnes' lemmas in an automated way. The resulting MB integrals are at most two-fold, and all of them can easily be done by closing the contours to the right and summing up residues in terms of nested harmonic sums defined recursively by [93] S i (n) =  Integrating over the angles of the particles three and four yields two hypergeometric functions, and we introduce an MB representation for each of them. Performing the integration over the last angle, we obtain a four-fold MB representation for F 10 . Two integrations can immediately be perfumed using Barnes' lemmas, and we obtain After resolving the singularities in ǫ and expanding under the integration sign, all the twofold integrals can be reduced to onefold integrals using Barnes' lemmas and their corollaries. The remaining onefold integrals are trivial to compute by closing the contour and summing up residues. We find,

The master integral F 5
We start by replacing the sums in the denominator of the integrand of F 5 by three MB integrals, We insert the energies and angles parametrization for the final-state particles and introduce MB integrations for the angular integrals. Note that the first angular integral necessarily involves three massless propagators and can thus not be done in closed form as a 2 F 1 , so we insert the three-fold MB representation (5.17) for it. We then arrive at a representation for F 5 as a six-fold MB integral convoluted with two angular integrations. It turns out that after performing the change of variables z 3 → z 3 + z 5 and z 6 → z 6 + z 1 + z 2 , the integrals over z 1 , z 2 and z 5 can be done in closed form using Barnes' first lemma. The remaining two angular integrations can easily be performed in terms of hypergeometric functions, and all but one MB integration can be performed using Barnes' lemmas. We thus arrive at a one-fold MB representation for F 5 , One might be tempted to sum up the residues at z 1 = n and z 1 = −2ǫ + n, n ∈ N × , for finite values of ǫ to obtain an expression for F 5 as a combination of two hypergeometric functions at 1 valid to all orders in ǫ. The two hypergeometric functions are however separately divergent (even for finite values of ǫ) and only their sum is finite. We therefore only compute a Laurent series for F 5 . Resolving singularities in ǫ and summing up residues in terms of harmonic sums we obtain,  We start by introducing two MB integrations in order to remove the sums in the denominator of the integral F 8 . The energies are integrated out in terms of Γ functions, and the angular integrations over the particles four and five are performed using eq. (5.16). At this stage we have three integrations left to do: the two MB integrations and the integral over cos θ 3 = 2y − 1 In order to proceed, we first apply the identity and then insert an MB representation for each hypergeometric function in the integrand of eq. (8.30). The reason to apply eq. (8.31) before inserting the MB integrations comes from the fact that in this way one of the four MB integrations can be performed using Barnes' first lemma. We then arrive at the following three-fold MB representation for F 8 ,

(8.32)
In the rest of this section we show how we can compute a Laurent expansion for this integral. We proceed in the standard way and resolve singularities in ǫ. At the end of this procedure, we have a collection of MB integrals of dimensionality at most three with integration contours that are straight vertical lines. These integrals can then be safely expanded in ǫ under the integration sign. In the following we discuss the computation of the two and three-fold integrals.
Three-fold MB integrals. There is one three-fold integral contributing to F 8 , We start by closing the z 3 contour to the right and take residues, +i∞ −i∞ .

(8.35)
The first two sums can be performed in terms of Γ functions and their derivatives. We illustrate this on the first term (the second term is similar),

(8.36)
In this way, the first two terms can effectively be reduced to the computation of two-fold integrals, and we will therefore not discuss them any further in this section. The third term can also be summed up in closed form. However, unlike the first two terms, the sum cannot be expressed in terms of Γ functions and their derivatives alone, but the sum evaluates to a 4 F 3 function. We therefore arrive at the following single three-fold MB integral, × Γ (z 2 + z 4 + 1) 4 F 3 (1, 1, 1, z 2 + 1; 2, 2, 1 − z 4 ; 1) .

(8.37)
Although it looks like we have managed to reduce the three-fold integral to a two-fold integral, it is still secretly three-fold, except that we have 'hidden' one integration inside the Re(z 2 ) = 1 3 and Re(z 4 ) = 1 5 .

(8.38)
We then insert an Euler integral representation for the 4 F 3 function as an integral over a 3 F 2 . The 3 F 2 function turns out to be reducible, and so we finally arrive at Next, we would like to exchange the MB and the Euler integration and sum up the residues of the poles of the Γ functions in the integrand. However, we are only allowed to do so if the Euler integration does not produce any new poles whose residues need to be taken into account. It is easy to see that in our case the Euler integral converges whenever Re(z 2 ) and Re(z 4 ) are positive. We can thus close both contours to the right and exchange the Euler and MB integrations and then sum up the residues coming from the Γ functions. We start by taking residues in z 2 and then the residues in z 4 . There are several cases to be considered separately: 1. The poles at z 2 = n 2 ∈ N × give rise to poles in z 4 at n 4 ∈ N × . We hence obtain double sums of the form ∞ n 2 ,n 4 =1 Sums of this type can be performed using XSummer [94,95], and give rise to complicated multiple polylogarithms whose arguments are rational functions of t and (1−t).
Using symbols (see Appendix D), all these complicated multiple polylogarithms can be reduced to harmonic polylogarithms with indices 0 and 1 in t.

2.
Taking the residues at the poles at z 2 = −z 4 +n 2 , n 2 ∈ N × gives rise to the expression Next we want to close the z 4 contour and sum up the corresponding residues. From the previous discussion, we are forced to close the contour to the right, and the summand obviously only has poles at z 4 = n 4 ∈ N × inside the integration contour. There is however a subtlety, and we cannot just sum up the residues. Indeed, it is easy to see that • for n 4 < n 2 , there are double poles. Shifting the summation variable n 2 → n 2 +n 4 , these residues give rise to double sums similar to eq. (8.41) and can again be performed using XSummer. We obtain complicated multiple polylogarithms with rational functions of t as argument. Using symbols, they can again be simplified to harmonic polylogarithms with indices 0 and 1 in t.
• for n 4 = n 2 , there is a simple pole. This gives rise to a simple sum which is trivial to perform in terms of harmonic polylogarithms in t.
This sum is not of the type (8.41), and we therefore need a different way to sum up the series. This procedure will be discussed in the rest of this section.
One way to sum the series S(x, y) is to recognize that it is related to an integral over an Appell F 3 function. More precisely we can write S(x, y) = − 1 xy ∞ n 2 ,n 4 =0 n 2 ! n 4 ! (n 2 + n 4 + 1)!
x n 2 +1 n 2 + 1 y n 4 +1 n 4 + 1 = − 1 xy  The particular Appell F 3 function we obtain turns out to be reducible, and so we obtain a simple two-fold integral representation for S, This integral is trivial to perform in terms of multiple polylogarithms with rational functions in x and y as arguments (see Appendix D). It then follows that S(1 − t, 1 − 1/t) can be expressed in terms of multiple polylogarithms with rational functions in t as arguments.
Using symbols, we arrive at the following simple expression,

(8.47)
Putting everything together, we arrive at a representation for F 8,3 as a one-fold integral over harmonic polylogarithms. This integral is trivial to perform, and we obtain Two-fold MB integrals. There is only one two-fold MB integral contributing to F 8 that cannot be reduced to simpler integrals by Barnes' lemmas and their corollaries. This integral reads The integral could in principle be done by closing contours to the left and summing up residues. This leads to double sums of the form (8.41) -thousands of them due to the presence of multiple poles -but without any parametric dependence. We therefore took a different route, which we present in the following. We start by noting that the integrand of eq. (8.49) agrees, up to higher order terms of O(ǫ), with the function It would however be wrong to conclude that then necessarily the integrals, seen as a Laurent expansion in ǫ, are also equal to the same accuracy, because the Laurent expansion of eq. (8.51) around ǫ = 0 might require to shift the integration contours to avoid pinch singularities in the limit ǫ → 0. It is however easy to see that, for the contours given in eq. (8.50), no new pinch singularity is created for ǫ → 0 in eq. (8.51). We therefore consider from now on ǫ to be infinitesimal but finite: it is large enough to separate poles at, e.g., z 2 = 0 from z 2 = ǫ, but small enough to ensure that no poles change their nature, i.e., poles that were left (right) of the contour (8.50) in eq. (8.49) remain left (right) of the contour in eq. (8.51). If we chose ǫ in this way, we conclude that where F 8,2 is given by eq. (8.51) integrated over the straight vertical lines defined by eq. (8.50). Our aim will be to find an ǫ expansion for F 8,2 . We perform the change of variables (z 3 , z 4 ) → (−1 − z 3 , −z 4 ), and we close the z 3 contour to the right and take residues. There are two towers of poles we need to take into account: z 3 = n 3 ∈ N and z 3 = −1 − ǫ + n 3 , n 3 ∈ N × . (8.53) Two comments are in order: 1. At this stage our assumption that ǫ be infinitesimal but finite is vital, because otherwise the two towers of poles would glue together and thus give rise to double poles.
2. The second tower of poles runs over the set {−ǫ, −ǫ + 1, . . .}. For technical reasons that will become clear below, it is easier to explicitly take into account the residue at z 3 = −1 − ǫ, and to compensate for this by adding it back, where F 8,2 is the integral obtained from F 8,2 by deforming the z 3 contour such that the pole at z 3 = −1 − ǫ is now to the right of the contour. The residue at z 3 = −1 − ǫ can be computed in closed form and gives rise to (8.55) Next we compute F 8,2 by closing the z 3 contour to the right and summing up residues. The resulting sums can easily be performed in terms of hypergeometric functions, and we get In order to proceed, we insert an Euler integral representation for each of the 3 F 2 functions. It is then easy to see that the Euler integrals are convergent for Re(z 4 ) > 0 and ǫ infinitesimal but finite, and so we can exchange the Euler and MB integrations provided that we close the integration contour to the right. The important point is that the 2 F 1 functions appearing inside the Euler integrals are independent of z 4 , and so they can be pulled out of the MB integral. Summing up the residues in z 4 , we then arrive at the following integral representation for F 8,2 .
The terms involving a single 2 F 1 function in the integrand immediately evaluate to 3 F 2 functions, which can be expanded in ǫ using HypExp. In addition, the fourth term can be done in closed form as follows: We insert an MB representation for each 2 F 1 and perform the Euler integration as a Beta function. One of the two remaining MB integrals evaluates to a 2 F 1 evaluated at 1, which reduces to Γ functions through Gauss' identity. The remaining one-fold MB integral then immediately evaluates to a 4 F 3 function, which can be expanded in ǫ using HypExp. We were not able to find a closed form for the last remaining Euler integral. We therefore insert an Euler integration for each 2 F 1 , and we obtain the expression, (8.58) where I denotes the integral It is easy to see that I is finite as ǫ → 0, and so we can expand in ǫ under the integration sign and perform the integration over t, u and v recursively using the techniques of Appendix D. This is a trivial exercise that leads to We have thus obtained the ǫ expansion of F 8,2 , and thus of F 8,2 . We find, where γ E = Γ ′ (1) denotes the Euler-Mascheroni constant.
The result for F 8 . We have now computed all the two and three-fold integrals contributing to F 8 . The remaining one-fold integrals are trivial to compute and we obtain  We start in the usual way and derive an MB representation for F 9 by inserting the energies and angles parametrisation and inserting MB integrations for the angular integrals. After applying Barnes' lemmas and their corollaries several times, we arrive at the following MB representation for F 9 , F 9 (ǫ) = F 9,1 (ǫ) + F 9,2 (ǫ) , (8.64) where , We can resolve singularities for F 9,1 and F 9,2 and expand in ǫ up to and including O(ǫ). The result is a collection of high dimensional MB integrals which, after closing the contour and taking residues, result in multi-fold harmonic sums. While we were able to perform all the harmonic sums in terms of zeta values for all MB integrals up to O(ǫ 0 ), at O(ǫ) new polygamma functions appear in the integrand which make the combinatorics of the sums rather intricate. We therefore chose a different method to evaluate the integral F 9 , which we describe in the rest of this section. We start by noting that the F 9 is finite in D = 6 dimensions. This can easily be checked by replacing ǫ by ǫ − 1 in the MB representations (8.65) and (8.66) and resolving singularities. Our goal is to find a parametric integral representation for F 9 in D = 6 − 2ǫ dimensions and to expand under the integration and perform the parametric integrations recursively. The result in D = 6 − 2ǫ can then be related to the (divergent) result in D = 4 − 2ǫ using the dimensional recurrence relation for F 9 of Section 6.
It is easy to derive a parametric representation for F 9 using the technique described in Appendix C. We find with Several comments are in order about the parametric integrals we just defined. First, one can easily check that both I 9,1 and I 9,2 are individually finite as ǫ → 0. Second, at first glance our goal to integrate out the integration variables one-by-one seems rather hopeless due to the appearance of the huge polynomial factor. However, as we will see shortly, there is a sufficient condition that allows one to test whether a parametric integral can be performed in terms of multiple polylogarithms, and this criterion is fulfilled for the integrands of I 9,1 and I 9,2 . We very briefly summarize this criterion in the following, and we refer to ref. [96] or to Appendix D for more details. In order to understand the criterion, it is important to first understand multiple polylogarithms and their integration.
Suppose now that we are given some integral over some rational function. We would like to be able to perform the integrations one-by-one in terms of multiple polylogarithms. Unfortunately, such a procedure is not possible for a generic integral. In ref. [96] a sufficient condition was derived that allows one to test whether such a recursive integration is possible. For a detailed description of the criterion we refer to ref. [96], or to Appendix D, where (a version of) the criterion is reviewed. In a nutshell, as can easily be seen from the definition of multiple polylogarithms, a sufficient condition to perform the parametric integrations recursively using multiple polylogarithms is that at each integration step we can find an integration variable in which all the denominators appearing in the integrand are linear. More formally, let S (i) be the set of all polynomial factor in the integrand of I 9,i . Next, we define the sets S (i) (x,y,...) , where x, y, . . . are integration variables, as the sets of irreducible non-monomial polynomial factors in the integrand after integration over (x, y, . . .) in this order. For more details how to construct these sets we refer to ref. [96] and to Appendix D. Here it suffices to say that the sets S (i) (x,y,...) can be constructed solely from the knowledge of the sets S (i) . A sufficient criterion for a recursive integration in terms of multiple polylogarithms is that there is an ordering of the integration variables such that the sets S (i) (x,y,...) corresponding to this ordering contain at least one integration variable in which all the polynomials of the set are linear. Miraculously, this condition is satisfied for both I 9,1 and I 9,2 . We have We see that if we perform the integration in the order (t 1 , x 1 , x 2 , x 3 , t 2 ) then at each step all the polynomials are linear in the next integration variable. The actual integration can be carried out in an algorithmic way. The procedure is however rather lengthy, so we do not discuss it here in detail, but we refer to Appendix D for a detailed description of the integration algorithm. The result is

Conclusions and outlook
In this article, we presented a method for the expansion of phase-space integrals in kinematic parameters. We applied our method to perform a threshold expansion of triple-real radiation partonic cross-sections for Higgs production at hadron colliders. Our method reduces the phase-space integrals in the coefficients of the expansion to a small number of soft master integrals. These master integrals emerge not only in Higgs production processes but in all processes which are topologically equivalent, such as Drell-Yan production. We have computed the master integrals using a combination of techniques for Mellin-Barnes integration and the integration of real parametric integrals. We presented explicit expressions for the first two coefficients in the threshold expansion of all triple-real radiation partonic cross-sections which are necessary for the inclusive Higgs cross-section at N 3 LO, both in terms of soft master integrals and as a Laurent series expansion in the dimensional regulator ǫ.
We can apply our method to compute further coefficients in the threshold expansion. This requires the reduction of higher rank cut loop integrals in their powers of propagators and irreducible numerators. It also requires the evaluation of some new master integrals originating from diagrams which did not contribute to the first two terms of the threshold expansion. We believe that it is feasible to compute in this way a sufficient number of terms in the threshold expansion for phenomenology purposes. A second possibility is to apply the reverse-unitarity method to compute the same cross-sections without resorting to a threshold expansion. The work that we have presented here will be particularly important for solving the unexpanded master integrals. With reverse unitarity, one can derive differential equations for the master integrals in arbitrary kinematics. For their solution, we can use the soft limits presented here as boundary conditions. Subleading terms in the expansion can serve as a check of the solutions or as an aid to guess them.
A complete computation of the N 3 LO coefficient requires also virtual contributions integrated over the phase space of up to two real partons in the final state. The corresponding one-loop, two-loop and three-loop amplitudes are already known in the literature as a Laurent expansion in ǫ [99][100][101][102][103][104][105][106][107][108][109]. Integrating them over phase-space requires in addition their universal single-real and double-real infrared limits at higher orders in ǫ. A threshold expansion of mixed real and virtual corrections requires an extension of our method. While the threshold expansion for real emission matrix elements can be performed by uniformly rescaling all final-state parton momenta, p i →zp i , and expanding in the small parameter z, for real-virtual matrix-elements the scaling of the loop momentum cannot be determined without further consideration. To determine the correct expansion at threshold we consider the strategy of regions [110,111]. We applied this method to all 2-loop real-virtual master integrals contributing to Higgs production at NNLO and could reproduce the soft limits of the integrals. In the following we will illustrate this technique on the example of the integral We introduce the light-cone coordinates, e.g. for q µ , and q µ ⊥ denotes the remaining (D − 2) Euclidian transverse components. Furthermore we choose our reference frame such that We make the soft region manifest by redefining the loop momentum and the final state parton momentum in terms of the soft-expansion parameterz After expanding inz we find the leading soft integral to be .
All external kinematic scales are factored out of the integrand and the integral can be easily computed using Schwinger parameters. This yields in agreement with the result in ref. [67]. In full kinematics other regions, which are suppressed in the soft limit, contribute to the integral. Finally, the computation of the collinear counterterms for the partonic cross-section has been performed in refs. [112,113]. Alternatively, this task can easily be achieved numerically following the procedure of ref. [36] and using the results for the partonic cross-sections at NNLO through O(ǫ) of refs. [67,68].
We believe that the prospects for a complete computation of the N 3 LO coefficient for Higgs production and other processes are excellent, and we are looking forward to completing this task in future publications.
After eliminating the integration over the Higgs momentum using momentum conservation the measure becomes It is well known that rotationally invariant integrals of the form can be rewritten by making the integration over the scalar products manifest: (A.11) Here Gram(p 1 , . . . , p N ) denotes the Gram determinant of the vectors p i , which is defined as We observe that the measure as well as the matrix element are invariant under rotations in the transverse space. Therefore, we can use this relation to rewrite the D − 2 dimensional transverse integration Using the definition (A.7) of the s-channel invariants, we can rewrite the integration over the scalar products as The integration over the norm of the transverse components can be performed using the on-shell δ-functions, so that we find where we have written the Gram determinant as Thus we finally arrive at the phase space measure as in eq. (4.4) and

B. Derivation of the phase space volume
In this appendix we derive an expression for the 2 → H + (N − 2) phase space volume that is valid for all orders in z and ǫ. We start by recalling the general phase space factorisation: If we assume that we know the phase space volume at N k LO, i.e. for 2 → H + k, this relation allows us to rewrite the phase space volume at N k+l LO as We are specifically interested in rewriting the phase space volume at N k+1 LO as a convolution of the N k LO phase space with a two particle phase space. Specialising to l = 1, we can use this formula to inductively derive the phase space volume for arbitrary orders. We have, From eq. (A.17) we obtain an explicit parametrisation of the two particle phase space and we find We can perform the integral over x 1,k+3 and x 2,k+3 in terms of beta functions and make the transformation µ 2 = sx, m 2 = sz to obtain In the following we use eq. (B.6) to prove inductively the following result: First, eq. (B.7) correctly describes the phase space volume for n = 1 and n = 2. In order to derive Φ n+2 iteratively from Φ n+1 we use eq. (B.7) and find To solve the integral, we make the transformation x = 1 − (1 − z)y and find where we have factored out Next, we introduce a Mellin-Barnes representation for 2 F 1 and exchange the MB integration with the parametric integration. This allows us to perform the integration over y in terms of another 2 F 1 . Then we find Introduce another MB integral for the 2 F 1 , we arrive at (B.12) In the next step we perform the change of variables ξ → ξ + χ, such that (B.13) The integral over χ can now be performed using Barnes' first lemma, and we find (B.14) The integral over ξ is just the MB representation of a 2 F 1 so that we finally find We have therefore inductively shown the validity of (B.7) for all n.

C. From Mellin-Barnes integrals to parametric integrals
In this appendix we describe how one can derive an Euler-type integral from a Mellin-Barnes integral with a balanced integrand. Roughly speaking, a MB integral is said to be balanced if for each integration variable z i the number of Γ functions of the form Γ(. . . − z i ) is equal to the number of Γ functions of the form Γ(. . . + z i ). More precisely, the integral We assume in the following that the contours are straight vertical lines such that the real parts of the arguments of all the Γ functions are positive 6 . In that case we can always derive an Euler-type integral representation for the MB integral. We start by noting that if an integral is balanced, then we can always express its integrand as a product of Beta functions, The integral (C.2) is convergent whenever Re(x), Re(y) > 0. It is easy to convince oneself that this condition is satisfied whenever the real parts of all arguments of the Γ functions are positive in the original MB integral. We can therefore replace each Beta function by its integral representation (C.2) in the integrand of the MB integral and, because all the integral are convergent, we can exchange the MB integrations and the integrations coming from the Beta functions. This leaves us with an integral of the form where t = (t 1 , . . . , t N ) and the R k are ratios of products of the t i and 1 + t i . Next, we would like to perform the MB integrations. This can be done using the formula Indeed, parametrizing the contour as z = z 0 + it, we obtain We can solve the δ constraints, and the result is the desired parametric integral. In the following we discuss two very simple examples that illustrate the above procedure. We start by discussing Barnes' first lemma, i.e., we consider the integral We assume that the integration contour and a, b, c and d are such that the real parts of all Γ functions are positive. We rewrite the integrand in terms of Beta functions, where the last step follows from eq. (C.4). Solving the δ-function constraint leads to a one-fold integral that can immediately be recognized as a Beta function, and we recover the usual form of Barnes' first lemma, (C.9) The second example we are going to discuss is Gauss' hypergeometric function. We consider the integral We again assume that all conditions for convergence are satisfied. Rewriting the integrand in terms of Beta functions, we obtain (C.11) Solving the δ-function constraint with respect to t 2 and performing the change of variables t 1 → ξ/(1 − ξ), we immediately arrive at the usual integral representation for the 2 F 1 function, (C.12)

D. Symbolic integration of parametric integrals
In this section we describe an algorithmic approach to compute certain classes of parametric integrals. We stress that this approach is not genuinely new, but has already been successfully applied, in some variant or another, to the computation of Feynman integrals, e.g., [96,[114][115][116][117][118]. More precisely, consider an integral of the form where a k and b k are integers and the P k ({x i }; {y j }) are polynomials with integer coefficients, which we assume irreducible over Z, i.e., they cannot be factorized into a product of nonconstant polynomials of lower degree. We assume that the integration range is [0, ∞] for each integration variable x i . While many statements remain true for generic integration boundaries, in certain cases the algorithm we are going to describe breaks down, due to the appearance of boundary contributions. We can however always map a generic integration region x i ∈ [a, b] to y i ∈ [0, ∞] by the change of variable We furthermore assume that the integral is convergent for ǫ = 0, and so we can expand in ǫ under the integration sign. Note that in the simplest divergent cases where the singularities in the integrand factorize we can reduce the problem to a convergent integral by subtracting the divergencies.
Our goal is to compute the coefficients of the Taylor expansion of I({y j }; ǫ) by integrating out the integration variables recursively one-by-one in terms of multiple polylogarithms, eq. (8.71). Obviously, not every integral can be performed in this way. In the rest of this appendix we describe a sufficient condition, first obtained in ref. [96], to be satisfied by the integrand of I({y j }; ǫ) so that it can be integrated in terms of multiple polylogarithms. We then review an algorithm for integrating these classes of integrals, and illustrate the procedure explicitly on a simple example.

D.1 Denominator reduction
In this section we review (a variant of) the sufficient condition of ref. [96] to determine whether an integral can be performed recursively in terms of multiple polylogarithms. The main idea is that we need to determine an ordering of the integration variables such that at each step during the integration all the denominators are linear in the next integration variable. We start by defining the set S of all the polynomials that are not monomials and that appear inside the integrand of eq. (D.1), To start the integration, we have to assume that there is one integration variable, say x a such that all the element of S are linear in x a . In that case we may write where Q a k ({x i }; {y j }) and R a k ({x i }; {y j }) are polynomials that are independent of x a . Note that, after expansion in ǫ, the integrand may also contain logarithms of the P k ({x i }; {y j }), which can be rewritten in terms of multiple polylogarithms, where for clarity we suppressed the arguments of the polynomials. Furthermore, we can use the shuffle algebra of multiple polylogarithms to replace every product of multiple polylogarithms by a sum, G(a 1 , . . . , a n 1 ; z) G(a n 1 +1 , . . . , a n 1 +n 2 ; z) = σ∈Σ(n 1 ,n 2 ) G(a σ(1) , . . . , a σ(n 1 +n 2 ) ; z), (D. 6) where Σ(n 1 , n 2 ) denotes the set of all shuffles of n 1 + n 2 elements, i.e., the subset of the symmetric group S n 1 +n 2 defined by Σ(n 1 , n 2 ) = {σ ∈ S n 1 +n 2 | σ −1 (1) < . . . < σ −1 (n 1 ) and σ −1 (n 1 + 1) < . . . < σ −1 (n 1 + n 2 )} . (D.7) Thus we can assume without loss generality that the integration over x a takes the form (D.8) Partial fractioning the factors in the denominator, e.g., we obtain a sum of integrals that can be reduced to the recursive definition of multiple polylogarithms, eq. (8.71). We can thus easily compute a primitive with respect to x a , and then take the limits x a → 0 and x a → ∞. We will address the issue of limits in the next subsection.
We would like to iterate this procedure and integrate over the next variable. This is however only possible if inside the new integrand we can still find an integration variable in which all polynomials are linear. The polynomials appearing inside the integrand are Q a k and R a k , which have been introduced through eq. (D.5) and (D.9), as well as the combinations Q a k R a l − Q a l R a k introduced by partial fractioning. This last polynomial however need not necessarily be linear, even if the Q a k and R a k are. In order to proceed, it is therefore mandatory that all the Q a k R a l − Q a l R a k factor into polynomials that are linear in a certain variable.
In ref. [96] a criterion was given that allows one to determine a priori whether the above procedure terminates, i.e., whether there is an ordering of the integration variables such that all the denominators stay linear at each integration step. We start by defining the set S (xa) as the set of irreducible factors that appear inside the polynomials Q a k , R a k and Q a k R a l − Q a l R a k . Then, if we can find an integration variable x b such that all the elements of S (xa) are linear in x b , we can restart the above procedure and integrate over x b . If we iterate this procedure and are able to construct a sequence of sets of polynomials such that in each set all the polynomials are linear in at least one integration variable, then we have found an ordering of the integration variables such that we can recursively integrate out all the integration variables in terms of multiple polylogarithms. We stress that this condition is sufficient, but not necessary: even if we fail to find a suitable sequence (D.10), the integral might still be expressible in terms of multiple polylogarithms (e.g., after a suitable change of variables). In addition, note that for the last integration step it is not necessary for the polynomials to be linear: we can then factor the polynomials into linear factors, whose roots involve algebraic expressions of the parameters {y j }.

D.2 Symbolic integration
We have now a criterion to determine if a given integral of the type (D.1) can be integrated recursively in terms of multiple polylogarithms, and we have already explained how to perform the first integration step. We assume from now on that we have found a sequence (D.10) and the associated ordering of the integration variables, and that we have performed the integration over x a as described in the previous subsection.
In order to continue and perform the next integration, we need to address two issues 1. We need to take the limits x a → 0 and x a → ∞ of the primitive with respect to x a .
2. If we want to integrate next over x b , we have to face the problem that x b might be 'hidden' inside the arguments of multiple polylogarithms of the form i.e., x b appears inside the polynomials Q a k and R a k . In order to compute a primitive with respect to x b , we need to rewrite all multiple polylogarithms of the form (D.11) as G( a; x b ), where a is independent of x b .
The limit x a → 0 can easily by taken by using the fact that While eq. (D.12) is sufficient to compute the value at x a = 0 in many circumstances, it can happen that the primitive has spurious poles and / or logarithmic singularities at x a = 0.
In the presence of poles we need to expand the multiple polylogarithms around x a = 0. This task can easily be achieved using the series representation of multiple polylogarithms. We can thus easily obtain the expansion around x a = 0 in terms of Z-sums to arbitrary order. Next, we discuss how to take the limit x a → ∞. If we definex a = 1/x a , then the problem can be reduced to taking the limitx a → 0, We would like to use the inversion relations for the multiple polylogarithms in the righthand side to express them as linear combinations of multiple polylogarithms of the form G (. . . ;x a ), for which the limitx a → 0 is trivial. Unfortunately, the corresponding functional equations are often unknown. Note at this point that the problem of finding the inversion relations is formally equivalent to the second problem we have not yet addressed, namely how to bring the primitive with respect to x a into a form where x b only enters through multiple polylogarithms of the form G( a; x b ): in both cases we are looking for functional equations that bring a certain multiple polylogarithm into a 'canonical form' where a certain variable only appears as the explicit upper integration boundary of the multiple polylogarithms. In the rest of this section we argue that, under certain conditions which we will later show to be equivalent to the criterion derived in the previous subsection, plus the hypothesis that the integration ranges are [0, ∞], we can always bring multiple polylogarithms into the 'canonical form' for some variable x such that a i is independent of x and the coefficients c i involve only multiple polylogarithms that are independent of x. Note that this result is similar to the result obtained in ref. [96]. In the following we give a constructive algorithm that allows us to derive the canonical form (D.17).
In order to achieve a rewriting of our multiple polylogarithms in canonical form, we need to derive the corresponding functional equations. The natural language to discuss functional equations among multiple polylogarithms are symbols [119][120][121][122][123] and the Hopf algebra of multiple polylogarithms [98]. We start by giving a concise review of symbols.
(D. 19) As an example, the symbols of the classical polylogarithms and the ordinary logarithms are given by where ∐ ∐ denotes the shuffle product on tensors. We can make the following observation about the symbol of a multiple polylogarithm: if the a i are independent of x, then the symbol of G(a 1 , . . . , a n ; x) contains exactly one term hat contains x in all its entries, and this term can be chosen to be of the form S(G(a 1 , . . . , a n ; x)) = (a n − x) ⊗ . . . ⊗ (a 1 − x) + . . . .
= G(a 2 , . . . , a n ; x) d log(a 1 − x) + . . . , where the last step follows from and where the dots indicate terms in the total differential that are independent of dx. The statement then follows recursively. Note that this statement is independent of whether the a i are zero or not. Assume now that we are given an integrable symbol T (which will correspond later to the symbol of the multiple polylogarithm we want to bring into canonical form) which is of uniform weight w and has rational coefficients. If T does not satisfy this last condition, we deal separately with the contributions of different weight and / or different rational function prefactor. Let us suppose that the entries are drawn from a set S. Without loss of generality, we may assume S to consist of irreducible polynomials in some variables x i , 1 ≤ i ≤ n, and for simplicity we assume for now that the polynomials are linear in all the x i (we will see in the next subsection that the correct restriction is that S satisfies the reduction criterion of the previous subsection). Furthermore, we assume that we have fixed an ordering on the variables, which we will take in the following as (x 1 , . . . , x n ). For each variable x i , we define a linear map φ x i which acts on elementary tensors s by Morally speaking, the map φ x i assigns to T the combination of multiple polylogarithms which will give the same terms that have x i in all entries through eq. (D.22). Using the maps φ x i , we can now formulate an algorithm that assigns to T a multiple polylogarithm in the canonical form associated to the ordering of the variables (x 1 , . . . , x n ). We start by defining a new symbol by subtracting off the contribution from φ x 1 (T ), (D. 26) By construction, each term in the symbol T 1 has at least one entry that is independent of x 1 . Next, concentrate on the terms of the form where by hypothesis the b i k are independent of x 1 . As we have subtracted of the contribution from φ x 1 (T ), these terms cannot come from a multiple polylogarithm of the form G(. . . ; x 1 ) of weight w, but it can only arise from the product It is easy to convince oneself that the difference contains only terms for which at most (w − 2) entries depend on x 1 . We can now go on an recursively subtract contributions with different multiplicities of x 1 . Assume for example that we have subtracted all contributions where x 1 appears in more than (w − r) entries, and that the resulting symbol is T r . We can then concentrate on the terms (i 1 ,...,iw) It is easy to convince oneself that in the difference x 1 appear in at most (w − r − 1) entries. We continue this procedure until we reach T w , which is independent of x 1 , and we restart the algorithm with T w and φ x 2 . We then repeat this procedure until we have exhausted all the integration variables, and the algorithm stops. The result of this algorithm is, by construction, a function of the form (i 1 ,...,in) c i 1 ,...,in G( a in ; x n ) . . . G( a i 1 ; x 1 ) , (D.32) whose symbol is T and such that the a i k are sequences of rational functions that are in the variables x k , k > i k , i.e., the sought-for canonical form for T . Using this algorithm we can easily take the limit x a → ∞, by applying it to the righthand side of eq. (D.16) with the ordering (x a , x b , . . .). We can also use it to rewrite the integral over x a such that we can easily integrate over the next integration variable x b . However, the result of the algorithm will at this stage still deliver the wrong answer, as we have so far only constructed a function in canonical form whose symbol matches the symbol of a given function. To illustrate this, assume we run the algorithm on the function in the right-hand side of eq. (D.16) with the ordering (x a , x b , . . .). The result is a function G(x a , x b , . . .) in canonical form such that It would however be wrong to conclude that G a; 1 xa and G(x a , x b , . . .) are equal, because the symbol maps to 0 all terms proportional to multiple zeta values. In ref. [124] an algorithm was described that allows to reconstruct these zeta-valued terms using the full Hopf algebra structure of multiple polylogarithms, augmented by some ideas by Brown [125]. In the following we only give a very brief account on how to reconstruct the zeta-valued terms, and we refer to ref. [124] for a detailed description of the algorithm.
In the following we assume for simplicity that the function G a; 1 xa − G(x a , x b , . . .) is real 7 , so we do not need to worry about imaginary parts proportional to iπ. Next we act with ∆ 2,1,...,1 on the difference, where ∆ 2,1,...,1 is the component of the iterated coproduct where the first component has weight two and all other components have weight one. Without loss of generality we can write where the a i k are irreducible polynomials and A i 1 is a combination of multiple polylogarithms of weight two. Without loss of generality we can assume that we have collected all term that have the same 'tail' log a i 2 ⊗ . . . ⊗ log a i w−1 , i.e., we can assume log a i 2 ⊗ . . . ⊗ log a i w−1 = log a j 2 ⊗ . . . ⊗ log a j w−1 if (i 2 , . . . , i w−1 ) = (j 2 , . . . , j w−1 ) . (D. 35) As we know that the symbol of the function vanishes, eq. (D.33), we necessarily conclude that A i 1 is proportional to ζ 2 , i.e., A i 1 = k i 1 ζ 2 , for some rational number k i 1 . This rational number can easily be determined by evaluating A i 1 numerically at a single point using any of the standard libraries to evaluate multiple polylogarithms [84,[126][127][128][129][130][131], and running for example the PSLQ algorithm [132]. Equation (D.34) then takes the form Next we drop ζ 2 , i.e., we only keep the tail of each elementary tensor. If we also drop the log signs, we obtain a symbol associated with the terms proportional to ζ 2 , Running the algorithm described at the beginning of this section on the symbol in eq. (D.37) we obtain a function G 2 (x a , x b , . . .) of weight w − 2 in canonical form such that We have in this way determined all the contributions proportional to ζ 2 , and the result is by construction in canonical form. We then repeat exactly the same exercise by acting with ∆ 3,1,...,1 to determine the terms proportional to ζ 3 , and we continue in this way until we have exhausted all possibilities and the algorithms stops. As a result we obtain the expression of G a; 1 xa in canonical form at function level. This terminates the algorithm to bring multiple polylogarithms into the canonical form corresponding to a certain ordering of the variables. This operation is sufficient to take the limit x a → ∞, and to bring the integral into a form where the next integration can easily be performed. We have however not yet shown that this algorithm will always terminate on the class of integrals we consider. In particular, when formulating the algorithm, we assumed that all entries that appear in the symbol are linear in all variables. In the next subsection we show that this condition can be relaxed, and that the algorithm always terminates if the integrand satisfies the reduction criterion of the previous subsection and if the integration range is [0, ∞].

D.3 Symbolic integration and denominator reduction
In this subsection we show that the algorithm we just described always terminates for integrals of the type (D.1) such that the set S of polynomials satisfies the reduction criterion. From now on we assume that we have found an ordering of the integration variables, which we take as (x 1 , . . . , x n ) and we can find a sequence S, S (x 1 ) , S (x 1 ,x 2 ) , S (x 1 ,x 2 ,x 3 ) , . . .
(D. 39) such that all the elements of S (x 1 ,...x k ) are linear in x k+1 . We start by showing that under these hypotheses and after having integrated out (x 1 , . . . , x k−1 ), 1. the symbol of the primitive with respect to x k has all its entries drawn from the set S (x 1 ,...x k−1 ) = {x k , . . . , x n } ∪ S (x 1 ,...x k−1 ) .
We start by proving the first statement, the second then immediately follows by taking the appropriate limits. As in the following we constantly switch between polynomials of the sets we just defined, and polynomials as entries of symbols, we introduce the notation [P ] to refer to P as an entry of a symbol. Note that we then have the identities We proceed by iteration in the number k of variables we have already integrated out. We start by analyzing what happens after the first integration. It is obvious from eqs. (D.5) and (D.9) that after the first integration the primitive only involves multiple polylogarithms G(a 1 , . . . , a w ; It is then easy to see (e.g., from the polygon approach to the symbol of ref. [123]) that the symbol can only have the following entries: (D. 42) The polynomials that appear inside the symbol are precisely those that appear in S and S (x 1 ) , which finishes the proof of the first statement for the first integration. Next consider taking the limits x 1 → 0 and x 1 → ∞. By definition, R 1 k and Q 1 k are independent of the limit, so we only need to consider the limits of [x 1 ] and [P k ]. [x 1 ] will give rise to logarithmic singularities in the limit, and these terms must cancel if the integral is convergent. For [P k ] we have, The logarithmic singularity in the last line must again cancel for convergent integrals, and so we see that the only polynomials that appear in the limit are those in S (x 1 ) . This finishes the proof of the second statement for the first integration. We stress that it is important that the integration region is [0, ∞], because otherwise we have to take into account effects coming from the integration boundaries.
Let us now suppose that the two statements are true for the first r − 1 integrations. The set S (x 1 ,...,x r−1 ) then consists of polynomials of the formP k =Q r k x r +R r k . Let us now compute the primitive with respect to x r . We start by running the algorithm to bring the multiple polylogarithms in the integrand into canonical form. As this involves the application of the map φ xr , we obtain multiple polylogarithms of the form G(a 1 , . . . , a w ; x r ) with a i ∈ 0, −R r k /Q r k . (D.44) If we compute the primitive, we integrate over linear functions in x r from the set S (x 1 ,...,x r−1 ) , and it is easy to see that this does not change eq. (D.44). Using exactly the same argument as for the first integration, we see that these multiple polylogarithms only contribute terms to the symbols whose entries are drawn from S (x 1 ,...x r−1 ) . In addition we have multiple polylogarithms that are independent of x r , whose symbols involve those elements of S (x 1 ,...x r−1 ) that are independent of x r . As these functions do not change if we take the primitive with respect to x r , we see that they do not alter the conclusion. So the symbol of the primitive with respect to x r must have all its entries drawn from S (x 1 ,...x r−1 ) . Taking the limits x r → 0 and x r → ∞ just like for the first integration then finishes the proof.
Having proved the two statements, we can show that our algorithm always terminates for the class of integral we consider. More precisely, we have to show that our algorithm can always produce the canonical form for the next integration step. This is done by applying the map φ xr , which requires all the entries in the symbol to be either independent of x r or linear in x r . By construction, this condition is always fulfilled for the elements of S (x 1 ,...x r−1 ) . It is easy to check that the same argument shows that the map φ x r+1 , which is called recursively by φ xr is well-defined, and so we can always find a canonical form for the integrand.

D.4 A toy example
In this section we discuss an example that illustrates how the algorithm described in the previous subsection can be used to compute parametric integrals. In particular, we used this algorithm to compute the parametric integrals (8.68) and (8.69). While these integrals are a straightforward application of the algorithm, intermediate expression are rather long, so we prefer to use a simpler integral where we can explicitly show all the steps. The integral we are going to consider is (D. 45) The integral is finite as ǫ → 0, and we want to compute the first few terms in the Taylor expansion I(ǫ) = I 0 + I 1 ǫ + I 2 ǫ 2 + I 3 ǫ 3 + O(ǫ 4 ) . (D.46) The first coefficient is trivial to compute, (D.47) We will now illustrate our algorithm in detail on the coefficient I 1 . If we integrate out the integration variables in the order (x 1 , x 2 , x 3 ), we obtain S = {1 + x 1 , 1 + x 2 , 1 + x 3 , 1 + x 2 + x 3 + x 1 x 3 } , (D. 48) We see that all the sets are linear in all integration variables, so we perform the integrations one after the other. The coefficient I 1 is given by the integral where we have already written all logarithms in terms of multiple polylogarithms, e.g.,

(D.50)
It is easy to compute a primitive with respect to x 1 for the integrand of I 1 , e.g., In the following we only concentrate on this single term (which is in fact the most complicated one) to illustrate the procedure. All other terms can be dealt with in a similar way. Before we compute the limit of eq. (D.51) as x 1 → 0, ∞, let us comment about the symbol of the primitive. We have (D. 52) We see that the entries in the symbol are drawn from the set S ∪ {x 1 , x 2 , x 3 }, as expected.
The same is true for all other terms in the primitive.
The limit x 1 → 0 is trivial, (D.53) The limit x 1 → ∞ is obtained by letting x 1 = 1/x 1 and deriving the inversion relation for this multiple polylogarithm. This is equivalent to bringing G (−1, (−x 2 − x 3 − 1)/x 3 ; 1/x 1 ) into canonical form with respect to the ordering of variables (x 1 , x 2 , x 3 ). We start by using our algorithm to construct a function G(x 1 , x 2 , x 3 ) which is in canonical form and has the symbol given in eq. (D.52). We find It is easy to check that the symbol of G(x 1 , x 2 , x 3 ) agrees with eq. (D.52). It would however be wrong to conclude that the functions are equal -they might differ by a rational multiple of ζ 2 . Evaluating the difference numerically at a single point, we obtain Taking the limit x 1 → ∞ is now trivial, and we obtain

(D.56)
Note that the function has a logarithmic singularity for x 1 → ∞, which will cancel against similar contributions from other terms. Furthermore, note that the symbols of all the (finite) terms in eq. (D.56) have entries drawn form the set S (x 1 ) ∪ {x 2 , x 3 }, as expected.
The previous steps can easily be implemented into a computer code and repeated for all the terms appearing in the primitive with respect to x 1 . The result of the integration over x 1 is, by construction, already in canonical form with respect to (x 2 , x 3 ), and we can immediately compute the primitive with respect to x 2 , e.g., (D.57) The limit x 2 → 0 is again trivial, while the limit x 2 → ∞ can again be computed by letting x 2 = 1/x 2 and constructing a function G(x 2 , x 3 ) in canonical form with the same symbol. We find , 0;x 2 − G(0, −1;x 2 ) + G(0, 0;x 2 ) − G (0, −1; x 3 ) .

(D.58)
Numerical evaluation at a single point immediately shows that we do not need to add any term proportional to ζ 2 . Taking the limit is now trivial, and we get Next we have to determine terms proportional to ζ 2 . This is achieved by acting with ∆ 2,1 on the difference, Taking the limitx 3 is now trivial, and we finally get I 1 = −5ζ 3 + 2π 2 9 + 5 3 .
(D. 65) The higher terms in the ǫ expansion can be obtained in exactly the same way. For this particular integral we find for example (D. 66)