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

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.

The mainstream strategy to calculate divergent phase-space integration is to divide integrals into singular part and finite part, so that the first part can be calculated easily (either analytically or numerically) and the second part can be calculated purely numerically using Monte Carlo GehrmannDeRidder:2005cm; Currie:2013vh; Czakon:2010td; Czakon:2011ve; Boughezal:2011jf; Catani:1996vz; DelDuca:2015zqa; Harris:2001sx; Binoth:2000ps; Borowka:2017idc; Cacciari:2015jma; Catani:2007vq; Boughezal:2015eha; Boughezal:2015dva; Gaunt:2015pea; Caola:2017dug; Magnea:2018hab; Magnea:2018ebr. If the process in consideration is sufficient inclusive, one can map phase-space integrals onto corresponding loop integrals by using the reverse unitarity relation Anastasiou:2002yz; Anastasiou:2002qz; Anastasiou:2003yy

where $\mathcal{D}^{\mathrm{c}}_{i}=k_{i}^{2}-m_{i}^{2}$ can be interpreted either as mass shell condition or as inverse propagator on cut. In this way, techniques developed for calculating loop integration can be used, like integration-by-parts (IBP) relations Chetyrkin:1981qh, differential equations Kotikov:1990kg; Gehrmann:1999as, dimensional recurrence relations Tarasov:1996br; Lee:2009dh; Lee:2017ftw, and also methods developed by introducing auxiliary mass (AM) Liu:2017jxz; Liu:2018dmc; Zhang:2018mlo; Wang:2019mnn; Guan:2019bcx; Yang:2020msy; Bronnum-Hansen:2020mzk. Further more, loop integration and phase-space integration can be dealt with as a whole because they are not significantly different from each other for these techniques.

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^{+}\geq 0$, $L^{-}\geq 0$, $M\geq 1$ and $N\geq 1$ are understandable. We denote $Q$ as the total cut momentum which satisfies $Q=\sum_{i=1}^{M}q_{i}=\sum_{i=1}^{N}k_{i}$ if $\{q_{i}\}$ 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 $\vec{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\epsilon$ is introduced to regularize all possible divergences, $\mathcal{D}^{\mathrm{t}}_{\alpha}$ are inverse of tree propagators which do not depend on loop momenta, $\mathcal{D}^{+}_{\alpha}$ 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, $\mathcal{D}^{-}_{\alpha}$ 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 $\vec{\nu}\equiv({\nu}^{\mathrm{t}}_{1},{\nu}^{\mathrm{t}}_{2},\cdots,{\nu}^{+}_{1},{\nu}^{+}_{2},\cdots,{\nu}^{-}_{1},{\nu}^{-}_{2},\cdots,\nu^{\pm}_{11},\nu^{\pm}_{12},\nu^{\pm}_{21},\cdots)$ with $\nu_{ij}^{\pm}\leq 0$, the vector $\vec{\eta}\equiv({\eta}^{\mathrm{t}}_{1},{\eta}^{\mathrm{t}}_{2},\cdots,{\eta}^{+}_{1},{\eta}^{+}_{2},\cdots,{\eta}^{-}_{1},{\eta}^{-}_{2},\cdots)$ denotes the introduced AM terms, and $\text{ dPS}_{N}$ is the measure of $N$-particle-cut phase-space integration. For total cross section, we have

where $\Theta$ is the Heaviside function. ¹¹1Note that we use $\Theta\!(k_{i}^{0}-m_{i})$ instead of usual $\Theta\!(k_{i}^{0})$ here. Although they are equivalent in the definition of Eq. (3), the advantage of $\Theta\!(k_{i}^{0}-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. Differential cross sections can be obtained by introducing constraints into $\text{ 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 $\eta^{+}_{\alpha}$ and $\eta^{-}_{\alpha}$ 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 $\eta^{\mathrm{t}}_{\alpha}$, 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 $\eta^{\mathrm{t}}_{\alpha}$ 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. ²²2Positive (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 $\eta^{\mathrm{t}}_{\alpha}$ on the two sides of the cut and take it to $\mathrm{i}0^{+}$ or $-\mathrm{i}0^{+}$ respectively.

For our purpose, we can choose components of $\vec{\eta}$ to be either fully related to each others or completely independent. An extreme is to choose all components of $\vec{\eta}$ to be the same, and an opposite extreme is to choose a strong ordering for all components of $\vec{\eta}$. Although all these choices are workable, to be definite in this work we assume that components of $\vec{\eta}$ can only be either $0^{+}$ or $\eta$. In this way, $F$ only depends on one AM $\eta$, and we denote it as ${F}(\vec{\nu};\vec{s},\eta)$ in the rest of this paper.

In this work, we study the calculation of physical ${F}(\vec{\nu};\vec{s},0)$ based on the auxiliary-mass-flow (AMF) method originally proposed in Ref. Liu:2017jxz for pure loop integration, where flow of $\eta$ from $\infty$ to $0^{+}$ is obtained by solving differential equations w.r.t. $\eta$. 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 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^{-}\rightarrow\gamma^{*}\rightarrow t\bar{t}+X$ at NNLO. We also verify the correctness of our calculation by various methods. Finally, a summary is given in Sec.LABEL:sec:pssum. The calculation of basal phase-space integration without denominator in the integrand is given in Appendix LABEL:sec:pscut.

The advantage of introducing $\eta$ is that, by taking $\eta\to\infty$, ${F}(\vec{\nu};\vec{s},\eta)$ 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 $\eta$ 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 Smirnov:1990rz; Smirnov:1994tg imply that, as $\eta\to\infty$, linear combinations of loop momenta can be either at the order of ${|\eta|}^{1/2}$ or much smaller than it. Therefore one can do the following AME,

where we decompose $\mathcal{D}^{+}_{\alpha}=\widetilde{\mathcal{D}}^{+}_{\alpha}+K_{\alpha}$ with $\widetilde{\mathcal{D}}^{+}_{\alpha}$ including only the part at the order of ${|\eta|}$. 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 $\eta$ presents in the propagator and $\widetilde{\mathcal{D}}^{+}_{\alpha}=0$) or decouples some loop momenta at the order of ${|\eta|}^{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 $\eta\to\infty$, ${F}(\vec{\nu};\vec{s},\eta)$ 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 $\eta$ and take $\eta\to\infty$. Eventually, ${F}(\vec{\nu};\vec{s},\eta)$ is translated to the following form

where $c$ are rational functions of $\vec{s}$ and $\eta$, $F^{\mathrm{bub}}$ include only single-scale vacuum bubble integrals, and $F^{\mathrm{cut}}$ denote basal phase-space integrations with integrand being polynomials of scalar products between cut momenta. $F^{\mathrm{bub}}$ have been well studied up to five-loop order Davydychev:1992mt; Broadhurst:1998rz; Kniehl:2017ikj; Schroder:2005va; Luthe:2015ngq; Luthe:2017ttc. $F^{\mathrm{cut}}$ can be easily dealt with because the only nontrivial information are the cut propagators, which will be explicit studied in Appendix.LABEL:sec:pscut. With these information in hand, the next question is how to obtain physical integrals.

It was shown in Ref. Smirnov:2010hn 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 Laporta:2001dd; Anastasiou:2004vj; Smirnov:2019qkx; Maierhofer:2018gpa; vonManteuffel:2012np; Lee:2013mka. Alternatively, one can achieve IBP reduction using finite-field interpolation vonManteuffel:2014ixa; Peraro:2016wsq; Klappert:2019emp; Klappert:2020aqs; Klappert:2020nbg; Peraro:2019svx, module intersection Boehm:2018fpv, intersection theory Mastrolia:2018uzb, or AME Liu:2018dmc. ³³3 IBP reduction can be achieved by AME for the specific ${F}(\vec{\nu};\vec{s},\eta)$ with $\eta$ introduced to all inverse propagators (not for cut propagators), which is a generalization of the method for pure loop integration introduced in Ref. Liu:2018dmc. With this strategy, coefficients of the expansion are polynomials of kinematical invariants. In any case, the search algorithm proposed in Refs. Liu:2018dmc; Guan:2019bcx 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 Kotikov:1990kg; Bern:1992em; Remiddi:1997ny; Gehrmann:1999as; Henn:2013pwa; Lee:2014ioa. Especially, DEs w.r.t. the AM $\eta$ are given by

where $\vec{J}(\vec{s},\eta)\equiv\left\{{F}(\vec{\nu}^{\prime},\vec{s},\eta),{F}(\vec{\nu}^{\prime\prime},\vec{s},\eta),\cdots\right\}$ is a complete set of MIs and $M(\vec{s},\eta)$ is the coefficient matrix as rational function of $\vec{s}$ and $\eta$. Boundary condition of the DEs can be chosen at $\eta\to\infty$, which can be easily obtained by the AME discussed above. By solving the above DEs (usually numerically) one can realize the flow of $\eta$ from $\infty$ to $0^{+}$. In this way, we get a general method to calculate physical MIs $\vec{J}(\vec{s},0)$ with any fixed $\vec{s}$.

Furthermore, in the case of more than one kinematical invariant, the AMF method can also be combined with DEs w.r.t. $\vec{s}$ to obtain MIs at different values of $\vec{s}$.

As a simple but nontrivial example, we calculate MIs encountered in the NNLO correction for $t\bar{t}$ 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_{t}^{2}$) besides $\eta$. Thus we can introduce dimensionless integrals

where $s=Q^{2}$, $x=\frac{4m_{t}^{2}}{s}$, $y=\frac{\eta}{s}$ and $\nu$ is the summation of all components of $\vec{\nu}$. Because the problem is simple, we make the following unoptimized scheme choice:

with $\eta_{1}^{\mathrm{t}}\ll\eta_{1}^{+}$. More precisely, if $\{\mathcal{D}^{-}_{\alpha}\}$ or $\{\mathcal{D}^{+}_{\alpha}\}$ depend on $\vec{s}$, we choose $\eta^{\mathrm{t}}_{1}=0^{+}$ and $\eta^{+}_{1}=\eta$; Otherwise, we choose $\eta^{\mathrm{t}}_{1}=\eta$ (introduction of $\eta_{1}^{+}$ is unnecessary in this case). The publicly available systematic package FIRE6 Smirnov:2019qkx is sufficient to do all needed reductions in this simple exercise.

There are many subprocesses at the partonic level. We use $\mathrm{V}^{L^{+}}\mathrm{R}^{N-1}\mathrm{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.1. Calculation of MIs for 4-particle cuts at NNLO (RRR) will be presented in Sec.III.2. Calculation of MIs for 3-particle cuts at NNLO (VRR) will be presented in Sec.LABEL:sec:VRRR. And calculation of MIs for 2-particle cuts at NNLO (VVR or VRV) will be presented in Sec.LABEL:sec:VVRR. Verification of these results will be given in Sec.LABEL:sec:check. We note that all these results have already been calculated in literature using other methods (see, e.g., Refs. Bernreuther:2011jt; Bernreuther:2013uma; Dekkers:2014hna; Magerya:2019cvz) although they are not fully publicly available. We provide our results as ancillary file in the arXiv version.