Calculation of Feynman loop integration and phase-space integration via auxiliary mass flow

We extend the auxiliary-mass-flow (AMF) method originally developed for Feynman loop integration to calculate integrals involving also phase-space integration. Flow of the auxiliary mass from the boundary ($\infty$) to the physical point ($0^+$) is obtained by numerically solving differential equations with respective to the auxiliary mass. For problems with two or more kinematical invariants, the AMF method can be combined with traditional differential equation method by providing systematical boundary conditions and highly nontrivial self-consistent check. The method is described in detail with a pedagogical example of $e^+e^-\rightarrow \gamma^* \rightarrow t\bar{t}+X$ at NNLO. We show that the AMF method can systematically and efficiently calculate integrals to high precision.


I. INTRODUCTION
With the good performance of the LHC, the particle physics enters an era of precision measurement. To further test the particle physics standard model and to probe new physics, theoretical calculation at high order in the framework of perturbative quantum field theory is needed to match the precision of experimental data. One of the main difficulty for high-order calculation is the phase-space integration. On the one hand, usually there are soft and collinear divergences under integration which makes it impossible to calculate phase-space integration directly using Monte Carlo numerical method. On the other hand, in general it is hard to express results in terms of known analytical special functions. Significant progresses have been obtained in the past decades.
To be definite, a schematic cut diagram for a general process is shown in Fig.1, where L + is the number of loop momenta (denoting as {l + i }) on the l.h.s. of the cut, L − is the number of loop momenta (denoting as {l − i }) on the r.h.s. of the cut, L = L + + L − is the total number of loop momenta, M is the number of external momenta (denoting as {q i }) which contains not only initial state external momenta but also fixed and unintegrated final state momenta, and N is the number of cut momenta (denoting as {k i }) which are on mass shell with m i being corresponding particle mass. L + ≥ 0, L − ≥ 0, M ≥ 1 and N ≥ 1 are understandable. We denote Q as the total cut momentum which are labeled to flow into the diagram and {k i } flow out of the diagram. A complete set of kinematical invariants after performing loop integration and phase-space integration is denoted as s with Q 2 as a special component. A general phase-space and loop integration with AM to be studied in this work is where the non-integer spacetime dimension D = 4 − 2 is introduced to regularize all possible divergences, D t α are inverse of tree propagators which do not depend on loop momenta, D + α are inverse of loop propagators on the l.h.s. of the cut which do not depend on loop momenta on the r.h.s. of the cut, D − α are inverse of loop propagators on the r.h.s. of the cut which do not depend on loop momenta on the l.h.s. of the cut, the vector ν ≡ · · · ) denotes the introduced AM terms, and dPS N is the measure of N -particle-cut phase-space integration. For total cross section, we have where Θ is the Heaviside function. 1 Differential cross sections can be obtained by introducing constraints into dPS N . The corresponding physical integral can be obtained from the above modified integral by taking all AMs to zero, It is understandable to take η + α and η − α to zero from the positive side of their real parts, because this is exactly the rule of Feynman prescription for Feynman propagators which guarantees the correct discontinuity of Feynman propagators. While for η t α , we can take them to zero from any direction as far as all tree propagators are either positive-definite or negative-definite. Our choice is to take all η t α to zero from the positive side so that tree propagators on the l.h.s. of the cut can be combined with the same propagators on the r.h.s. of the cut. 2 For our purpose, we can choose components of η to be either fully related to each others or completely independent. An extreme is to choose all components of η to be the same, and an opposite extreme is to choose a strong ordering for all components of η. Although all these choices are workable, to be definite in this work we assume that components of η can only be either 0 + or η. In this way, F only depends on one AM η, and we denote it as F ( ν; s, η) in the rest of this paper.
In this work, we study the calculation of physical F ( ν; s, 0) based on the auxiliary-mass-flow (AMF) method originally proposed in Ref. [28] for pure loop integration, where flow of η from ∞ to 0 + is obtained by solving differential equations w.r.t. η. We will see that this method is not only systematical and efficient, but also can give high-precision result. The rest of the paper is organized as follows. In Sec.II, we describe the general strategy to 1 Note that we use Θ(k 0 i − m i ) instead of usual Θ(k 0 i ) here. Although they are equivalent in the definition of Eq. (3), the advantage of Θ(k 0 i − m i ) is that its derivative (being Dirac delta function) can be safely set to zero in dimensional regularization. It is guaranteed by the fact that Dirac delta function restricts all space components of k i to be at the origin, where is well regularized by dimensional regularization. Therefore, our choice is convenient to use inverse unitarity. 2 Positive (or negative) definiteness of tree propagators is alway satisfied in the narrow-width approximation, where particle production and decay are factorized and can be calculated separately. Otherwise, we should distinguish η t α on the two sides of the cut and take it to i0 + or −i0 + respectively. calculate integrals involving both loop integration and phase-space integration. In Sec.III, the method is explained in detail by a pedagogical example e + e − → γ * → tt + X at NNLO. We also verify the correctness of our calculation by various methods. Finally, a summary is given in Sec.IV. The calculation of basal phase-space integration without denominator in the integrand is given in Appendix A.

