Soft Expansion of Double-Real-Virtual Corrections to Higgs Production at N$^3$LO

We present methods to compute higher orders in the threshold expansion for the one-loop production of a Higgs boson in association with two partons at hadron colliders. This process contributes to the N$^3$LO Higgs production cross section beyond the soft-virtual approximation. We use reverse unitarity to expand the phase-space integrals in the small kinematic parameters and to reduce the coefficients of the expansion to a small set of master integrals. We describe two methods for the calculation of the master integrals. The first was introduced for the calculation of the soft triple-real radiation relevant to N$^3$LO Higgs production. The second uses a particular factorization of the three body phase-space measure and the knowledge of the scaling properties of the integral itself. Our result is presented as a Laurent expansion in the dimensional regulator, although some of the master integrals are computed to all orders in this parameter.


Introduction
The discovery of the Higgs boson [1] opens the path to the exploration of the sector responsible for the breaking of the electroweak symmetry and makes the Standard Model a fully predictive theory, because all couplings are now uniquely fixed. Precision measurements of the properties of the Higgs boson therefore represent a crucial test of the Standard Model, and any deviation, however small, will open a window to new physics scenarios.
The main production mechanism of a Higgs boson at hadron colliders is gluon fusion. Unfortunately, gluon-fusion suffers from large perturbative instabilities [2,3,4,5], and the uncertainty on the next-to-next-to leading order (NNLO) prediction still remains of the order of 10%. In order to fully exploit the potential of the LHC, it is therefore crucial to improve on the NNLO result by going to the next order in the perturbative expansion.
Recently, the next-to-next-to-next to leading order (N 3 LO) QCD correction to Higgs boson production has been computed as an expansion around threshold, and it was shown that the remaining QCD scale uncertainty is reduced to a small 2% [6]. Most of the ingredients and subprocesses that go into the computation of ref. [6] have been published separately over the last few years: the three-loop corrections to Higgs production in gluon fusion, as well as the corrections from the emission of an additional parton at one or two loops, have been computed in full generality in refs. [7,8,9,10,11]. In order to obtain a finite result, appropriate ultra-violet and infrared counterterms need to be included [12,13,14]. The contributions from the emission of two partons at one loop and three partons at tree-level, however, have not been computed for arbitrary kinematics so far, but they are only known as an expansion around threshold. In particular, in refs. [15,16,17,18,19] the first two terms in the threshold expansion of the triple-real and double-realvirtual corrections have been computed by reducing the corresponding phase-space and loop integrals to a small set of master integrals, all of which belong to a special class of integrals that we call soft integrals. These soft integrals are not only important when computing the first two terms in the threshold expansion of the cross section, but they also contribute to the full result for the cross section as boundary conditions to the differential equations satisfied by the master integrals in full kinematics. Despite their importance, the results for the individual soft integrals that contribute to the soft-virtual and next-to-soft corrections at N 3 LO of refs. [16,19] have never been published explicitly 1 .
The purpose of this paper is threefold: first, we close the aforementioned gap and we present all the technical details that went into the computation of the N 3 LO cross section in the soft-virtual and next-to-soft approximations of refs. [16,19]. We discuss in detail methods to perform the threshold expansion of one-loop matrix elements for the production of a heavy colorless state in association with two partons to any desired order, at least in principle. Second, we present techniques to compute the resulting soft master integrals analytically as a Laurent expansion in the dimensional regulator with coefficients that are polynomials in multiple zeta values. Finally, we give explicitly all the soft master integrals contributing to the results of refs. [16,19], as well as several new soft master integrals which only contribute at higher orders in the threshold expansion [6].
In a nutshell, our approach to threshold expansion can be described as follows: we start by expanding both the phase-space measure and the interference amplitudes into a power series in some small parameter quantifying the deviation from threshold. At every order in the expansion, the phase-space integrals are mapped onto (cut) Feynman integrals by virtue of reverse unitarity [3,20], which identifies on-shell phase-space constraints with cuts of Feynman propagators. The resulting cut integrals are then reduced to a small set of soft master integrals using integration-by-parts (IBP) techniques [21]. Similarly, the expansion of the one-loop matrix elements in the integrand is performed using the strategy of regions [22], which allows to exchange the threshold expansion and the loop integration, provided that contributions from all 'regions' in loop momentum space that lead to singularities in the soft limit are included. Similar to the findings of ref. [8], we observe that in the present case the relevant regions correspond to regions where the loop momentum can be either 'hard', 'soft' or 'collinear' to one of the initial-state momenta. While the IBP reduction of the hard and soft regions can be dealt with using standard techniques, the collinear regions require some special attention, because the resulting integrals contain propagators which are non-linear in the Lorentz invariants, and they are hence not amenable to standard techniques. We find that, due to the analytic structure of these propagators, it is always possible to recast the collinear-type integrals into a more standard form using partial fractioning, and we observe that the ensuing master integrals are always soft integrals, independently of the region.
In a second part of the paper we discuss techniques to evaluate the soft master integrals as Laurent expansions in the dimensional regulator, and we present two complementary approaches to achieve this. The first technique is reminiscent of the method of ref. [15] and allows one to obtain a multifold Mellin-Barnes representation with poles at integer values for each soft master integral. The second method uses a particular factorisation of the three-body phase-space measure, similar to the method presented in ref. [23] for purely real emissions. This method allows to arrive at a simple parametric integral representation, whose simplicity can be understood via the scale invariance properties of soft integrals.
The paper is organised as follows. In Section 2 we describe the setup of the problem and introduce some notation. In Section 3 we outline the methods used for the expansion around threshold and the technique of reverse unitarity. We discuss the reduction to master integrals and list the complete set of master integrals for each region in Section 4. In Section 5 we describe the methods used for their evaluation and we present analytical results, either as an expansion in the dimensional regulator or exactly in terms of generalised hypergeometric functions. We draw our conclusions in Section 6.

