On the calculation of soft phase space integral

The recent discovery of the Higgs boson at the LHC attracts much attention to the precise calculation of its production cross section in quantum chromodynamics. In this work, we discuss the calculation of soft triple-emission phase space integral, which is an essential ingredient in the recently calculated soft-virtual corrections to Higgs boson production at next-to-next-to-next-to-leading order. The main techniques used this calculation are method of differential equation for Feynman integral, and integration of harmonic polylogarithms.


Introduction
To maximize the physics outcome of the LHC program, theoretical predictions for signal and background processes have to match the unprecedented accuracy of experimental measurements. For QCD initiated processes, a standard approach to improve the accuracy of theoretical prediction is by calculating perturbative radiative corrections beyond Leading Order (LO) approximation. At Next-to-Leading Order (NLO), such an calculation consists of virtual correction, which involves virtual parton circulated in the loop, and real correction, which involves real on-shell parton emission. When considering corrections beyond NLO, mixed real-virtual correction is also required. In the past few decades, significant efforts have been made to develop analytic techniques for loop amplitude calculation. In particular, a variety of methods have been developed for the analytic calculation of loop integral, cf. [1] and references therein. For phase space integral, techniques for analytic calculation are considerably less developed. An obvious reason is that experimental measurements often involve complicated phase space cuts, which are very difficult to implement analytically. The standard approach is to perform the phase space integral numerically, keeping in mind that appropriate subtraction terms have to be added in order to render the integral finite [2,3]. In that case analytic inclusive phase space integrals can be used in the construction of infrared subtraction terms for Next-to-Next-to-Leading Order (NNLO) calculation [4]. Moreover, when appropriate subtraction method was not yet available, analytic integral was the only method for obtaining cross section or distribution at NNLO at hadron collider [5,6,7,8,9].
Recently, an important step towards the calculation of Higgs production at Next-to-Nextto-Next-to-Leading Order (N 3 LO) has been taken in Refs. [10,11], in which the soft-virtual corrections are obtained. The soft-virtual corrections are cross section in the threshold limit, where radiation in the final state is constrained to have low energy. Since the radiation energy is low, calculation can be done in the eikonal limit, where QCD partons are emitted from light-like Wilson lines which parameterize the directions of incoming partons. Calculation of soft-virtual corrections at N 3 LO requires the calculation of pure virtual corrections through to three loops, and soft phase space integrals with one, two, and three partons in the final state. For soft phase space integrals, the ones with three partons in the final state are the most challenging one. Such triple soft phase space integrals are first calculated in the pioneering work [12], using various powerful techiques for Feynman integrals, ranging from Melin-Barnes transformation to symbols and coproduct.
In Ref. [11], a completely different approach is taken. Auxiliary integrals are introduced in the intermediate steps, as well as nontrivial dependence on internal auxiliary scales. The auxiliary integrals are then solved by the method of differential equation for Feynman integral. Finally the auxiliary scales are integrated over and the desired triple soft phase space integrals are obtained. This approach makes heavy use of the property of scale invariance for eikonal integral. It is worth mentioning that the triple soft phase space integrals are single scale integrals with trivial dependence on the single scale, therefore prohibit direct application of the method of differential equation. By introducing the auxiliary integrals, additional nontrivial scales are introduced into the problem, and nontrivial system of differential equation can be derived. While differential equation approach to single scale integrals have been employed before in different context [13,14] by introducing external scales, we find that auxiliary scales being introduced internally in the current calculation is particularly interesting, and significantly simplifies the calculation at all stages. In this paper we document the detail of such a calculation, in the hope that it will be useful in other problems. This paper is organized as follows. In Sec. 2, we explain the method in detail, using the NNLO double soft emission phase space integral as example. In Sec. 3, we apply this method to triple soft phase space integrals. In Sec. 4 we give a brief conclusion and outlook. Detail for solving the system of differential equation for triple soft integrals is given in Appendix. A.
2 Differential equation approach to single scale soft phase space integral

