Five-Particle Phase-Space Integrals in QCD

We present analytical expressions for the 31 five-particle phase-space master integrals in massless QCD as an $\epsilon$-series with coefficients being multiple zeta values of weight up to 12. In addition, we provide a computer code for the Monte-Carlo integration in higher dimensions, based on the RAMBO algorithm, that has been used to numerically cross-check the obtained results in 4, 6, and 8 dimensions.


Introduction
Nowadays, perturbative calculations play the key role in describing data from high-energy particle colliders, such as the LHC, as well as in improving the precision of numerical parameters in the Standard Model and other models. It is clear now that higher-order calculations will play an even more crucial role in processing data from future high-luminosity colliders, like the FCC or the ILC, where theoretical errors will dominate over experimental statistical errors. These arguments motivate us to make one step forward beyond available fully-inclusive phase-space integrals for a four-particle decay [2] and calculate a set of yet unknown integrals that corresponds to a five-particle decay of a color-neutral off-shell particle in Quantum Chromodynamics.
A particular application of these integrals we have in mind is the extraction of NNLO time-like splitting functions [3] from a semi-inclusive one-particle decay process, as for example discussed in [4,5]: the integrals we will be looking at correspond to a fully inclusive cross section, and can be used to determine integration constants when calculating exclusive quantities using the method of differential equations.
In this article, we focus on the calculation of master integrals that can be used to express any other integral of the corresponding topology provided a set of integration-by-parts rules (IBP) [6] is known. Our approach is based on techniques for solving dimensional recurrence relations (DRR) [7] described in [8,9]. In particular, we use the DREAM package [8] to obtain numerical results for the desired integrals with 2000-digit precision, and restore their analytical form in terms of multiple zeta values (MZV) [10,11,12] up to weight 12 using the PSLQ method [13] as implemented in Mathematica. We also present a Monte-Carlo code, based on the RAMBO algorithm [14], for numerical integration of the phase-space integrals in arbitrary (integer) number of dimensions that has been used to check consistency of the obtained results.
This article is organized as following. In Section 2 we introduce our notation and describe our calculational method in more detail. In Section 3 we provide complete results for four-particle integrals and discuss numerical cross-checks using Monte-Carlo integration. In Section 4 we make our final remarks.
Accompanying to this paper are the auxiliary files on the arXiv 1 containing the complete master integrals with MZV weight up to 12, as well as the Monte-Carlo integration routines with the corresponding results.

The Method
We start by identifying a set of five-particle phase-space master integrals using two different approaches for consistency. As the first approach we exploit the equivalence of IBP rules for cut and ordinary propagators, and obtain the complete basis of phase-space master integrals by taking all five-particle cuts of the 28 master integrals for four-loop massless propagators found in [15], discarding those that do not correspond to the squared matrix elements of the 1 → 5 process (i.e. only leaving graphs which are bipartite), and reducing the remaining integrals with Laporta-style IBP reduction [16,17] as implemented in FIRE5 [18]. As an alternative approach, we construct the complete expression for the total cross section of the 1 → 5 process in QCD using QGRAF [19] 1 https://arXiv.org/abs/1803.09084   Table 1: Cut diagrams for five-particle phase-space master integrals in QCD. Dashed lines represent cut propagators and carry final-state momenta p 1 , . . . , p 5 . Labels represent propagators, so that "123" corresponds to p 1 + p 2 + p 3 and "012" to q − p 1 − p 2 (where q is the initial-state momentum, i.e., q = p 1 + · · · + p 5 ).
and FORM [20], and then reduce it with the help of FIRE5. Both methods give 31 master integrals listed in Table 1, with each having up to 6 unique propagators. Our notation for these integrals is j are propagators that take the form of invariant scalar products s kl...q = (p k + p l + · · · + p q ) 2 , (2.2) dPS 5 is a five-particle phase-space element in D dimensions and S Γ is a common normalization factor chosen for convenience 2 to be With this normalization and knowing the volume of the complete N-particle phase space 3 we can already fix the value of F 1 as: Next, with the help of LiteRed [21] and FIRE5 [18] we derive a set of lowering dimensional recurrence relations which express master integrals in D + 2 dimensions in terms of master integrals in D dimensions: In the general case M(D) is expected to have a block-triangular structure, but in our case it can be shuffled into triangular form. The general structure of our M(D) can be visualized as: In other words, each F i only depends on itself and master integrals from lower sectors: As an example, for F 28 we have: The general solution for system (2.9) can be written as: where is a homogeneous solution, which can be found just from the diagonal entries, M ii (D); is a partial solution, which only depends on the integrals from subsectors, and can be constructed analytically as an indefinite nested sum, or evaluated numerically with DREAM, as long as F 1 (D), the lowest integral, is known; • ω i (D) is an arbitrary periodic function, such that ω i (D + 2) = ω i (D), which cannot be determined from the DRR relations alone, and needs to be fixed separately.
In the general case this solution can be evaluated as a series in ε around any dimensionality using the DRA method [22], which requires the analysis of singularities and the asymptotic behaviour of F i (D) in the limit of D going to imaginary infinity. Our case contains two important simplifications, and will not require such analysis.
The first simplification is that since M(D) is triangular, the homogeneous part of eq. (2.9) decouples into a set of first order difference equations: For rational M ii (D) written in the following form (compare to eq. (2.10)): the homogeneous solution can immediately be found as: (2.14) Explicitly, for H 28 (D) we have: The second simplification, is that as we will argue further, all ω i (D) for i > 1 are zero. To see this, first let us look at the asymptotic behavior of F i (D) at large D. Rewriting eq. (2.1) as an integral over invariants s i j gives and Ω k is the surface area of a unit hypersphere in k-dimensional space If ∆ N (s i j ) has a unique global maximum inside the integration region, we can apply Laplace's method to eq. (2.16) and find its asymptotic as where C i is a constant that depends on the location of the maximum and the denominators D (i) j , but not on D.
The global maximum of ∆ N is reached when all s i j (i = j) are identical and equal to 2 N(N−1) . Geometrically this configuration corresponds to the vectors p i pointing to the vertices of a regular N-hedron embedded into Euclidean space of (N − 1) dimensions. The maximum value is then and explicitly we get It follows that all F i (D) have identical asymptotic behavior up to a constant C i . As a confirmation, it can be shown that eq. (2.21) is asymptotically the same expression as we had for F 1 in eq. (2.6). Now we can compare this asymptotic for F i (D) to the asymptotic for H i (D) from eq. (2.14). Indeed, for H 28 (D) from eq. (2.15) we have: This grows asymptotically exponentially faster than eq. (2.21). In fact, all H i (D) for i > 1 grow exponentially faster than F i (D), and since F and H are connected via eq. (2.11), this can only happen if the corresponding periodic functions ω i (D) are zero.
Thus, to find F i (D) we only need to find R i (D), the inhomogeneous solutions to eq. (2.9). We compute them as a series in ε = (4 − D)/2 using DREAM with 2000-digit accuracy, and then restore the analytical form of the series coefficients in terms of MZVs using the PSLQ method [13]. This way we obtain the analytical result for all master integrals up to MZVs of weight 12 using the corresponding bases from [10] and the SummerTime package [12] for their numerical evaluation. Corresponding expressions are presented in the auxiliary files on the arXiv.