Double-real-virtual corrections to Higgs production
In this paper we consider the production of a Higgs boson (H) in association with two massless partons (k and l) via the scattering of two massless partons in the initial state (i and j), (2.1) Here the partons i, j, k and l can be gluons or massless quarks of N F different flavours, which are not directly coupled to the Higgs boson. We work in the large top-mass limit, and we assume the gluons to couple directly to the Higgs boson via the effective operator At N 3 LO all of these partonic processes enter the inclusive cross section in the large topmass limit through the interference of the relevant tree and one-loop amplitudes. To be more concrete, let us write the contribution to the cross section as where N ij;kl contains all averaging and symmetry factors and s ≡ s 12 denotes the square of the partonic center of mass energy. The mass of the Higgs boson will be denoted by m H , and for later convenience we introduce the following notations, and The dimensionless (coefficient) function C (1) ij→klH is defined as If we set the scale introduced by dimensional regularisation µ equal to m H (a convention which we shall adopt throughout this paper), then the coefficient function only depends on z and the dimensional regulator ǫ, related to the dimension of space-time D via ǫ = ij→Hkl denotes the L-loop amplitude for the the process ij → Hkl and the differential three-particle phase-space measure is defined as where δ + (q 2 − m 2 ) ≡ δ(q 2 − m 2 ) Θ(q 0 ), with Θ the Heaviside step function. In the rest of this paper we present methods to compute the coefficient function C ij→klH (z) as an expansion around threshold.