Introducing the method
The method of differential equation [15,16,17] is powerful in the calculation of Feynman integrals with multiple scales. For single scale problem, the dependence on the scale is trivial from dimensional analysis. The nontrivial part of the integral is a function of dimensional regularization parameter only 1 . Thus the method of differential equation can not be applied directly. In some cases, one can introduce external auxiliary scales such that nontrivial differential equation for Feynman integral can be derived. The integrals with additional external auxiliary scales can then be related to the original integrals by taking limit of the auxiliary external scales [13,14]. Unfortunately, the calculation of Feynman integral becomes more complicated as the number of scales increased.
In this work, we solve single scale soft phase space integral by introducing internal auxiliary scales. Differential equations with respect to these internal scales are derived and solved. Then the internal auxiliary scales are integrated over, and the desired single scale integral is recovered. There are several advantages in this approach. First, the internal auxiliary scales will be introduced by inserting a delta function. Some of the momentum integral can be integrated out trivially by the delta function. Effectively, the number of loops involved in the calculation is reduced by one. Reducing the number of loop significantly simplifies the complexity of Intergration-By-Parts (IBPs) reduction [18,19,20]. Second, the property of scale invariance for eikonal integral can be fully exploited after the introduction of internal auxiliary scales. The consequence is that auxiliary integral only depends on a scale invariant combination of the auxiliary scales. Therefore the functional dependence of auxiliary integral on the auxiliary scales has very restricted form. Third, the internal auxiliary scales are introduced in a way that integrating over them is also trivial. In fact, only simple interative integrals of harmonic polylogarithms (HPLs) [21] are needed for this calculation, in contrast to the extensive uses of multiple polylogarithms in Ref. [12].
The single scale soft phase space integrals discussed in this work are related to the soft-virtual corrections for color singlet particle production at hadron collider. Perhaps the most important example is the soft-virtual corrections for Higgs production [10,11]. Consider the process (2.1) Thanks to QCD factorization, the hadronic cross section for this process can be written schematically as a convolution of parton distribution functions and partonic cross section The partonic cross sectionσ(z) depends on partonic threshold variable z = M 2 H /ŝ 2 , where M H is the mass of Higgs boson, and √ŝ is partonic center-of-mass energy. The soft-virtual corrections are referred to the leading term ofσ(z) in the expansion of z around 1. Since the energy of final state radiation is , calculation of soft-virtual corrections is equivalent to calculation of differential distribution in E X , in the limit where E X is much smaller than all other scales in the problem. In this limit, the leading power contributions in (1 − z) come from the partonic process In QCD perturbation theory, X contains various combination of massless QCD parton. For example, at NNLO, the double-emission contribution contains the process while the triple-emission contribution at N 3 LO contains the process where p 2 i = k 2 i = 0. The phase space integrals in the limit E X ≪ M H are single scale integrals, whose dependence on E X can be determined uniquely by dimensional analysis. The main purpose of this work is the calculation of these integrals.
We describe in detail the method mentioned above, taking double-emission soft phase space integrals as an example. After application of IBPs reduction 3 , the result for double-emission contribution to NNLO soft-virtual corrections can be expressed in terms of two nonzero two-loop master integrals,

