Real-Virtual contributions to the inclusive Higgs cross-section at N3LO

We compute the contributions to the N3LO inclusive Higgs boson cross-section from the square of one-loop amplitudes with a Higgs boson and three QCD partons as external states. Our result is a Taylor expansion in the dimensional regulator epsilon, where the coefficients of the expansion are analytic functions of the ratio of the Higgs boson mass and the partonic center of mass energy and they are valid for arbitrary values of this ratio. We also perform a threshold expansion around the limit of soft-parton radiation in the final state. The expressions for the coefficients of the threshold expansion are valid for arbitrary values of the dimension. As a by-product of the threshold expansion calculation, we have developed a soft expansion method at the integrand level by identifying the relevant soft and collinear regions for the loop-momentum.


Introduction
The theoretical prediction of the Standard Model inclusive Higgs boson cross-section is an important reference for the experimental verification of the model at the energies probed by the Large Hadron Collider. The uncertainty which is associated with the truncation of the perturbative expansion at NNLO is currently estimated at the order of ±10% (see, for example, ref. [1]). With the accumulation of more data by ATLAS and CMS, this perturbative uncertainty will be ultimately one of the largest systematic uncertainties entering the extraction of the coupling strengths of Higgs boson interactions. A calculation of the N 3 LO cross-section may reduce it to about ±5% [2].
A computation of the N 3 LO cross-section is therefore very much desired. Many important steps have been taken towards this objective. The three-loop gg → h amplitude has been computed in refs. [3][4][5]. Contributions which are associated with collinear/ultraviolet counter-terms and the partonic cross-sections at lower orders have been computed in refs. [2,[6][7][8]. N 3 LO corrections due to processes with triple real radiation were computed in ref. [9] as a threshold expansion around the limit of soft-parton emissions. An important ingredient for the threshold expansion of the Higgs boson cross-section with both real and virtual radiation, the two-loop soft current, was presented in refs. [10,11].
In this article, we focus on a different contribution to the N 3 LO cross-section, the integration over phase-space of the squared one-loop amplitudes for partonic processes for Higgs production in association with a quark or gluon in the final-state. We denote these squared real-virtual contributions as (RV) 2 . We have been able to perform many independent calculations for the (RV) 2 cross-sections by employing different methods.
First, we have two different methods that are suited for a threshold expansion. In method (Ia), we reduce the one-loop amplitudes to bubble and box master integrals and find appropriate representations of the box master integrals which allow for a trivial integration over phase-space order by order in the threshold expansion. Method (Ib) is also a threshold expansion method and relies on expanding the amplitudes at the integrand level in the regions where particles in the loops are soft or collinear to an external particle. We have identified all such regions and we have proven that the expansion by regions yields the same result as our first method for an arbitrary order in the threshold expansion. The coefficients of the threshold expansion are exact in the dimensional regulator ǫ = 2 − d 2 . Methods (II) and (III) yield expressions of the partonic cross-sections which are valid for arbitrary values of the partonic center of mass energy. Method (II) treats the combined loop and phase-space integrations as a three-loop forward scattering amplitude, which is reduced to master integrals with the reverse-unitarity method [12][13][14][15]. The master integrals themselves are evaluated by solving the differential equations they satisfy in terms of generalized hypergeometric functions. The solution of the differential equations has been achieved by observing recursion patterns for the hypergeometric integrals, similar to the ones for in multiple polylogarithms [16]. Due to the hypergeometric representation of our master integrals a threshold expansion is easily performed, thus providing an independent verification of the results of Methods (Ia) and (Ib). To obtain the master integrals with full kinematical dependence we solve the differential equations as an ǫ expansion in terms of harmonic polylogarithms. Method (III) identifies counterterms for combined loop and phase-space integrals and proceeds with an expansion in ǫ followed by a direct integration over phase-space variables. This is made possible by embedding the classical polylogarithms that arise in the ǫ expansion into a larger space of multiple polylogarithms where the integration over phase space is trivial. This embedding is achieved by deriving the required functional equations using symbols [17][18][19][20][21] and the Hopf algebra structure of multiple polylogarithms [22][23][24].
This article is organized as follows. In Section 2 we introduce the partonic crosssections we are computing and define our notation. In Section 3 we present our results for the partonic cross-sections. In Section 4 we detail the methods employed for the purposes of our calculation. We present our conclusions in Section 5.

The real-virtual squared cross section
We consider the partonic-processes for associated Higgs production, where g, q,q are symbols for gluon, quark and anti-quark partons correspondingly and H for the Higgs boson. The brackets refer to the momenta of the particles. We define the kinematic invariants, s ≡ 2p 1 · p 2 + i0, t ≡ 2p 2 · p 3 − i0, u ≡ 2p 1 · p 3 − i0, where we indicated explicitly the small imaginary parts carried by them. The partonic cross-sections for these processes are given by where X ∈ {gg → Hg, gq → Hq, qq → Hg} labels the different subprocesses 1 and the sum symbol denotes a summation over colors and polarizations of the initial and final state particles. We work in conventional regularization in d = 4 − 2ǫ for both the phase space and the matrix element. The d-dimensional phase-space measure is given by N X denote the averaging over initial state spins and colors in d dimensions, 5) with N and V ≡ N 2 − 1 the number of quark and gluon colors respectively. In the following it will be convenient to parametrize the invariants as s = M 2 h z , t = s δ λ, u = s δ (1 − λ) , (2.6) with δ = 1 − z. Note that a physical scattering process corresponds to s > 0 and 0 < z, λ < 1, and the limit δ → 1 corresponds to the threshold region where the additional final state parton is soft. Using the parametrization (2.6) the phase-space measure becomes where Θ(x) is the Heaviside step function. Equation (2.3) then reads In the rest of this paper we compute in an effective theory of the Standard Model where the top-quark is decoupled. The effective unrenormalized Lagrangian reads: where the first term corresponds to an effective QCD Lagrangian with N f = 5 flavors. The Wilson coefficient C can be cast as a function of the QCD coupling, the bare heavy-quark masses and the Higgs field vacuum expectation value [25][26][27][28][29]. Performing a loop-expansion of the amplitudes A X = ∞ j=0 A (j) X in the effective theory with j being the number of loops, we have: (2.10) The first two terms of the above expansion enter the already known inclusive Higgs boson cross-section through NNLO [12,[30][31][32][33][34][35]. The third term in square brackets contributes to the N 3 LO coefficient. In this article, we compute the part of the partonic cross-section due to the square of the one-loop amplitudes, namely, (2.11)

Results
The computation of the cross-sections (2.11) for the different subprocesses is the main subject of this paper. We evaluated the phase-space and loop integrals in different ways that are detailed in Section 4. In this section, we summarize first our main findings. We find for the partonic cross-sections: where g 2 s ≡ 4πα s , and we have In eq. (3.1) we have introduced the quantity The function Σ gg (z; ǫ) has a pole as z → 1. This pole constitutes the soft singularity of the gluon initiated cross-section and it will only be remedied by integrating the partonic cross-section with parton distribution functions. We separate this singular part manifestly from the remainder and write: where the z → 1 singular part is given by The cross-sections of the qq and qg channels are regular in the limit z → 1: (3.6) We have calculated the regular functions Σ reg X (x; ǫ) as an expansion in the dimensional regulator ǫ through O(ǫ 0 ). The expressions are composed of multiple polylogarithms (MPLs) [16] with weights up to five. Specifically, the MPLs are defined as iterated integrals via G(a 1 , . . . , a n ; z) = z 0 dt t − a 1 G(a 2 , . . . , a n ; t), where G(z) = 1 and z, a i ∈ C and the weight is the number of indices a i . All our MPLs have argument z and all a i are elements of the set {−1, 0, 1}. The last property allows to relate our MPLs to a more restrictive class of functions, the so-called harmonic polylogarithms (HPLs) [36]. The relation is given explicitly by H(a 1 , . . . , a n ; z) = (−1) σ(a 1 ,...,an) G(a 1 , . . . , a n ; z), (3.8) where σ(a 1 , . . . , a n ) denotes the number of elements in the set {a 1 , . . . , a n } that are equal to 1. Software to evaluate HPLs numerically in a fast and accurate way (at least up to weight four) is publicly available [37][38][39][40][41]. Due to the magnitude of the expressions obtained for the functions Σ reg X (x; ǫ) we refrain from stating them here explicitly and make them publicly available together with the arXiv submission of this article and on the web-page [42]. In the file results.txt we present the formulae for the Σ reg X formatted in Maple input form. Given that the triple-real contributions to the inclusive Higgs boson cross-section at N 3 LO [9] have only been computed as an expansion around z = 1 − δ → 1, we also provide the same expansion for the (RV) 2 cross-section functions Σ X (z, ǫ) of this article. We have discovered a characteristic structure for Σ X (z, ǫ) in the δ → 0 limit; all logarithmic contributions of the form log δ exponentiate into factors of δ −aǫ , where a is an integer in the interval [2,6]. Namely, where the functions η (a) X (δ; ǫ) are meromorphic functions of δ. We further decompose for a = 4, while for a = 4 we write, (3.12) The structure of Eqs. (3.9)-(3.11) originates from loop integrations over distinct kinematic configurations, where the loop momentum can be either soft (s), or collinear to the first incoming parton (c 1 ), or collinear to the second incoming parton (c 2 ), or, otherwise, hard (h). Every term in the squared amplitude may be thought of as associated to a product of two regions r 1 and r 2 , r i ∈ {s, c 1 , c 2 , h}, one from the amplitude itself and the other from its complex conjugate. In the following we denote the contribution from such a product of regions by (r 1 , r 2 ). The coefficients φ X (δ; ǫ) are in one-to-one correspondence with these regions. The correspondence is explicitly given by (3.13) A derivation of the above decomposition will be given in Section 4. The analytic form of the double soft (s, s) termsη (6;1) X (δ; ǫ) is rather simple. We find: .
The remainingη a;i X terms are more complicated combinations of generalized hypergeometric functions, which can be readily cast as a Laurent series in δ

Methods
In this section, we discuss how we evaluated the loop and phase-space integrals that contribute to the real-virtual squared cross-sections (3.1). We employed various methods that each have their own strengths and weaknesses. We checked that we obtain consistent results when comparing the different approaches. These methods are: 1. Threshold expansion of the cross section: In this approach, we derive a representation of the cross-section as an expansion close to threshold where δ → 0. We first expand the loop amplitude in the limit where the final state parton is soft, and then perform the phase-space integration order-by-order in the expansion. The threshold expansion of the loop amplitude is obtained in two different ways: first by finding a suitable representation in terms of convergent hypergeometric functions within the entire phase-space, and, second, by expanding around the relevant soft, collinear and hard regions of the loop momentum. We will detail our threshold expansion techniques in Section 4.1.
2. Differential equations for master integrals: Using a duality of loop and phasespace integrals we reduce them simultaneously to a minimal set of master integrals. The master integrals satisfy differential equations that can be solved in two ways: either order-by-order in dimensional regularization in terms of harmonic polylogarithms, or in terms of generalized hypergeometric functions. The boundary conditions for the differential equations are obtained by matching to the leading term of the threshold expansion, which we compute with one of our threshold expansion methods mentioned above. We present our approach based on differential equations in Section 4.2.
3. Direct integration using multiple polylogarithms: It is possible to derive analytic results of the loop integrals entering our amplitudes in terms of polylogarithmic functions. These expressions are singular in soft and collinear limits of the phasespace and we render all integrals convergent by constructing appropriate counterterms. Then we perform the two-body phase-space integration by embedding the polylogarithmic functions into a larger class of multiple polylogarithms for which the integration is trivial. In the end, we recast the final result in terms of harmonic polylogarithms only. This method is explained in Section 4.3.

Expansion of the cross-section around threshold
We start by computing the one-loop amplitudes A 1 X with X ∈ {gg → Hg, qq → Hg, qg → Hq} in dimensional regularization and for arbitrary values of the regulator ǫ. We generate the Feynman diagrams with QGRAF [43] and perform the spin and color algebra using FORM [44]. Using the methods of refs. [45,46] for tensor integrals and well established reduction techniques for scalar integrals [47,48], the amplitudes are reduced to the one-loop scalar bubble and box master integrals, where the q i are considered light-like and ingoing and s ij = (q i + q j ) 2 .
In a next step, we construct the squared one-loop amplitudes and cast them in terms of three functions, in the form: Explicit expressions for the functions A qqg , A ggg , B ggg are given in Appendix C.
Threshold expansion using hypergeometric functions. Loop integrals in dimensional regularization can be expressed, to all orders in the dimensional regulator, as (generalized) hypergeometric functions. For example, the box integral 2 defined in eq. (4.1) admits the representation (see, e.g., ref. [49]), Our goal is to insert the parametrization (2.6) and then to perform the integration over λ term-by-term in the series representation of the hypergeometric functions. The result is a power series in δ, i.e., the desired expansion of the cross-section close to threshold.
While all the hypergeometric series in eq. (4.5) are convergent in the Euclidean region where s ij < 0, the one-loop amplitude in the physical region involves the functions Box(−t, −u, s), Box(s, −t, −u) and Box(s, −u, −t). It is easy to check that the corresponding hypergeometric series are no longer convergent in the physical scattering region. It is, however, always possible to analytically continue the 2 F 1 function such that arguments lie inside the unit disc, yielding another representation in terms of 2 F 1 functions. While this approach is adequate to find a meaningful expansion around ǫ in terms of polylogarithms, it does not allow one to find (convergent) hypergeometric series expansions around δ = 0. Instead, one needs at least a double sum representation to achieve this task. It turns out that such a representation is known in the literature [49], The corresponding result for Box(s, −u, −t) is obtained from Box(s, −t, −u) by exchanging t and u. The generalized Kampé de Fériet function S 1 is defined as 8) and the Appel function F 2 is defined as Using these expressions for the box functions, we can easily exchange the phase-space integration and the infinite summations, and all the integrals can be performed in terms of Euler's Beta function, Threshold expansion of hard, soft and collinear regions. It is possible to derive representations such as the ones of (4.6) and (4.7) with a more physical method, performing Taylor expansions around soft, collinear and hard regions of the integrand of loop integrals in momentum space. The method of expansions by regions [50] promises to hold in general, although its generality has only been stated as a conjecture and a verification of the validity of the approach is necessary in specific cases 3 . For the production of the Higgs boson near threshold, the partonic center of mass energy is close in value to the Higgs boson mass, and thus we have a kinematic variable which is small, δ = 1 − z ∼ 0. From eq. (2.6) we infer that the external momenta scale as For a particle propagating in the loop, we find four types of non-trivial scalings of its momentum k: where all propagators in the loop are off-shell, where the loop integrand is singular at the point k µ = 0, where the integrand has a singular surface as k µ ∝ p µ 1 , where the integrand has a singular surface as k µ ∝ p µ 2 .
In the above, the transverse momentum k ⊥ of the particle is defined via: A scaling of the loop momentum is called a region. In a given region, we can perform a systematic expansion of the integrand around δ = 0. This yields multiple new integrals which are simpler than the unexpanded integral. For some regions, we are able to compute analytically all (i.e., an infinite number of) terms of the expansion. For the remaining regions, we limit ourselves to a finite number of terms in the expansion and perform an algebraic reduction [47,48] (after expansion) to master integrals. The soft and collinear regions of our loop integrals correspond to the singular surfaces which solve the Landau equations [53] while the loop momentum scalings can be identified with the scalings of the coordinates which are normal to the singular surfaces [54]. In the following we discuss how we can reproduce the hypergeometric function representations given in eq. (4.6). We start by discussing the asymptotic expansion of Box(s, −t, −u), which we find convenient to parametrize as (4.14) We find that the full integral is reconstructed by two regions: After Taylor expanding the loop integrand in every region in the small variable δ, we can use integration-by-parts identities to reduce the coefficients of the Taylor expansion to a small set of master integrals. In the (c 2 )-region we find that all the coefficients are proportional to the one-loop bubble integral Bub (−t), and this region reconstructs, order by order in δ, the first term of the hypergeometric representation for Box(s, −t, −u) given eq. 4.7, and we have verified this statement explicitly up to O δ 10 . The (h)-region yields the second and last term of eq. (4.7). In this region, we have been able to calculate all terms in the expansion around δ = 0 with an analytic integration. We see that the sum of the (h) and (c 2 ) regions is equal to the correct expression for the one-loop box. All other soft and collinear regions are zero, as we can readily verify. Next, we turn to the asymptotic expansion of the Box(−t, −u, s), given by (4. 16) We find that the expression for Box(−t, −u, s) given in eq. (4.6) is reconstructed entirely from the following regions: 1. (s)−region, where the k ∼ δ, yielding the first term of eq. (4.6). It is interesting that the (s)−region consists of a single term without any subleading terms in the expansion in δ.
2. (c 1 )-region where the momentum k is collinear to p 1 . This region reconstructs the second term of eq. (4.6) as we have verified explicitly up to O δ 10 .
3. (c 2 )-region where the momentum k − p 3 is collinear to p 2 . This region reconstructs the third term of eq. (4.6) as we have verified explicitly up to O δ 10 .
4. (h)-region, where k is hard. This region reconstructs the last term of eq. (4.6) as we have verified explicitly up to O δ 10 .
All other soft and collinear regions are zero.
We have seen that an expansion in hard, soft and collinear regions yields series representations for the one-loop master integrals of the required amplitudes which converge in the entire phase-space, and thus we can immediately perform the phase-space integration in terms of Beta functions order-by-order in the expansion. While in our case the strategy of expansion by regions is only an alternative method for deriving the threshold expansions of Section 4.1, it can be the method of choice for the phase-space integration of more complicated one-loop amplitudes. Here we have presented expansions by regions at the level of master integrals. We would like to remark that such expansions can also be performed at the integrand of loop-amplitudes before any reduction to master integrals has taken place. Combined with the method of reverse-unitarity [9] we have a powerful algebraic technique for the simultaneous threshold expansion of integrals over loop and external momenta.