Threshold expansion for real emissions at one loop
Our goal is to obtain the expansion of the coefficient function C (1) ij→klH (z, ǫ) close to z = 1, where all the emitted final-state partons are soft. In other words, we want to expand eq. (2.7) in the soft momenta q 3 and q 4 . To be more concrete, let us define a set of rescaled momenta p i as After this rescaling both the phase-space measure and the matrix element in eq. (2.7) depend on the scaling parameterz. It will be convenient to introduce a notation for how a given quantity scales withz, e.g., we will write q 1 ∼ 1 and q 3 ∼z to refer to the scaling behaviour defined by eq. (3.1). Our goal is to expand both the phase-space measure and the amplitudes in this scaling parameter. In the rest of this section we discuss the expansion of each of these objects in turn. Let us start by discussing the threshold expansion of the phase-space measure. Changing variables to the rescaled momenta (3.1), it is obvious that the only quantity in eq. (2.8) that depends non-trivially on the scaling parameter is the on-shell condition for the Higgs boson, At leading order in the threshold expansion we can simply ignore the term linear inz appearing inside the δ function. If we want to study the subleading terms in the expansion, however, we need to take into account this term, i.e., we need to expand the δ function around threshold. This can be achieved using reverse unitarity [3,20] to interpret phase-space integrals as Feynman integrals with cut-propagators, i.e., where some of the propagators have been replaced by on-shell δ-functions, The subscript c is a reminder that this propagator is cut. Cut propagators can be differentiated just like normal propagators, but satisfy the extra condition The soft expansion of the phase-space measure therefore reads, (3.7) Note that this procedure necessarily introduces cut propagators raised to higher powers in the subleading terms in the expansion. In Section 4 we will see that integrals involving these additional powers of cut propagators can always be reduced to the case n = 0. As a result, all phase-space integrations will be performed against the soft phase-space measure defined by Let us now turn to the threshold expansion of the integrand. This requires expanding inz both the tree-level and one-loop amplitudes appearing inside the matrix element in eq. (2.7). The expansion of the tree-level amplitude can easily be obtained by introducing the rescaled momenta (3.1) and expanding the resulting rational function inz.
The expansion of the one-loop amplitude, however, is more subtle. In order to understand why, let us mention already now that the coefficient function is not meromorphic at z = 1, but it can nevertheless be written in the form where the functions C (1,r) ij→klH (z) for r ∈ {h, c, s} are holomorphic at z = 1 and therefore admit a Taylor expansion around this point, (3.10) A naive expansion of the loop integrand inz fails to reproduce this analytic structure, because the expansion of the loop integrand gives rise to a holomorphic function, and the phase-space measure can only account for the term proportional toz −4ǫ (see eq. (3.6)). In other words, the threshold expansion does not commute with the loop integration. Indeed, it is well-known that for some kinematic limits the loop integration does not commute with expansions of the loop integral, and hence a naive Taylor expansion of the loop integrand is bound to fail. The strategy of regions instead allows one to recover the correct expansion of the integral [22]. In this approach one sums up the expansions around all loop momentum regions which can be identified to lead to singularities in the limit z → 1. In ref. [8] it was argued that in the present case these different 'regions' can be classified as 'hard' (h), 'soft' (s) and 'collinear' (c), depending on whether the loop momentum is hard, soft or collinear to one of the hard partons. The three different terms in eq. (3.9) are then associated with these three regions. The collinear region is further decomposed into two different collinear regions, 'collinear-1' (c 1 ) and 'collinear-2' (c 2 ), where the loop momentum is collinear to p 1 or p 2 . Each region is associated with a specific scaling of the (components of the) loop momentum k, similar to the rescaling (3.1) used to define the threshold expansion of the phase-space measure. In order to define these scalings, we parametrise the momentum flowing through a suitable propagator by and the measure can then be written in this parametrisation as The different regions are then associated with the following scalings inz, The overall scaling of a given region in eq. (3.9) can then be understood by combining the scalings of phase-space measure, eq. (3.6), and loop momentum measure in the different regions.
At this stage we have to address a critical point in this procedure. The scalings we have just defined are not invariant under the shift symmetry of the loop integral, i.e., the rescaling can only be applied once an adequate shift of the loop momentum has been identified. We solve this issue in the following way: first, we note that, up to swaps of the external momenta p 3 , p 4 , all the loop integrations appearing inside the one-loop amplitude can be mapped to one of the following four pentagon topologies, Using the loop momentum scalings defined in the previous paragraph we have checked explicitly that the first few terms of thez expansions, obtained by summing up the expansions of non-vanishing hard, collinear and soft regions of various pentagons, boxes, triangles and bubbles defined through eq. (3.14), reproduce the same expansions as one can obtain from an expansion in Feynman parameter space by the methods of ref. [24]. Combined with the expansion of the phase-space measure discussed at the beginning of this section, we have therefore obtained a machinery to compute the threshold expansion of the coefficient functions C (1) ij→klH (z, ǫ). Let us conclude this section by making a technical comment about the expansion in the collinear regions. The procedure outlined above will effectively lead to an expansion in √z in the collinear region, because the transverse components of the loop momentum scale like √z . In other words, there seems to be a contradiction to the statement that the coefficient function in the collinear region, C (1,c) ij→klH (z, ǫ), is holomorphic at z = 1. In the following we show that this coefficient function in the collinear region is indeed holomorphic at z = 1, and therefore all terms proportional to powersz n/2 , with n odd, must vanish. In order to prove this statement, consider an integral of the form, where F is a Lorentz invariant function which only depends on the loop momentum k through the scalar products k 2 , k · p 1 and k · p 2 . Using the scaling properties of the collinear regions and the relations of eq. (3.12), it can be seen that all integrals which appear in thē z expansions of the collinear regions of the four pentagon topologies defined in eq. (3.14) indeed fall into this category. As the only source of √z in the expansion are the transverse components of the loop momentum k ⊥ , all terms in the expansion proportional toz n/2 with n odd must be proportional to an integral of the type (3.16). A sufficient condition to prove that C (1,r) ij→klH (z, ǫ) is holomorphic at z = 1 is therefore that all integrals of the type (3.16) vanish for odd values of n. In Appendix A we prove that this is indeed the case. More precisely, we prove the following result, where is the metric on the space transverse to the vectors p 1 and p 2 , and We note that eqs. (3.17) and (3.11) can be used efficiently to reduce any one-loop integral, which appears in the collinear region, containing arbitrary powers of the Lorentz invariants k·p 3 and k·p 4 in the numerator, to integrals containing only the Lorentz invariants k·p 1 , k·p 2 and k 2 in numerator and denominator.