II. AUXILIARY MASS EXPANSION AND FLOW
The advantage of introducing η is that, by taking η → ∞, F ( ν; s, η) can be reduced to linear combinations of simpler integrals. As scalar products among external momenta and cut momenta are finite, we have the auxiliary mass expansion (AME) for tree propagators which removes a tree propagator from the denominator if η has been introduced to it. Because loop momenta can be any large value, one cannot naively expand loop propagators in the same way as tree propagators. However, the standard rules of large-mass expansion [35,36] imply that, as η → ∞, linear combinations of loop momenta can be either at the order of |η| 1/2 or much smaller than it. Therefore one can do the following AME, 1 where we decompose D + α = D + α + K α with D + α including only the part at the order of |η|. Similarly we can do the expansion for loop propagators on the r.h.s. of the cut. The AME of loop propagators either removes some propagators from the denominator (if η presents in the propagator and D + α = 0) or decouples some loop momenta at the order of |η| 1/2 from kinematical invariants. The later effect results in some single-scale vacuum integrals multiplied by integrals with less number of loop momenta.
We find that, as η → ∞, F ( ν; s, η) is simplified to linear combination of integrals with less inverse propagators in the denominator (maybe multiplied by single-scale vacuum integrals). If the simplified integrals still have inverse propagators in the denominator (except propagators in single-scale vacuum integrals), we can again introduce new AM η and take η → ∞. Eventually, F ( ν; s, η) is translated to the following form where c are rational functions of s and η, F bub include only single-scale vacuum bubble integrals, and F cut denote basal phase-space integrations with integrand being polynomials of scalar products between cut momenta. F bub have been well studied up to five-loop order [37][38][39][40][41][42]. F cut can be easily dealt with because the only nontrivial information are the cut propagators, which will be explicit studied in Appendix.A. With these information in hand, the next question is how to obtain physical integrals. It was shown in Ref. [43] that, for any given problem, Feynman loop integrals form an finite-dimensional vector space, with basis of which called master integrals (MIs). The step to express all loop integrals as linear combinations of MIs are called reduction. With reverse unitarity relation in Eq.(1), one can map phase-space integrations onto corresponding loop integrations. Therefore, integrals over phase-space and loop momenta defined in Eq. (2) can also be reduced to corresponding MIs.
Reduction of a general integral to MIs can be traditionally achieved by using IBP relations based on Laporta's algorithm [44][45][46][47][48][49]. Alternatively, one can achieve IBP reduction using finite-field interpolation [50][51][52][53][54][55], module intersection [56], intersection theory [57], or AME [29]. 3 In any case, the search algorithm proposed in Refs. [29,32] can significantly improve the efficiency of reduction, which makes the reduction of very complicated problems a possibility. Reduction can not only express all integrals in terms of MIs, it can also set up differential equations (DEs) among MIs [23,24,[58][59][60][61]. Especially, DEs w.r.t. the AM η are given by where J( s, η) ≡ {F ( ν , s, η), F ( ν , s, η), · · · } is a complete set of MIs and M ( s, η) is the coefficient matrix as rational function of s and η. Boundary condition of the DEs can be chosen at η → ∞, which can be easily obtained by the AME discussed above. By solving the above DEs (usually numerically) one can realize the flow of η from ∞ to 0 + . In this way, we get a general method to calculate physical MIs J( s, 0) with any fixed s. Furthermore, in the case of more than one kinematical invariant, the AMF method can also be combined with DEs w.r.t. s to obtain MIs at different values of s.
As a simple but nontrivial example, we calculate MIs encountered in the NNLO correction for tt production in e + e − collision mediated by a virtual photon to demonstrate the validity of AMF method. For the purpose of total cross section, there are only two kinematical invariants (Q 2 and m 2 t ) besides η. Thus we can introduce dimensionless integralsF where s = Q 2 , x = 4m 2 t s , y = η s and ν is the summation of all components of ν. Because the problem is simple, we make the following unoptimized scheme choice: we choose η t 1 = 0 + and η + 1 = η; Otherwise, we choose η t 1 = η (introduction of η + 1 is unnecessary in this case). The publicly available systematic package FIRE6 [46] is sufficient to do all needed reductions in this simple exercise.
There are many subprocesses at the partonic level. We use V L + R N −1 V L − to distinguish processes with different number of independent loop momenta and cut momenta. The presence of N −1 in stead of N is a result of momentum conservation which reduces one independent cut momentum. For completeness, we first provide MIs at NLO in Sec.III A. Calculation of MIs for 4-particle cuts at NNLO (RRR) will be presented in Sec.III B. Calculation of MIs for 3-particle cuts at NNLO (VRR) will be presented in Sec.III C. And calculation of MIs for 2-particle cuts at NNLO (VVR or VRV) will be presented in Sec.III D. Verification of these results will be given in Sec.III E. We note that all these results have already been calculated in literature using other methods (see, e.g., Refs. [62][63][64][65]) although they are not fully publicly available. We provide our results as ancillary file in the arXiv version.