Reverse unitarity and differential equations
In this section we evaluate the real-virtual squared cross-sections using the reverse-unitarity approach [12][13][14][15]. Reverse unitarity establishes a duality between phase-space integrals and loop integrals. Specifically, on-shell and other phase-space constraints are dual to "cut" propagators A cut-propagator can be differentiated similarly to an ordinary propagator with respect to its momenta. It is therefore possible to derive integration-by-parts (IBP) identities [55,56] for phase-space integrals in the same way as for loop integrals. The only difference is an additional simplifying constraint that a cut-propagator raised to a negative power vanishes: In this approach, we are not obliged to perform a strictly sequential evaluation of the loop integrals in the amplitude followed by the nested phase-space integrals. Rather, we combine the two types of integrals into a single multiloop-like type of integration by introducing cutpropagators and then derive and solve IBP identities for the combined integrals. We solve the large system of IBP identities which are relevant for our calculation with the Gauss elimination algorithm of Laporta [47]. We have made an independent implementation of the algorithm in C++ using also the GiNaC library [57]. In comparison to AIR [48], which is a second reduction program used in this work, the C++ implementation is faster and more powerful, storing all identities in virtual memory rather than in the file system. All integrals that appear in the real-virtual squared cross section are reduced to linear combinations of 19 master integrals, which we choose as follows: Single solid lines represent scalar massless propagators. The phase-space integration is represented by the dashed line and the cut-propagators are the lines cut by the dashed line. The cut propagator of the Higgs boson is depicted by the double-line. Every master integral has a one-loop integral on the left-and a complex-conjugated one-loop integral on the right-hand side of the cut. In each side of the cut, we find scalar bubble, box or triangle integrals, where the latter is defined by with q 2 i = 0, p 2 i = 0 and (p 1 + p 2 ) 2 = 0. The scalar bubble and box integrals have been defined in (4.1). A comment is in order about the appearance of the triangle integrals in this approach, which seems to be at odds with the fact that in the expression of the oneloop amplitude presented Section 4.1 only bubble and box integrals appeared. Indeed, it is well-known that eq. (4.38) can be expressed as a linear combination of bubble integrals, (4.39) These relations however introduce new denominators which need to be taken into account in the reduction of the phase-space integrals. We therefore prefer not to use eq. (4.39), but we work directly with the triangle integrals instead.
To evaluate the master integrals we employ the method of differential equations [12,58,59]. Differentiating the corresponding cut-propagator with respect to the square of the Higgs mass, results in another phase-space integral. This new integral can again be reduced by IBP identities to our basis of master integrals. Proceeding in this way we obtain a system of linear first order differential equations for the master integrals, The system is triangular where y i (δ) only depends on master integrals that can be solved for independently of M i (δ). In other words, the system can be solved hierarchically, starting from the differential equations with vanishing or known functions y. Every time we solve such an equation, its solution serves to determine the y function of a next equation. In this way, at any stage of this procedure the y function is a linear combination of already evaluated master integrals that can be integrated in order to determine the integral M j (δ). The coefficients A ij (δ) are rational functions in δ and ǫ and have isolated singularities in δ only at δ = 0, 1, 2. The first step to solving this type of differential equation is to find a solution for the homogeneous part. The general homogeneous solution associated to the differential equation (4.42) is given by and is determined up to an integration constant M i (0). We determine this integration constant by calculating the soft limit of the master integral explicitly following the methods discussed in Section 4.1. We find that only 7 of our 19 master integrals have non-trivial boundary conditions. Interestingly, with our choice of basis of master integrals, the nontrivial boundary conditions are in one-to-one correspondence to the leading terms of the 7 regions of the soft expansion of the squared amplitude of eqs. (3.10)-(3.11). The non-trivial boundary conditions are: Only the real part of the boundary conditions is presented here, given that the imaginary part does not contribute to the cross-section.
Once the homogeneous solution is found we can compute a particular solution to the inhomogeneous equation by The full solution for the master integral is then given by We perform the integration in the equation above with two different approaches.
Solving differential equations as an expansion in ǫ. One well established strategy is to expand the differential equations in powers of the dimensional regulator [58,59]. After expanding the integral of (4.46) a solution is naturally given by iterated integrals leading to multiple polylogarithms [16] of the form G(a 1 , . . . , a n ; δ), with a i ∈ {0, 1, 2}. Expressing the functions in terms of the variable z = 1 − δ recasts the solutions in terms of more familiar harmonic polylogarithms [58].
Solving differential equations in terms of hypergeometric functions. The integrand of eq. (4.46) takes the form where c 1 , c 2 and c 3 are linear polynomials in ǫ. This structure reminds of the Euler-type integral representation of hypergeometric functions. Inspired by the large variety of techniques available for the solution of iterated integrals in terms of multiple polylogarithms [16] we define an iterated integral with integration kernel δ a−1 (1 − xδ) −b . The n th iterated integral is then recursively defined by F an,..., a 1 (x n , . . . , where we have abbreviated a i = a i b i . We find that these iterated integrals interpolate between multiple polylogarithms and hypergeometric functions. For example, in this framework the multiple polylogarithm is given by where 0 = 0 0 and 1 = 1 1 . The Gauss hypergeometric function is given by with v = we find an explicit sum representation for this type of iterated integrals.
F an,..., a 1 (x n , . . . , Equation (4.53) is valid whenever the sums are convergent. Further properties and derivations are discussed in more detail in Appendix A. The solution of differential equations using iterated integrals is illustrated with an example in Appendix B.
The iterated integrals defined in this section are a powerful tool and enable us to solve all 19 master integrals in terms of hypergeometric functions. The results is valid to all orders in ǫ. The iterated integrals can be written as multiple sums eq. (4.53) from which it is very convenient to extract a threshold expansion in δ.
Results for the master integrals. The master integrals that we have computed in this section are useful for the evaluation of any cross-section for a 2 → 1 process at N 3 LO. Due to the length of their expressions, we provide them in terms of harmonic polylogarithms in the form of the text file masters.txt enclosed with the arXiv submission of this work and on the web-page [42]. We have computed all 19 master integrals to all orders in ǫ in terms of hypergeometric functions and as an expansion in ǫ in terms of harmonic polylogarithms as described above.

Direct integration using multiple polylogarithms
We present here an alternate method to compute the (RV) 2 Higgs boson cross-section, based on subtraction terms. The phase-space integral over the squared amplitude can be written schematically as, (4.54) In this expression M i denote the one-loop master integrals and N i,j are rational functions, all of which depend on the invariants s 12 , s 23 and s 13 . Since the results for the required one-loop master integrals are known to all orders in ǫ [49], the integrals are well defined in dimensional regularization. Our goal is to expand the integrals in ǫ under the integration sign and to perform the integration order by order in ǫ. After expansion, however, the integrals may develop soft and collinear divergencies. The strategy is to subtract the singular limits of the integrand before expansion, and to perform the remaining (finite) integration in terms of multiple polylogarithms. The construction of the counterterms that render the integration finite proceeds in two steps. First, we analytically continue all the hypergeometric functions that appear in the all order expressions of the one-loop master integrals such that they are convergent in the whole phase-space. This is achieved by using the well-known identities Second, the soft and collinear counterterms are easily constructed by expanding the integrand around the collinear limits, i.e., s 13 → 0 or s 23 → 0. The counterterms can be trivially integrated to all orders in the dimensional regulator in terms of Γ functions. At the end of this procedure we are left with finite one-dimensional integrals. We expand the hypergeometric functions appearing in the integrand in ǫ using HypExp [60], resulting in a representation for the integrand in terms of classical polylogarithms up to weight four. More specifically, we are left with integrals of the form with n + m ≤ 4 and where P i is a polynomial and R i,k are rational functions. Note that, while individual terms in the sum are singular for λ → 0, 1, the sum is finite by construction, and so the integral is well defined. In order to perform the integration over λ, we rewrite the classical polylogarithms in terms of multiple polylogarithms [16] of the form G (a 1 (z), . . . , a n (z); λ), where a i (z) are rational functions of z using symbols [17][18][19][20][21] and the Hopf algebra structure of multiple polylogarithms [9,[22][23][24]. All the integrals can then easily be performed using the recursive definition of the multiple polylogarithms, G(a 2 , . . . , a n ; λ) = G(a 1 , . . . , a n ; 1) . (4.57) Finally, we observe that the results of the integration can also be expressed in terms of harmonic polylogarithms [36] of weight up to five, and we checked that the results are in agreement with the differential equation approach.

Conclusions
In this article, we have taken one more step towards the evaluation of the Higgs boson cross-section at N 3 LO in perturbative QCD. Specifically, we have analytically evaluated the one-loop contributions to the partonic cross-sections from 2 → 2 processes. This has required us to perform single real emission phase space integrals over the square of oneloop amplitudes. We performed these integrations using various independent methods. The motivation for this was, besides having a number of cross checks for the correctness of our results, to try out the most promising state of the art techniques for the evaluation of multiloop integrals. The (RV) 2 -corrections are simpler than other, yet unknown, contributions to the N 3 LO Higgs boson cross-section. However, they share many of the complexities that appear in corrections with mixed real and virtual radiation. Therefore, our calculations in this paper serve as a perfect testing ground for our computational techniques. In this publication, we have achieved two goals. We have obtained analytic expressions for the (RV) 2 partonic cross-sections which are valid for arbitrary values of the Higgs mass and energy as an expansion in ǫ through the finite part. Given that a calculation of the full hadronic Higgs boson cross-section at N 3 LO will most likely be achieved first as an expansion around threshold, we have obtained here a threshold series expansion for the (RV) 2 partonic cross-sections. The coefficients of the threshold expansion are valid for arbitrary values of ǫ. We performed our calculations with various methods and techniques.
Our first method is based on reverse unitarity and integration by part identities in order to reduce the integrals of the partonic cross-sections to 19 master integrals. For these master integrals we derived a system of differential equations which we could solve to all orders in the dimensional regulator in terms of iterative functions of the kind of hypergeometric functions which we defined explicitly for this purpose. Their series representations yield a threshold expansion for the master integrals and the partonic cross-sections. The newly introduced hypergeometric functions share many properties with multiple polylogarithms, for example they satisfy a shuffle algebra. They also appear to cover a wide range of known multi-dimensional hypergeometric functions. It will be an exciting new direction for further research to establish further properties of these functions and their usefulness in higher order perturbative corrections. We have also solved the differential equations for the master integrals order by order in ǫ in terms of harmonic polylogarithms. This yields the main result of this publication, which is the (RV) 2 partonic cross-sections for arbitrary values of the Higgs mass and the partonic center of mass energy as an expansion in ǫ through O ǫ 0 .
We have also followed a direct integration approach to obtain the same results. For this task, we introduced counter-terms at the level of the squared matrix element in order to subtract its collinear and soft divergences. Having rendered the phase-space integrand finite, we then expand in the dimensional regulator and perform a direct integration in terms of harmonic polylogarithms. The integration is made possible by embedding the classical polylogarithms that result form the expansion of the hypergeometric function in a larger space of multiple polylogarithms and exploiting their Hopf algebra structure.
We have also followed a different strategy, with a more restricted objective to obtain a threshold expansion of the cross-sections. Experience from NNLO has shown that very good approximations to the inclusive Higgs boson cross-section can be obtained with the first few terms of this expansion. We obtain the series around threshold by expanding the integrand of the loop amplitudes in soft, collinear and hard regions and perform the phase-space integration term by term. We have successfully applied this method to compute a sufficiently large number of terms in the threshold expansion to all orders in the dimensional regulator. Unlike some of the other methods we have used, the strategy of regions is generic, which makes it a particularly attractive option for computing the (RV) and (RV) 2 contributions in more complicated processes, where a direct evaluation in terms of polylogarithms may be unfeasible or just extremely difficult.
Our results are rather lengthy and we have only typeset parts of them in this document. Instead, we provide in the source of the arXiv submission two text files results.txt and masters.txt with the expressions for the partonic cross-sections and the master integrals respectively. The two files can also be downloaded from [42].
We believe that with our calculation and the methods that we have developed in this paper to be closer in our objective to compute the Higgs boson cross-section at N 3 LO. We look forward to applying the techniques presented here towards this objective.

A. Hypergeometric functions through iterated integrals
In this appendix we define a class of iterated integrals as also introduced in Section 4.2. First let us define the integral where we have abbreviated for later convenience a = a b . We have made use of Gauss' hypergeometric function with the third argument being the first argument increased by one. Next, we define recursively the n th iterated integral by F an,..., a 1 (x n , . . . , x 1 ; δ) = δ 0 dt t an−1 (1 − x n t) −bn F a n−1 ,..., a 1 (x n−1 , . . . , x 1 ; t) . (A. 2) The integration kernel t a−1 (1−ct) −b has the same form for every iteration step with indices a, b and argument c changing. Next, we derive a hypergeometric series representation for these iterated integrals. To simplify the expressions we rewrite eq. (A.1) and introduce a function f that is implicitly given by In the next step we integrate over the integration kernel and the F a (c; t) Using the identity (a + n) m = (a) n+m Γ(a) Γ(a + n) , we can write F a 2 , a 1 (c 2 , c 1 ; δ) = δ a 1 +a 2 (a 1 + a 2 )a 1 ∞ n,m=0 (a 1 + a 2 ) m+n (a 1 + a 2 + 1) m+n We now proceed iteratively, and find the following series representation for the iterated integrals F an, a n−1 ,..., a 1 (c n , . . . , with the abbreviations Following the the same procedure as for the sum representation we can derive a general Mellin barnes representation for our iterated integrals by utilizing the Mellin-Barnes representation of the Gauss Hypergeometric function This leads to F an, a n−1 ,..., a 1 (c n , . . . , These iterated integrals interpolate between multiple polylogarithms [16] and hypergeometric functions. In the framework of the above definitions the multiple polylogarithm is given by Here the indices of the iterated integrals only take the form 0 = 0 0 and 1 = 1 1 . Even for general indices we discover further similarities of these iterated integrals with multiple polylogarithms. As in the case of polylogarithms this class of hypergeometric functions may be written as multiple nested sums F an, a n−1 ,..., a 1 (c n , . . . , where k 0 = 0. The representation (A.10) may be useful to expand the iterated integrals in terms of the dimensional regulator (see, e.g., ref. [61]). The definition of these function as iterated integrals implies that they form a shuffle algebra, where Σ(i, n−i) denotes the set of all shuffles of n elements, i.e., the subset of the symmetric group S n defined by Σ(i, n − i) = {σ ∈ S n | σ −1 (1) < · · · < σ −1 (i) and σ −1 (i + 1) < · · · < σ −1 (n)}. (A. 12) To illustrate an application of the shuffle-product for generalized iterated integrals, let us look at the following example. We would like to integrate an iterated integral over a non-standard integration kernel. c 1 ; t). (A.13) To simplify the integral we make use of .14) and find Next, we apply the shuffle product and find Further identities among iterated integrals can be derived using integration-by-parts or by partial fractioning products of integration kernels. Further properties and parallels of generalized iterated integrals and generalized poly-logarithms are under investigation.

B. NNLO RV Master Integrals as hypergeometric functions
In this appendix we demonstrate how certain differential equations for master integrals appearing in physical cross-sections can be solved using iterated integrals as introduced in Appendix A. We consider the example of the master integrals contributing to the RV Higgs boson cross-section at N N LO. The master integrals were introduced and evaluated as an expansion in the dimensional regulator ǫ in ref. [12] and evaluated to even higher order in ref. [7]. Here we solve them to all orders in ǫ in terms of hypergeometric functions. The master integrals and the corresponding differential equations are given by We have abbreviated ∂ δ = ∂ ∂δ . Solid lines represent scalar propagators. The phase-space integral is represented by the dashed line and the cut-propagators are the lines cut by the dashed line. The cut propagator of the Higgs boson is depicted by the double-line. The complex conjugated one-loop integral is on the right-hand side of the phase-space cut. The differential equations were obtained using the methods described in Section 4.2.
The system of differential equations is decoupled and can be solved as described in Section 4.2. To solve the differential equations we require the following boundary conditions, which can be obtained from ref. [7], and all other cases vanish. These boundary conditions are given by the leading (soft) term of three master integrals in the limit of δ → 0. They are obtained using the methods described in Section 4.1. For convenience we will from now on set s = 1.
The first two differential equations are homogeneous and can be easily solved to give We find that the homogeneous solution to the differential equation of master Y 3 is vanishing and the inhomogeneous solution is according to eq. (4.46) given by To obtain this result we made use of eq. (A.1). As we proceed to solve the remaining two master integrals we find that the inhomogeneous solution to their differential equation is in turn dependent on Y 3 . We are able to find solutions to the inhomogeneous equation making use of the definition of our iterated integrals in eq. (A.2).
The iterated integrals with two indices contributing to Y 4 and Y 6 can be written as (1, 1, δ) = δ a 1 +a 2 (a 1 + a 2 )a 1 F 1,2 1,1 a 2 + a 1 a 2 + a 1 + 1 where we introduced the Kampé de Fériet function, Note that the generalized Kampé de Fériet function S 1 of Section 4.1 is a special case of eq. (B.21),

D. Results: Threshold Expansion
In this section, we present the results of our calculation for the functions η a;j X of Section 2.