IBP reduction and master integrals
In the previous section we argued that the threshold expansion of the coefficient function C (1) ij→klH receives contributions from hard, soft and collinear regions, see eq. (3.9). The integrals appearing in the expansion are, however, not independent, but can be reduced to a small set of master integrals using integration-by-parts (IBP) identities [21]. While IBP identities have been introduced for Feynman integrals rather than phase-space integrals, we can use reverse unitarity to interpret phase-space integrals as Feynman integrals with cut-propagators, to which IBP identities are known to apply. Here we go one step further, and we apply IBP identities combined with reverse-unitarity to the individual hard, soft and collinear regions. As a result, we obtain for each region a small set of hard, soft or collinear master integrals in terms of which the coefficients in the threshold expansion can be expressed. Independently of the region, we will see that the master integrals will fall into the class of so-called soft integrals, i.e., integrals with respect to the soft phase-space measure (3.8) of a function that is independently homogeneous under a rescaling of the initial momenta p 1 and p 2 , as well as under a simultaneous rescaling of all the final-state soft momenta. In the remainder of this section we review the IBP reduction in each region, and we present the analytic results for the master integrals in each region. Details about the computation of the master integrals will be given in Section 5.

The hard region
We start by discussing the master integrals coming from the hard region. Since the loop momentum is hard, we can immediately expand in the soft final state momenta under the integral sign, and the IBP reduction for the combined loop and phase-space integral follows the same lines as for the purely real soft emission discussed in ref. [15]. We find that the whole contribution from the hard region can be expressed in terms of only two master integrals, The double line denotes the Higgs boson, and the dashed line represents the phase-space cut. All other internal uncut lines are scalar propagators. Moreover, Bub(p 2 ) denotes the usual one-loop bubble integral, where we defined the usual loop factor The computation of the master integrals will be discussed in detail in Section 5. Note that Bub(s) = Bub(m 2 H ) + O(z), and so the hard master integrals can be written as a bubble integral Bub(m 2 H ) multiplying one of the two soft master integrals appearing in the double-real soft corrections at NNLO.