6)
2 It also depends on Higgs mass M H , top-quark mass M t , renormalization scale µ R and factorization scale µ F , which are irrelevant to our discussion. 3 We use both REDUZE 2 [22] and LiteRed [23] for IBP reduction.
where γ E = 0.577216 . . .. The integrals are defined in d = 4 − 2 ǫ dimension to regulate infrared singularities. The dependences of I 1 and I 2 on the kinematical variables E X , M H , p 1 and p 2 are highly constrained. From dimensional analysis, we have where we have also made use of the fact that the scale transformation for the soft integrals is not anomalous, i.e., anomalous term of the form (p 1 ·p 2 ) aǫ is always absent.
In the rest of this paper we will set 2E X M H = 1 for simplicity because their dependence are known exactly. We will also let p 1 ·p 2 = 2, and introduce the standard notation As explained above, I 1 and I 2 are single scale integrals. Therefore the method of differential equation for Feynman integral can not be applied directly. The solution is to introduce internal scales by inserting an unit operator into the original integrals, where we have introduced the auxiliary integrals .
Integrals of the form in Eq. (2.12) are first introduced in the context of fully exclusive soft function in ref. [24], and have been used in the calculation of threshold soft function for direct photon production [25] and boosted top-quark pair production [26]. We will see later that they are particularly convenient for evaluating soft phase space integral. An important feature of the auxiliary integrals in Eq. (2.12) is that they transform in a simple manner under Eq. (2.8). The scale symmetry and dimensional analysis determine that where functions f i (y) depend on a single kinematical variable y = l 2 /(l + l − ). To determine the functional form of f i (y), one only needs to calculate the integrals in Eq. (2.12) with the constraints 14) The full functional dependence can easily be recovered using Eq. (2.13). For this reason, we will also denote the auxiliary integrals as This is an important simplification because the number of independent scale have been significantly reduced, but not reduced too much such that there is still nontrivial dependence on the scale. The auxiliary integrals in Eq. (2.12) have been calculated to all order in dimensional regularization parameter in Ref. [24].
Results in Eq. (2.16) can be obtained by parameterizing the phase space integral in light-cone coordinate and performing the phase space integral explicitly. However, due to the quadratic constraint imposed by the delta function δ((l − k 1 ) 2 ), the resulting integral is not easy to computed. Alternatively, they can be obtained by dispersion method [27], as along as the discontinuities of the integrals are well understood. While the dispersion approach is intrinsically interesting, we prefer not to elaborate it here and leave it for future work. Instead, we use the method of differential equation, explained in the next section, for the calculation of such integrals.

Solving the auxiliary integrals by differential equation
Method of differential equation for Feynman integral is originally proposed for loop integral only. It was then realized that the same method can also be applied to phase space integral with small modification [7,9]. The reason is that on-shell condition for phase space integral can be regarded as Feynman propagator for the purpose of IBP reduction or calculating derivative, Cut propagators are then recovered as delta function after IBP reduction. The only difference of phase space integral IBP reduction with conventional loop integral reduction is that integral with zero or positive power of p 2 should be set to zero, because they do not give rise to discontinuity.
The derivation of differential equation for phase space integral can then be proceeded similar to loop integral. Taking derivative w.r.t to x can be expressed as The expression on the right-hand side of Eq. (2.18) can be reduced using IBP identities. We obtain It is a straightforward exercise to solve the system order by order in ǫ, up to undetermined integration constants. The differential equation for I 1 (x) is homogeneous. Its integration constants have to be fixed by explicit calculation. Using the result given in Eq. (2.16) [24], I 1 (x) expanded in ǫ is given by where we have only expand through to O(ǫ 4 ) for simplicity. In Eq. (2.20), ζ n is Riemann's zeta value, H w (x) are HPLs, first introduced in physics in Ref. [21]. In our case, the required HPLs are defined recursively by If w consists of 0's only, then it is defined to be H 0n (u) = ln n u/n!. As a standard notation, we replace (k − 1) 0's followed by a 1 with k in the weight vector. For example, H 0,0,1 (u) ≡ H 3 (u). Substituting I 1 (x) into Eq. (2.19) and solving the differential equation for I 2 (x), we obtain 2 are integration constants yet to be specified. An obvious way to fix these constants are calculating them on the boundary explicitly, x = 0 and x = 1. However, I 2 (x) is singular in either of these two limits, as can be seen from the known result, Eq. (2.16). Setting naively x = 0 or x = 1 in the integrand doesn't work and one has to perform expansion carefully along the line of the method of expansion by region [28]. This however requires further study on the application of the method of expansion by region to phase space integral, which is beyond the scope of current work. Instead, we fix these constants by going to d = 6 − 2ǫ dimension and requiring that I 2 (x) in 6 − 2ǫ dimension vanishes in the limit of x → 0. Using dimensional recurrence relation [29] 4 , I 2 (x) in 6 − 2ǫ dimension is given by The condition lim x→0 I 6−2ǫ 2 is the (i + 1)-th element of the vector b 2 . We have therefore determined the nontrivial auxiliary integral I 2 (x) in almost an algebraic way. In particular, the integration constants are fixed by regularity condition alone.