A. RR and VR
For RR, i.e., γ * → ttg process, inverse propagators can be written as The obtained 2 physical MIs are F cut 2,3,1 and F cut 2,3,2 , which are calculated in Appendix.A. For VR, inverse propagators can be written as We obtain 1 family 4 with ν of characterized MIs being {1, 1, 0}. Physical MIs are given by where Φ 2 = dPS 2 is defined in Eq.(A5) and the remaining 1-loop integrals are easy to calculate analytically.

B. RRR
For RRR, i.e., γ * → ttgg or ttqq processes, inverse propagators can be written as where Q 2 = s. Then phase-space integrals can be expressed aŝ By using FIRE6, we find that there are 37 MIs for finite η and the number is reduced to 15 as η vanishing.
which can also characterized by Feynman diagrams Fig.2 (a) and (b), respectively. Using our method, we can calculate As usual, we use the corner integral, which has no inverse propagator in the numerator and has the maximal kinds of inverse propagators in the denominator with power of each inverse propagators being unit, to represent the family. 5 Note that the number of MIs in this work may not be the minimal value. It does not matter as far as it is finite. + 0.00589205324016960475844728350032 4 + 0.0286458793127880349046701435743 5 + 0.118038608602644851031978155328 6 + 0.425508191586298756000765241113 7 where we truncate the expansion to order 8 with about 30-digit precision for each coefficient. Physical MIs with other values of x can be calculated similarly. In the following, let us take a sub-family {0, 0, 0, 0, 0, 0, 1} shown in Fig.2 (c) as an example to illustrate the calculation procedure of MIs. For brevity, we define MIs for this family as: with ν t 1 , ν t 2 ≤ 0 and ν t 7 ≥ 0. This family contains 5 MIs for finite η (or y) The boundary condition of these DEs is given bŷ By solving the DEs (23) with boundary condition we obtain, e.g., With the boundary condition at x = 1/2 given in Eq. (25), solving the above DE can also obtainF ({0, 0, 1}; x, 0) at any value of x.
The most complicated family is the first family, which can also be characterized by Feynman diagram Fig.3 (a). Using In the following, let us take a sub-family {1, 0, 1, 1, 0, 0, 0, 0} shown in Fig.3 (b) as an example to illustrate the calculation procedure of MIs. For brevity, we define MIs for this family as: To calculate the last 2 physical MIs, we set up DEs for corresponding MIs with x = 1/2 and finite η, which gives The boundary condition of these DEs is given bŷ where and dPS 3 can again be calculated by the method in Sec.III B or be obtained from RR in Sec.III A. Knowing the first 5 MIs, by solving the DEs (34) with boundary condition we obtain, e.g., Finally, the DEs w.r.t. x is given by With the boundary condition at x = 1/2 given in Eq. (37), solving the above DEs can also evaluateF ({1, 1, 1}; x, 0) andF ({1, 1, 2}; x, 0) at any value of x.
The most complicated family is the first family, which can also be characterized by Feynman diagram Fig.4 (a). We use the same method (note, F bub up to 2-loop is needed to calculate boundary conditions) in VRR to calculate MIs for these families. For VRV, inverse propagators can be written as in addition to which there is a scalar product l + 1 · l − 1 . The obtained family is: which can also be characterized by Feynman diagram Fig.4 (b). Because loop integrations in MIs of this case are factorized, their calculation is as simple as one-loop case.