The soft region
Next, we discuss the contribution from the soft region, where k ∼z. It is straightforward to extend the methods of ref. [15] to the soft region at one loop. Indeed, since the loop momentum scales in the same way as the final-state soft momenta, IBP reduction proceeds in the same way as for purely real soft emissions at tree-level, except that the loop momentum is not constrained to be on shell. At the end of this procedure, we find the following set of soft master integrals, , , , , , , , .
The number of propagators of the soft master integrals depending on the loop momenta ranges from two to five. Since we are working in the soft region, the loop integrals are not the complete momentum space integrals, but they correspond to the eikonal approximation. Note that the eikonal approximation of the loop integrals depends on the specific orientation with which it enters the soft phase-space integrals, or equivalently, which invariants become soft. We have explicitly evaluated all loop integrals in the eikonal approximation by determining the scaling of the Feynman parameters in the limit using the package asy.m [24]. We find that in all cases, except for one pentagon integral, the remaining parametric integrations can be perfomed in closed form to all orders in the dimensional regulator ǫ. For the remaining pentagon we are able to obtain Mellin-Barnes integral representation valid to all orders in ǫ. In the following we give a brief summary of the results for the soft virtual integrals. The bubble type integrals in M S 3 and M S 7 can just be computed as in eq. (4.1), and will not be discussed any further. Note that we present the results in the Euclidean region where all invariants are negative, and the analytic continuation to the scattering region is given by Note also that the loop integrals appearing inside the master integrals have the correct homogeneity properties to turn the master integrals into soft integrals. The computation of the master integrals will be discussed in Section 5.
Soft Triangle Integrals. The soft master integrals M S 2 and M S 8 contain a virtual integral with three propagators. −p 1234 The integral corresponds to a completely off-shell triangle where two of the virtualities of the external momenta (s 24 = p 2 24 and s 13 = p 2 13 ) scale asz.
Soft Box Integrals. We observe that there are four configurations of loop integrals with four propagators. Two of them arise from box master integrals with one off-shell leg and two from configurations with two off-shell legs.
1. The loop integral appearing in the definition of M S 1 is given by 2. The soft master integrals M S 4 , M S 11 and M S 12 contain the loop integral where 2 F 1 denotes Gauss' hypergeometric function, 3. In M S 5 and M S 10 we find a limit of a two-mass-easy box integral. 4. The four propagator loop integral appearing in master integral M S 6 is a limit of a two-mass-hard box integral.
The loop integral in M S 9 is related to the integral above via the permutation of external momenta (p 1 → p 2 , p 2 → p 1 , p 4 → p 3 , p 3 → p 4 ).
Soft Pentagon Integrals. We also need the soft limits of two configurations of the one-loop pentagon with one off-shell leg.
1. Unlike the other one-loop soft integrals we require, the pentagon integral appearing in M S 13 does not permit a closed expression in terms of simple 2 F 1 hypergeometric functions. Nonetheless, we are able to derive a compact Mellin-Barnes representation.

The pentagon integral required for M S
14 is given by