From auxiliary integral to soft phase space integral
With the auxiliary integrals at hand, we are ready to solve the soft phase space integrals in Eq. (2.6) and (2.11), in which we are actually interested. To this end, we need to restore the full dependence of I i (l, p 1 , p 2 ) on l, l + and l − . Using dimensional analysis and scale symmetry, we have (2.27) Since neither the integrand nor the delta function constrain depends on l ⊥ , we can integrate the transverse angular components out explicitly where Introducing the parameterization

30)
l + and l − can be easily integrated out in closed form, because I i (y) are functions of u only, For I 1 the remaining integral in u can be readily carried out order by order in ǫ. However, additional complexity arises for I 2 , for which the explicit integral can be written as We note that the integrand contains an unregulated singularity at u = 0. To regulate the integral, one approach is to extract the ln u dependence of I 2 (u) to higher order in ǫ and resum it as u aǫ , and then apply the method of sector decomposition [30] for the integral. Alternatively, the singularity becomes a spurious one in 6 − 2ǫ dimension, and we can calculate it directly without any regularization procedure. We can then use dimensional recurrence relation for the soft phase space integral to obtain the result in 4 − 2ǫ dimension. The advantage of this approach is that there is no power divergence in the integrand in 6 − 2ǫ dimension, and the integral can be straightforwardly carried out order by order in ǫ. Similar strategy has been employed in the direct calculation of soft phase space integral [12], in the calculation of multi-box diagrams [31], and in searching for quasi finite master integral basis [32]. In d = 6 − 2ǫ dimension, Eq. (2.32) is now given by where we have kept the factor e ǫγ E unchanged under dimensional shift by convention. Eq. (2.33) still contains the singular factor 1/u in the integrand. However, since the function I 6−2ǫ 2 (u) vanishes in the limit of u → 0, the integral is finite. We can then expand the integral order by order in ǫ, and integrate in the variable u using the definition of HPLs in Eq. (2.21) and integration-by-parts relations. In this way we obtain results for the original soft phase space integral in d = 6 − 2ǫ dimension, 6 16 , where we have kept the expansion through to transcendental weight 6. Then we can convert the results back to d = 4 − 2ǫ dimension ones using the following dimension recurrence relations  6 16 , These are the correct results for the NNLO soft phase space integral. We have seen in this section that the combined use of differential equation, integral of HPLs, and dimensional recurrence relation leads to an almost algebraic way for the determination of soft phase space integral. We will apply this method to the calculation of soft triple-emission phase space integrals in next section.

Soft triple-emission phase space integral
For soft-virtual corrections at N 3 LO, we need to calculate phase space integrals with three soft partons emitted in the final state. After reduction by IBP identities, we find a set of 7 master integrals need to be calculated [11], , where we have used the short-hand notation The calculation of these integrals follows closely the method given in Sec. 2.

