NNLO phase space master integrals for two-to-one inclusive cross sections in dimensional regularization

We evaluate all phase space master integrals which are required for the total cross section of generic 2 ->1 processes at NNLO as a series expansion in the dimensional regulator epsilon. Away from the limit of threshold production, our expansion includes one order higher than what has been available in the literature. At threshold, we provide expressions which are valid to all orders in terms of Gamma functions and hypergeometric functions. These results are a necessary ingredient for the renormalization and mass factorization of singularities in 2 ->1 inclusive cross sections at NNNLO in QCD.


Introduction
Signals of novel physics, such as the recently discovered Higgs boson by the ATLAS [1] and CMS [2] collaborations, as well as standard candles, such as Drell-Yan production of electroweak gauge bosons [3,4], have been at the epicenter of experiments at hadron colliders. The analysis and interpretation of data relies crucially on high precision theoretical estimates, based on perturbation theory, for the rates of 2 → 1 processes. Most importantly, higher order corrections due to Quantum Chromodynamics (QCD) are generally large at LHC energies, altering the central value and lowering the theoretical uncertainty of many observables.
The state of the art in fixed-order perturbative QCD for 2 → 1 processes is next-tonext-to-leading-order (NNLO). The first inclusive cross section computation at this order was performed more than two decades ago for Drell-Yan production [5]. This calculation was followed by several computations for inclusive Higgs production cross sections [6,7,8,9,10,11,12,13,14,15] with the emergence of powerful computational techniques [6,7].
The precision of NNLO cross sections is quite high as indicated by the study scale variations. For example, the uncertainty due to scale variations in Higgs production via gluon fusion is less than 10% and in Drell-Yan production less than 2% (see, for example, ref. [16] and ref. [17]). However, a more reliable estimate of the uncertainty due to the truncation of perturbation theory would be the evaluation of yet another term in the series expansion. Now that a Higgs boson has been discovered, measuring its properties and couplings precisely becomes a priority. With this prospect of precision physics in the Higgs sector and the paramount importance of Higgs coupling measurements for our understanding of physics at high energies, it is evident that many important observables will need to be validated at an even higher order in perturbation theory, i.e. at next-tonext-to-next-to-leading order (NNNLO or N 3 LO).
In this paper, we present one of the ingredients needed for 2 → 1 inclusive cross sections at N 3 LO. In particular, we compute the double-real and real-virtual master integrals for a generic 2 → 1 scattering amplitude at one order higher in the dimensional regularization parameter ǫ = (4 − D)/2, with D the space-time dimension, than previously known [7]. The paper is organized as follows. In section 2, we motivate our computation by explaining where they enter the calculation of a N 3 LO observable. In section 3 we provide some details about the method we use. Sections 4 and 5 contain the results and we conclude in section 6.