The collinear region
In this section we discuss the reduction to master integrals in the collinear region. Unlike the soft and the hard regions discussed in the previous sections, the phase-space integrals in the collinear region cannot be reduced to master integrals using standard techniques. This is due to the appearance of propagators which are non-linear in the combined phase space and loop Mandelstam invariants. Effectively, the IBP relations do not close on certain topologies, and we are not aware of a method capable of dealing with IBP identities relating integral across different topologies. However, as we already pointed out in Section 3, all the one-loop integrals appearing in the collinear region have the property that their denominators only depend on the loop momentum k through the Lorentz invariants k 2 , k·p 1 and k · p 2 . This brings as a consequence that there exist relations among the denominators, which can be used to find partial fraction identities to reduce all pentagons and boxes to a number of different triangle topologies. After integration over the loop momentum these triangles can be identified as linear combinations of only two types of bubbles. As a consequence, we can write the contribution from the collinear region in the form where A and B are expanded as power series inz, whose coefficients are rational functions of the Mandelstam invariants and the dimensional regulator ǫ. We dropped the real part of the loop integrals, because they are real for s 13 , s 14 < 0. We emphasise that eq. (4.5) is true to all orders in the expansion parameterz and can be proved simply by inspecting the propagators in the collinear region. Next, we express bubble integrals in terms of tadpoles via the identity Since the single denominator in eq. (4.6) is clearly linear in all Lorentz invariants, this representation allows for a straightforward and efficient application of the IBP reduction technique to the combined soft phase space and (now trivial) one loop integral. We find that in the collinear region we can reduce all integrals to just four soft master integrals,  where the dotted lines represent numerator factors. We stress that the master integrals appearing in the collinear region are again soft integrals, and hence they can be computed using the same techniques as the master integrals in the soft and hard regions. This will be discussed in Section 5.

Evaluation of soft one-loop integrals
In the previous section we have seen that, independently of the region they originate from, the master integrals appearing in the threshold expansion are always soft integrals, i.e., they all take the form where dΦ S 3 denotes the soft phase-space measure of eq. (3.8), and the integrand is homogeneous with respect to a rescaling of the soft momenta as well as of the initial-state momenta p 1 and p 2 , , F I i (p 1 , λ p 2 , p 3 , p 4 ; ǫ) = λ β iI F I i (p 1 , p 2 , p 3 , p 4 ; ǫ) , F I i (p 1 , p 2 , λ p 3 , λ p 4 ; ǫ) = λ γ iI F I i (p 1 , p 2 , p 3 , p 4 ; ǫ) , for some α iI , β iI , γ iI .
In this section we present methods to evaluate soft integrals analytically. We discuss two methods: the first is based on the derivation of a Mellin-Barnes representation for soft integrals, while the second exploits a particular phase-space factorisation as well as the homogeneity of the soft integrals to derive a parametric integral representation for soft integrals. For an additional method, using differential equation techniques, we refer to ref. [23]. We note that although we concentrate exclusively on the case of double-real emissions at one loop, all these techniques can in principle be generalised to arbitrary numbers of loops and legs in an obvious way.

Mellin-Barnes representations for soft integrals
In this section we describe a general method to write a given soft integral as a (possibly multi-fold) Mellin-Barnes (MB) integral with poles at integer values of the integration variables. In ref. [15] such a procedure was introduced for soft integrals of purely real emissions at tree-level. In this section we briefly review the procedure of ref. [15], and in the end we argue that it can easily be extended beyond tree-level. As the procedure is essentially identical to the purely real case, we will be brief and refer to ref. [15] for details.
The procedure of ref. [15] starts from the observation that at tree level one may assume without loss of generality that the integrand F is a product of powers of two-particle Mandelstam invariants. Indeed, this can always be achieved by replacing every sum of two-particle invariants in the denominator by its MB representation, using the well-known formula 1 where the contour separates the poles at z = n from those at z = −λ − n, n ∈ N. Next, we parametrise the soft phase space using the energies and the angles of the soft momenta in the center-of-mass frame of the initial-state system. We write, with s 12 = 1, where β i is the four-velocity in the direction of p i . The soft phase-space measure becomes, where Ω D−1 i parametrises the solid angle of the soft momentum p i . Due to our assumption that the integrand is a product of powers of two-particle Mandelstam invariants, the integration over the energies is simply a Beta function. The remaining angular integrals can be written as MB integrals with poles at most at integer values of the integration variables [25,26].
If we follow this procedure, every tree-level soft integral can be written as a multifold MB integral with poles at integer values of the integration variables. If the integrand contains loop integrals which evaluate to complicated special functions, this claim is no longer necessarily true. It does, however, stay true if the loop integral itself admits an MB representation of the same type. In our case, the integrand contains at most hypergeometric functions, which admit the MB representation An exception is one pentagon integral, which we did not express in terms of simple hypergeometric functions, but it still admits a multifold MB representation of a similar type, see eq. (4.5). We can therefore obtain an MB representation of the desired type for all soft integrals considered in this paper. The MB integrals obtained in this way can be evaluated using standard techniques. In some cases it is possible to perform all the MB integration in closed form without expanding in ǫ, and one obtains generalised hypergeometric functions, Whenever we are not able to obtain a closed expression in terms of hypergeometric functions, we resolve the poles in ǫ using standard techniques [27]. The result is a Laurent series in ǫ whose coefficients are MB integrals whose contours are straight vertical line. In all cases these integrals can be evaluated numerically in a fast and efficient way. Alternatively, one can close the integration contours at infinity and sum up the residues of the poles of the Gamma functions in terms of nested harmonic sums [28] that evaluate to multiple zeta values.