E. Verification of the results
On the one hand, the AMF method can be used to calculate MIs at any given value of x. On the other hand, MIs at different values of x can be related by DEs w.r.t. x. We have verified that results of MIs for different values of x, e.g., x = 1/2 and x = 2/5, obtained by the two strategies are consistent with each other. This provides a highly nontrivial self-consistent check because DEs w.r.t. η and that w.r.t. x are significantly different.
Our numerical values forF ( ν; 0, 0), i.e. MIs for massless QCD, are in full agreement with the known analytical results of massless MIs in literature [66,67].

IV. SUMMARY AND OUTLOOK
In this paper, we extend the AMF method originally developed for Feynman loop integrals [28] to calculate MIs involving also phase-space integration. As a pedagogical example, we use it to calculate MIs encountered in e + e − → γ * → tt + X at NNLO. Our results agree with results obtained by using other methods (our fully results are available as ancillary file in the arXiv version). Although AMF method depends on reduction procedure to decompose all integrals to MIs and to set up DEs of MIs w.r.t. η, the efficiency of reduction has been significantly improved thanks to the recently proposed search algorithm [29,32].
It is clear that the AMF method can be used to calculate MIs of any process, as systematic as the sector decomposition method. However, comparing with the later method, AMF method is much more efficient and can provide very high precision. The high-precision nature makes it possible to obtain analytical results with proper ansatz. For problems with two or more kinematical invariants, where differential equation method works, the AMF method can not only systematically provide as many as needed boundary conditions for DEs w.r.t. kinematical invariants, but also provide highly nontrivial self-consistent check for the obtained results. All the above advantages make the AMF method very useful for perturbative calculation at high orders. r = 0, 1, and 2. MIs for general cases can be studied as following. By using optical theorem (see, e.g., Ref. [66]), the calculation of basal phase-space integral is translated to the calculation of imaginary part of the corresponding sunrise pure loop integral. The later can be calculated by the AMF method for loop integrals [28]. Furthermore, a one-dimensional-integral representation of all these MIs can be obtained from Ref. [68].
For r = 0 and N ≥ 2, theres is only one MI, which is given by For r = 1, in general there are two MIs, which are given by (n = 1, 2) where p F q are generalized hypergeometric functions which can be evaluated by using the public program HypExp [69]. In the special case of N = 2, there is only one MI F cut 1,2,1 . For r = 2, in general there are three MIs, which are given by (n = 1, 2, 3) In the special case of N = 2, there are only one MI, which means F cut 2,2,2 and F cut 2,2,3 can be reduced to F cut 2,2,1 . In the special case of N = 3, there are only two MIs, which means F cut 2,3,3 can be reduced to F cut 2,3,1 and F cut 2,3,2 . The number of MIs obtained here is consistent with the expectation in Ref. [70]. The results (A2) and (A3) agree with partial results in literature [see Eq.(B.7) in Ref. [71] and Eqs. (3.14-3.16) in Ref. [62]].
To obtain (A2) and (A3), we decompose the measure of phase-space integration to two parts: where Ω is a subset of {1, 2, . . . , N } and # Ω is the number of elements of Ω. Specifically, the Ω is chosen to be the set of massless particles for r = 1, and is chosen to be the set of massive particles for r = 2. The phase-space volume of two-particle cut is used: (A8)