Auxiliary integrals for soft triple-emission phase space integrals
As in the double-emission case, we introduce internal auxiliary scale by inserting into Eq.
We can then write the auxiliary integrals as , Topo(1, 1, 1, 0, 0, 0, 1, 1 , J 4 (x) = Topo(1, 1, 1, 0, 1, 1, 0, 1, 1) , Topo(1, 1, 1, 1, 1, 0, 0, 1, 1 Topo(1, 1, 1, 1, 1, 1, 1, 0 , (3.5) where they correspond to the original integrals J 1−6 . The calculation of J 7 requires additional care and will be addressed later in Sec. 3.2. The calculation of the auxiliary integrals in Eq. (3.5) follows closely the procedure illustrated in Sec. 2. We present the results here and leave the details to appendix A, where we have given the results through to transcendental weight 4 only. The corresponding expression through to weight 6 in electric form can be requested from the author. Our next step is to integrate these auxiliary integrals over the auxiliary scales to obtain soft triple-emission phase space integral. As was explained in Sec. 2, to further integrate the auxiliary integrals, we need to go do d = 6 − 2ǫ dimension. Using dimension recurrence relations generated automatically by LiteRed, we obtain the auxiliary integrals in 6 − 2ǫ dimension. After restoring the full kinematic dependence of the auxiliary integrals using scaling symmetry and dimensional analysis, we can straightforwardly calculate the 6 − 2ǫ dimension soft phase space integrals using Eq. (2.27). Through to weight 6, they are given by Finally, using dimension recurrence relations for soft phase space integral generated by LiteRed, we obtain the results in 4 − 2ǫ dimension, It is interesting to note that the rational part of J 5 and J 6 equals up to an overall normalization.

The factorizable integral
We still have the last integral J 7 to calculate. While it is straightforward to solve the corresponding auxiliary integral using differential equation, it is however difficult to fix the integration constants. In particular, the regularity condition does not fix the integration constants. To calculate this integral, we observe that the set of momentum {k 1 , k 2 } and {k 3 } are entangled together only by the delta function, It is well-known in resummation study that such entanglement can be factorized into a doubleemission integral and a single-emission integral via Laplace/Fourier transformation. Even more straightforwardly, we can consider a one parameter integral . (3.10) The original integral is simply J 7 (1). To factorize the integral, we rewrite it as We note that by dimensional analysis, J 7,12 (ω) is proportional to I 2 in Eq. (2.36). In d = 4 − 2ǫ dimension, we have and J 7,3 (ω) is the simple soft single-emission phase space integral, (3.14) Substituting Eqs. (3.13) and (3.14) into Eq. (3.11), we obtain Eqs. (3.8) and (3.15) are the main results of this work. Similar integrals have also been computed in Ref. [12] using completely different method and our results agree perfectly with theirs after taking into account the difference in the normalization.

Conclusion and outlook
In this work we have presented the calculation of soft triple-emission phase space integrals, which are relevant for soft-virtual corrections to Higgs production at N 3 LO [10,11]. The major techniques used in this calculation are differential equation for auxiliary integrals, and integration of harmonic polylogarithms. The essential idea is to introduce internal scales into single scale integrals. Since integral in the eikonal limit enjoys the property of scale invariance, the internal scales should be introduced in a way that doesn't spoil the scale invariance. In this way the resulting auxiliary integrals obey strong constraints from scale invariance, and therefore simplify significantly the solution of differential equations. A particular nice feature of the system of differential equation is that, all the nontrivial integration constants can be determined by regularity condition in 6 − 2ǫ dimension. Furthermore, the solution of differential equations, as well as the subsequent integration of auxiliary scales, only involves the well-known generalization of logarithms and polylogarithms, the harmonic polylogarithms. Thanks to all these techniques, we are able to calculate the soft phase space integrals in an almost algebraic way. The meaning of such a calculation is two folds. First, it provides an alternative method to the one used in the calculation of Ref. [12]. Second, it provides the first independent confirmation of the integrals calculated in Ref. [12].
• If the master integrals in the x → 0 limit are free of logarithmic singularity, then they vanish in that limit.
These conditions can be inferred by studying the asymptotic expansion of corresponding auxiliary integrals, treating the cut propagators as normal propagators. For the integration constants of J 1 (x) and J 2 (x), we calculate directly using the dispersion method, to which we wish to come back in the future.