Soft integrals from phase-space factorisation
We present in this section an alternative way to compute the soft integrals of Section 4, based on a factorisation of the phase space that separates the soft part of the phase space from the phase space for the emission of the Higgs boson. We start by writing the phase space for the production of a Higgs bosons in association with two massless partons as a convolution dΦ 3 = dµ 2 2π dΦ 2 (m 2 H , µ 2 ; p 12 ) dΦ 2 (0, 0; Q) , (5.8) where dΦ 2 (m 2 1 , m 2 2 ; q) denotes the phase-space measure for the decay of a heavy state with momentum q into two particles with masses m 1 and m 2 , Note that for now we work with the full phase-space measure, and we expand inz at the end. A soft integral of the type (5.1) can then be written as where we defined F I i (p 1 , p 2 , Q) = dΦ 2 (0, 0; Q) F I i (p 1 , p 2 , p 3 , p 4 ; ǫ) . (5.11) Let us first concentrate on the computation of the integral F I i . The homogeneity of the integrand of the original soft integral, eq. (5.2), combined with the fact that the integration measure does not depend on p 1 and p 2 , implies that F I i is homogeneous under a rescaling of any of its arguments, Lorentz invariance then implies that the non-trivial functional dependence of F I i can only be through the 'cross ratio' .
Without loss of generality, we may write We can also give a kinematical meaning to the cross ratio u: it is related to the angle θ 12 between p 1 and p 2 in the rest frame of Q, u = 1 − cos θ 12 2 .
withβ i = (1, − n i ). Here we used the fact that in the rest frame of Q the final-state massless particles are back-to-back and their energy is Q 2 /2. The remaining integral over the solid angle can in all cases be performed using the formula [25,26] dΩ D−1 (5.18) The scalar product appearing in the right hand side can take the following values, β 1 · β 2 = β 1 · β 2 = 1 − cos θ 12 = 2u , β 1 · β 2 = β 1 · β 2 = 1 + cos θ 12 = 2(1 − u) . (5.19) Note that in some cases the integrand may contain hypergeometric functions depending on angular variables. We can reduce the problem to the angular integral (5.18) by introducing an MB representation for these hypergeometric functions. Eventually, we can in this way explicitly determine the function f I i (u) in eq. (5.14). We stress that although the computation was done in the rest frame of Q, the result is independent of the frame.
Inserting eq. (5.14) into eq. (5.10), we are left with the computation of the remaining phase-space integral in eq. (5.10). We write Q = α p 1 +β p 2 +Q ⊥ , with p 1 ·Q ⊥ = p 2 ·Q ⊥ = 0 and we work in the center-of-mass frame of the collision, and the on-shell conditions for Q and the Higgs boson enforce α + β = 1 + O(z) and Q 2 ⊥ = α(1 − α) − µ 2 . Note that this implies 0 < α < 1 and 0 < µ 2 < α(1 − α). In terms of this parametrisation the invariants are 20) and the phase-space measure can be written as where Ω D−2 Q parametrises the solid angle of Q in the center-of-mass frame. Putting everything together we see that the integrals over α and over the solid angle are trivial, and so we get, with δ = (γ − α − β)/2, We obtain in this way a representation for the soft integral as a simple integral over the function f I i . This last integral is usually easy to carry out. In many cases it can be performed using the Euler-integral representation of the generalised hypergeometric function, a 2 ; b 1 ; u z) .