NNLO contributions to N 3 LO observables
In this paper we are concerned with the computation of the fully inclusive cross section for the production of a massive state X. In perturbation theory this cross section can schematically be written as with S being the total center of mass energy, s the partonic center of mass energy carried by the momenta of the incoming partons, p 1 and p 2 , and m X the (on-shell) mass of particle X. The bare strong coupling constant is denoted as α are the bare parton distribution functions (PDF's) andσ ij,n is the partonic cross section for X + n jets, which in turn admits a perturbative expansion in the number of loops.σ ij,n may be written as a phase space integral over the squared X + n jets amplitude as follows, where we define the measure of the phase space volume for a massive particle X of mass m X and n − 1 massless particles by where (2.5) Using this notation, the Lorentz invariants appearing in the phase space integrals are defined through While the inclusive cross section σ X must obviously be finite order by order in perturbation theory, the individual pieces contributing to a given loop order are divergent. Final state infrared (IR) divergences cancel mutually between the real and virtual corrections, whereas ultraviolet (UV) and initial-state IR divergences have to be dealt with by replacing the bare coupling and PDF's by their renormalized counterparts, which requires the introduction of explicit counterterms proportional to poles in the dimensional regularization parameter ǫ multiplying lower order coefficients of the cross sections. As an example, ifσ (ℓ) denotes the ℓ-th order correction to the cross section, the N 3 LO UV/PDF counterterms can schematically be written as where the × may indicate a convolution in the case of the PDF counterterms. As a consequence, the poles of the counterterms produce finite contributions to the N 3 LO cross section from the higher orders in the ǫ expansion ofσ (ℓ) . The aim of this paper is to provide the NNLO master integrals which are required inσ (2) to one order higher in the ǫ expansion. At this order in perturbation theory three different contributions need to be taken into account. As an example, the NNLO correction to the partonic cross section for gg → H gets contributions from the following three types of interference diagrams, • double-virtual: the two-loop production amplitude interfered with the Born amplitude as well as the square of the one-loop amplitude for gg → H, e.g.,

⊗ ⊗
• real-virtual: the one-loop production amplitudes for gg → Hg, gq → Hq and gq → Hq interfered with the corresponding Born amplitudes, e.g.,

⊗ ⊗
First, the double-virtual corrections decouple from the phase space integration, reducing the computation of the two-loop phase space integrals effectively to the computation of the two-loop QCD form factor. The master integrals for the two-loop QCD form factor have been evaluated in ref. [18,19,20] to all orders in ǫ in terms of Γ functions and hypergeometric functions, which can easily be expanded in ǫ using the HypExp package [21]. For this reason we will not consider these integrals any further, but we simply list them in appendix A.1 for completeness. Next, in ref. [7] it was shown how to compute the phase space integrals describing the real-virtual and double-real corrections as a Laurent expansion in the dimensional regulator up to order ǫ 0 . More precisely, it was shown that all the relevant integrals can be reduced to a small set of master integrals, which were computed to the order in ǫ required for NNLO computations. In the rest of this paper we compute all the master integrals to one order higher in the ǫ expansion, thus preparing the ground for using them in computations beyond NNLO.

Phase space integrals from differential equations
In this section we give a brief account on how to compute the master integrals for the realvirtual and double-real phase space integrals. Our computations are based on the method introduced in ref. [7], which we shortly review in the following. In ref. [7] all real-virtual and double-real topologies were expressed as a two-loop forward scattering amplitude with, respectively, two and three cut propagators. Such cut integrals are related to the desired phase space integrals via Cutkosky's rule, Eq. (3.1) allows us to interpret a phase space integral as a loop integral, and thus we can apply results developed for the computation of Feynman integrals. In particular, these two-loop integrals can be reduced to a set of master integrals by virtue of the Laporta algorithm [22], using integration-by-parts (IBP) and Lorentz invariance identities, as implemented in the program AIR [23], with the additional constraint that integrals with some of the cut propagators not present vanish and need not to be considered in the reduction process [7]. The master integrals themselves are calculated using the method of differential equations [24,25]: we differentiate each master integral under the integration sign with respect to a kinematic invariant, in our case the mass m 2 X or equivalently the variable z. During the differentiation process, loop integrals with higher powers of propagators are produced. They can again be expressed in terms of the master integrals themselves via the IBP identities and the system of differential equations closes upon itself. If we denote the master integrals by F i (z, ǫ), we arrive at a set of first order ordinary linear differential equations for the master integrals, where the c ij (z, ǫ) are rational functions of z and ǫ which have poles in z at most at z = 0 and/or z = ±1. In the case where the set of differential equations is triangular order by order in ǫ, the set of equations may be solved order by order using standard techniques for ordinary differential equations. In particular, if ω(z) is a homogeneous solution to the linear differential equation then the solution to the inhomogeneous differential equation is given by the integral where α ∈ R is an integration constant. The poles appearing inside the integral in eq. (3.4) are determined by the poles in the coefficient A(z). Thus in our case the solution for the master integrals can be expressed through (iterated) integrals with poles at most at z = 0 and z = ±1. Integrals of this form lead naturally to harmonic polylogarithms [26], defined recursively through the iterated integrals H(a 1 , . . . , a n ; z) = z 0 dt f (a 1 , t) H(a 2 , . . . , a n ; t) , where a i ∈ {−1, 0, 1} and In the case where all the a i are zero, the integral (3.5) is divergent, and we define instead In this way we can express all the master integrals, up to the desired order in ǫ, in terms of harmonic polylogarithms up to weight four, which can be evaluated numerically in a fast and accurate way [27,28,29,30,31]. In order to fully determine the solutions of the differential equations, we need to fix the integration constants. This can be achieved by requiring the master integrals to take particular values, computed separately by other means, for some special value of z. Since 0 < z < 1, there are two natural choices for such special values. The case z = 0 physically corresponds to the situation where the produced particle X is massless. However, the limit z → 0 is in general not smooth, as the amplitude might develop new infrared poles in this limit. For this reason we choose the point z = 1 as an initial condition to the differential equation, which corresponds to the soft limit in which the momenta of all the final-state partons vanish. As this limit is important not only as an initial condition for the differential equations, we review the structure of the phase space integrals in the soft limit in the next section.

Soft limits of master integrals
In the previous section we argued that the point z = 1 is a good initial condition for the differential equations satisfied by the master integrals. There is however another reason why the point z = 1 deserves special attention. Indeed, some of the master integrals contain simple poles at z = 1, and without loss of generality we can always write such integrals in the form where n is an integer and the function F n (z, ǫ) is finite at z = 1. We can isolate the divergence associated to the soft singularity by expanding and the + indicates the common plus-prescription, It then becomes apparent that, if we want to know F(z, ǫ) up to a certain order in the ǫ expansion, we need to know F n (1, ǫ) at one order higher in the expansion in ǫ than F n (z, ǫ), since F n (1, ǫ) is multiplied by the pole in eq. (3.9). More generally, even if a given master integral has no pole as z approaches 1, when appearing inside an amplitude it might get multiplied by a coefficient such that the product develops at most a simple pole at z = 1. Hence, we have to use eq. (3.9) to perform the ǫ expansion and therefore, by the same logic, we always need to know the values of the master integrals at z = 1 to one order higher in the ǫ expansion. Thus, there is an independent physics motivation to study the soft limits of the master integrals in some detail. We define the soft limit of a master integral F(z, ǫ) as the (unique) function We stress that the right-hand side of eq. (3.13) is valid to all orders in ǫ. In addition we define the hard part F H (z, ǫ) of a master integral by (3.14) In sections 4 and 5 we will determine the soft limits of the master integrals directly to all orders in ǫ, without the help of the differential equations. Our knowledge of the soft limits combined with eq. (3.13) as an initial condition for the differential equations then allow us, at least in principle, to compute the analytic expressions for the master integrals up to any order in ǫ. In the rest of this paper we perform this task explicitly for all master integrals up to transcendental weight 4.

Definitions and conventions
In this section we compute the analytic expressions for the master integrals for the doublereal emission contributions. In ref. [7] it was shown that all phase space integrals for the production of a massive particle and two massless particles can be reduced to a linear combination of 18 master integrals, which we choose as follows: The solid lines denote the massive particle and dotted lines represent numerator factors. In all integrals we have pulled out an overall dimensionful scaling factor as well as an overall normalization P(ǫ) = 1 2 Note that the integral X 1 is just the phase space volume, which can be written to all orders in ǫ as For the other integrals no similar all-orders formulae are known to us. We therefore compute these integrals using the differential equation technique. It turns out that the system of differential equations for the double-real master integrals are triangular order by order in ǫ, except for X 1 and X 11 , which can be decoupled order by order by replacing X 11 by the linear combination After this change of variables the system is triangular order by order in ǫ, and we can immediately solve for the master integrals, provided that we know the soft limits of the master integrals which serve as an initial condition. The computation of the soft limits will be the topic of the rest of this section.

Soft limits of double-real master integrals
To evaluate the soft limits of the double-real master integrals we use the following exact parametrization for the 2 → 3 phase space, In this parametrization the propagators of the massless partons read where ξ =z (y 1 y 4 + y 2 y 3 + 2 cos(x 2 π) √ y 1 y 4 y 2 y) , (4.24) and This parametrization may be derived from the "energies and angles parametrization" of ref. [32] via the following transformation We now take the soft limit, z → 1. It is easy to see that in the soft limit we have the relations If we now perform the change of variables It is then easy to see from eqs. (4.23) and (4.29) that in the soft limit all the invariants, except for s 34 , take a fully factorized form. In particular the terms s 13 + s 14 and s 23 + s 24 effectively reduce to t 1 andt 1 respectively. This construction therefore allows one to derive all the soft limits of almost all master integrals, with the exception of X 17 and X 18 , in terms of simple Beta-functions, The integral X 17 may be trivialized in a similar way in the "hierarchical parametrization" of ref. [32]. Unfortunately, we are not aware of any parametrization which allows one to factorize all denominators of X 18 simultaneously. We proceed by expressing the integral in terms of energies and angles, i.e., by writing p i = E i (1, n i ), where n i is a unit vector in the direction of the spatial component of p i . The integral over the energies of p 3 and p 4 yields just another Beta-function, and we arrive at where dΩ (d−1) i is differential volume element of the d − 1 dimensional solid angle of p i and the angles are defined as cos θ ij = n i . n j . We first use the result [33] The integral then becomes , eq. (4.33) may be written as where we used the recursive definition of the p F q function. Hence we arrive at expressions valid to all orders in ǫ for the soft limits of all the double-real master integrals. We observe that in all cases, except for X S 18 , the results are proportional to the soft limit of the phase space volume, the constant of proportionality being a rational function of z and ǫ. We therefore define In this normalization the results for the soft limits of the master integrals read (4.48) One might wonder whether our result for S 18 , eq. (4.48), can be expressed to all orders in ǫ through ratios of Γ functions only. Without loss of generality we can assume that in that case the 3 F 2 function entering S 18 could be written in the form, where R(ǫ) is a rational function. In the following we argue that the form (4.49) can be excluded for S 18 , i.e., that 3 F 2 (1, 1, −ǫ; 1 − ǫ, 1 − 2ǫ; 1) cannot be reduced to a ratio of Γ functions only. Indeed, if we assume that the 3 F 2 function in eq. (4.48) can be written in the form (4.49), then, because of the well-known relation its ǫ expansion will to all orders in ǫ only involve zeta values of depth one, i.e., it will be possible to write its ǫ expansion without multiple zeta values (MZVs). Multiple zeta values are known to satisfy many identities among themselves, and it is conjectured that all the identities among MZVs are generated by the double shuffle relations. Using the double-shuffle relations, it can be shown that the first time a MZV cannot be written as a polynomial of MZVs of depth one is at weight 8, and this irreducible MZV can be chosen to be ζ 5,3 . Using HypExp we can perform the expansion of 3 F 2 (1, 1, −ǫ; 1 − ǫ, 1 − 2ǫ; 1) up to O(ǫ 8 ). The result obtained by HypExp is a complicated combination of MZVs of weight 8.
Using the PSLQ algorithm we fit this combination of MZVs to a basis of MZVs of weight 8, which we choose to be {ζ 4 2 , ζ 2 ζ 2 3 , ζ 3 ζ 5 , ζ 5,3 }. We obtain We see that at O(ǫ 8 ) the expansion of 3 F 2 (1, 1, −ǫ; 1 − ǫ, 1 − 2ǫ; 1) involves an irreducible MZV of weight 8. Conversely, we checked that if we remove this MZV from our basis at weight 8, the PSLQ algorithm fails to converge. We thus conclude that the 3 F 2 function appearing in S 18 cannot be written as a ratio of Γ functions of the form (4.49).
Having at our disposal both the differential equations and the soft limits for all the master integrals, we can in principle solve the differential equations to any order in ǫ.
We have explicitly performed this task and computed the ǫ expansion of all the master integrals up to transcendental weight 4. As the results are rather lengthy, we do not show them here explicitly but collect them in appendix A.3. The results only involve harmonic polylogarithms up to weight 4. In ref. [34] it was shown that up to weight 4 almost all harmonic polylogarithms can be expressed through classical polylogarithms only, and that only three new functions are needed. In our case only one new function appears, which we choose to correspond to the harmonic polylogarithm Li 2,2 (−1, z) = −H(0, 1, 0, −1; z) . (4.52) We have checked that our results agree with the results of ref. [7] up to terms in the ǫ expansion with transcendental weight at most three. In addition we have checked all our results numerically for a few nontrivial values of z. For this purpose we used the method of non-linear mappings introduced in ref. [32].

Master integral definitions
In ref. [7] six master integrals for the real-virtual contributions to a 2 → 1 partonic cross section were identified. Defining the prefactor as well as the one-loop integrals the real-virtual master integrals may be defined as the following two-body phase space integrals,

Soft limits of real-virtual master integrals
For the evaluation of the soft limits of the real-virtual master integrals we employ the strategy presented in ref. [35]. This strategy involves the following steps: 1. Substitute known expressions for the one-loop master integrals which are valid to all orders in ǫ.

Apply analytic continuation formulae to the hypergeometric functions which express the box integrals, if required.
3. Evaluate the integrals at z = 1.
We find that the soft limits of all master integrals can then be trivially evaluated as Betafunctions: Note that only the real part of the master integrals is shown, as the imaginary part never enters the computation of a physical observable. As in the double-real case, we list the results for generic z that we obtained by solving the differential equations in appendix A.2. The expressions agree with the results of ref. [7] up to terms in the ǫ expansion with transcendental weight at most three. The new orders were checked numerically using the approach of [35]. This strategy is as follows: 1. Substitute known expressions for the one-loop master integrals which are valid to all-orders in ǫ.

Apply analytic continuation formulae to the hypergeometric functions which express the box integrals, if required.
3. Expand the hypergeometric in terms of polylogarithmic functions.
4. Expand the real emission singularities (these are factorized here) in terms of deltaand plus-distributions.

Conclusions
We have given results for all double-real and real-virtual master integrals needed in the calculation of a generic 2 → 1 inclusive scattering cross section. The analytic results for the master integrals are provided in Maple input form as ancillary material to the arXiv submission. The previously known expressions have been supplemented with one more order in the ǫ expansion. Together with the previously published double-virtual results [20] listed in appendix A.1, all contributions to a N 3 LO cross section involving NNLO master integrals can thus be computed, as outlined in section 2.
We have divided the master integrals into a soft and a hard part, each of which is given separately. This separation is motivated in a twofold way. First, the soft limits provide a natural boundary condition to solve the differential equations for the master integrals. Second, divergences due to soft radiation have to be subtracted, leaving poles at z = 1 and thus requiring the knowledge of the master integral to one order higher in the ǫ expansion, as explained in section 3.1.
We have found the method of differential equations to be very adequate for this type of task, and see no obvious problems in applying the same techniques to integrals of higher complexity. Such integrals could include the triple-real master integrals which constitute another part of a 2 → 1 N 3 LO cross section. In that case, the main difficulty might be the vastly increased amount of master integrals, as well as a much more intertwined set of differential equations due to more complicated IBP identities for the triple-real topologies. Once these difficulties overcome, one should be able to perform analytically all the phase space integrals appearing in 2 → 1 cross sections at N 3 LO, and thus, by combining them with the recently computed three-loop QCD form factor [36,37], arrive at the first prediction for an LHC observable computed perturbatively up to N 3 LO.

A. Results for the master integrals
In this appendix we collect the analytic results for the hard part of the real-virtual and double-real master integrals omitted throughout the main text. As the imaginary part of the master integrals does not enter the computation of any physical observable, we only show the real part. For completeness we also include the analytic expressions for the twoloop virtual master integrals, which can easily be obtained from the expressions for the two-loop form factor master integrals of ref. [18,19,20].