Four-Particle Integrals
As the first consistency check of our method we reproduce results for four-particle phase-space integrals reported in [2]. We perform all the steps described in Section 2. Generating the IBP rules with the help of LiteRed and then proceeding with DREAM we obtain the final result with 2000-digit accuracy. The series reconstructed with PSLQ, and truncated to MZV weight 6 are (using the original notation, and omitting S Γ and q 2 factors): These values are in a complete agreement with the known results, and we have included the series up to MZV weight 12 among the auxiliary files on the arXiv.

Numerical Verification
As another cross-check we have calculated the leading terms of the ε-expansion of F i numerically using the direct way: through Monte-Carlo integration of eq. (2.1) over the phase-space. While such a technique cannot be easily applied to divergent integrals, we can sidestep this issue by noting that our master integrals only suffer from IR divergences that disappear already at D = 6. This way we can check several leading terms of the expansion at D = 4 − 2ε by calculating the corresponding integrals in D ≥ 6 since both are connected by dimensional recurrence relations.
To calculate a finite integral of the form eq. (2.1) we choose a uniform mapping from a hypercube into momentum coordinates using an algorithm similar to RAMBO [14] but extended into arbitrary (integer) D. Then we calculate the integrand from scalar products of the momenta, and finally we integrate over the hypercube using the Vegas [23] implementation from CUBA [24].
Note that although the integrals we are calculating are finite, the integrands are not. Exposing an integration algorithm like Vegas to such infinities may lead to unpredictable behaviour, so as a precaution we choose to regulate these infinities by adding a small parameter α to the denominator of the integrand, and then to calculate the integral with progressively smaller values of α (from 2 −30 to 2 −100 ), checking if convergence was reached afterwards.
The results of this method are summarized in Table 2, and show good agreement between numerical and analytic results. Our integration program is written in C using the GNU Scientific Library [25] and CUBA. Its source code can be found at https://hg.tx97.net/rambo, and also in the auxiliary files on the arXiv. With a requested accuracy of 0.1% the complete integration takes less than two days on a 12-core machine, with each integration taking between a minute and two hours.
An alternative to the RAMBO-like algorithm, we have also implemented the generation of phase space points by a sequence of two-particle decays. This method is computationally faster,  but does not give a flat distribution of points over the phase space. Still, the convergence of Vegas integration with this method appeared to be as good as with the RAMBO version. We do not supply a corresponding table of results, because it is very similar to Table 2.
We would like to note that such numerical integration is commonly done differently, e.g. via the method sector decomposition [26]. To apply that method to phase space integrals one needs to parameterize the phase space by a suitable rational mapping from a hypercube, and for five particles we found this to be challenging, whereas a four-particle phase-space mapping is known from [2].

Conclusions and outlook
In this article we present analytical expressions for five-particle phase-space integrals expressed in terms of multiple zeta values up to weight 12. The results are calculated using dimensional recurrence relation method with a 2000-digit accuracy using the DREAM package, and analytically restored via the PSLQ algorithm. We also present a computer code for the numerical integration of phase-space integrals in a higher-number of dimensions that has been used to cross-check the obtained results. The approach presented here shows excellent performance for calculating singlescale integrals without ultra-violet divergences and can be easily applied to other problems of this kind.
Further work in this direction includes the calculation of multi-scale phase-space integrals needed for the extraction of time-like splitting functions, with the presented integrals used as boundary conditions of differential equations. Also of interest is the calculation of the remaining unknown integrals with 2-, 3-and 4-particle cuts appearing in the optical theorem for the four-loop propagator.