(5.23)
In those cases where the function f I i is more complicated, we can either insert an MB representation for it, or alternatively perform the integral over u order by order in ǫ in terms of harmonic polylogarithms [29].

Analytic results for the master integrals
In this section we present the analytic results for all the master integrals introduced in Section 4. As only the real part of the master integrals enters the final result for the cross section, we only present result for the real part, which we normalise according to where c Γ is defined in eq. (4.2) and Φ S 3 (ǫ) denotes the soft phase-space volume, We have computed all the master integrals of Section 4 using the two different approaches described in the previous section, and we found complete agreement between the approaches. With the exception of the pentagon integrals M S 13 and M S 14 we present the results for the master integrals as a Laurent expansion in the dimensional regulator up to terms of transcendental weight six for the soft and the hard regions and up to transcendental weight five for the collinear region. For M S 13 and M S 14 we were only able to obtain the Laurent expansion up to terms of weight three and five respectively, which is however sufficient to compute the Higgs boson cross section up to finite terms 2 . Whenever we were able to do so, we also include results for the master integrals valid to all orders in ǫ in terms of generalised hypergeometric functions, which can easily be expanded in ǫ using the HypExp package [30]. Note that all integrals presented in this section satisfy recursion relations with respect to the space-time dimension [31]. We have checked that our results satisfy these dimensional recurrence relations. techniques to compute the coefficients in the expansion, in principle to any desired order, and to express the result in terms of a small set of soft master integrals. These results are the missing pieces which went into the computation of the inclusive gluon-fusion Higgs production cross section as an expansion around threshold [6,16,19].

Analytic results for the soft region
Our main tool to reduce the coefficients appearing in the threshold expansion to a linear combination of soft master integrals is reverse unitarity, which allows one to map phase-space integrals to cuts of Feynman integral. We perform this expansion separately for the phase-space measure and for the interference diagrams, and we observe that the phase-space integrals always reduce to integrals against the soft phase-space measure. The one-loop matrix elements are expanded in the soft limit using the strategy of regions, and the relevant regions are identified with the regions where the loop momentum is either hard, soft or collinear to one of the initial-state momenta. In each region, we combine the expanded interference diagrams with the corresponding phase-space measure and use IBP identities to reduce them to soft master integrals, independently of the region.
The soft master integrals themselves are evaluated using two different approaches. The first method allows to derive a Mellin-Barnes representation for the soft integrals with poles at most at integer locations, provided that a similar Mellin-Barnes representation for the loop integration can be obtained. In our case most of the Mellin-Barnes integrals can be done in closed form to all orders in the dimensional regulator. In the remaining cases we were able to obtain a representation of the integral as a Laurent expansion in the dimensional regulator. The second method builds upon a specific factorisation of the phase space and exploits the knowledge of the scaling behavior of the integral with the external momenta to arrive at a one-fold parametric integral representation for the soft integral, which can be performed using modern integration techniques.
The results for the soft master integrals presented in this paper are sufficient to obtain, at least, the first 37 terms in the threshold expansion of the N 3 LO gluon-fusion cross section, and conjecturally they provide the full set of boundary conditions to compute all the phasespace master integrals in general kinematics appearing in inclusive N 3 LO cross sections for the production of a heavy colorless state. We therefore anticipate that the results of this paper will have a substantial impact on future results for hadron collider cross sections at N 3 LO, e.g., Drell-Yan production or Higgs production via bottom-quark fusion.