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 which also involve phase-space integration. The 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 the traditional differential-equation method, providing systematic boundary conditions and a highly nontrivial self-consistency check. The method is described in detail using 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.
  • 加载中
  • [1] A. Gehrmann-De Ridder, T. Gehrmann, and E. W. N. Glover, JHEP 09, 056 (2005), arXiv:hep-ph/0505111
    [2] J. Currie, E. Glover, and S. Wells, JHEP 04, 066 (2013), arXiv:1301.4693
    [3] M. Czakon, Phys. Lett. B 693, 259-268 (2010), arXiv:1005.0274 doi: 10.1016/j.physletb.2010.08.036
    [4] M. Czakon, Nucl. Phys. B 849, 250-295 (2011), arXiv:1101.0642 doi: 10.1016/j.nuclphysb.2011.03.020
    [5] R. Boughezal, K. Melnikov, and F. Petriello, Phys. Rev. D 85, 034025 (2012), arXiv:1111.7041 doi: 10.1103/PhysRevD.85.034025
    [6] S. Catani and M. H. Seymour, A General algorithm for calculating jet cross-sections in NLO QCD, Nucl. Phys. B 485, 291-419 (1997) [hep-ph/9605323] [InSPIRE], [Erratum: Nucl. Phys. B 510, 503 (1998)]
    [7] V. Del Duca, C. Duhr, G. Somogyi et al., JHEP 04, 036 (2015), arXiv:1501.07226
    [8] B. W. Harris and J. F. Owens, Phys. Rev. D 65, 094032 (2002), arXiv:hep-ph/0102128
    [9] T. Binoth and G. Heinrich, Nucl. Phys. B 585, 741-759 (2000), arXiv:hep-ph/0004013 doi: 10.1016/S0550-3213(00)00429-6
    [10] S. Borowka, G. Heinrich, S. Jahn et al., Comput. Phys. Commun. 222, 313-326 (2018), arXiv:1703.09692 doi: 10.1016/j.cpc.2017.09.015
    [11] M. Cacciari, F. A. Dreyer, A. Karlberg et al., Fully Differential Vector-Boson-Fusion Higgs Production at Next-to-Next-to-Leading Order, Phys. Rev. Lett. 115, 082002 (2015), arXiv: 1506.02660 [InSPIRE], [Erratum: Phys. Rev. Lett. 120, 139901 (2018)]
    [12] S. Catani and M. Grazzini, Phys. Rev. Lett. 98, 222002 (2007), arXiv:hep-ph/0703012 doi: 10.1103/PhysRevLett.98.222002
    [13] R. Boughezal, X. Liu, and F. Petriello, Phys. Rev. D 91, 094035 (2015), arXiv:1504.02540 doi: 10.1103/PhysRevD.91.094035
    [14] R. Boughezal, C. Focke, X. Liu et al., Phys. Rev. Lett. 115, 062002 (2015), arXiv:1504.02131 doi: 10.1103/PhysRevLett.115.062002
    [15] J. Gaunt, M. Stahlhofen, F. J. Tackmann et al., JHEP 09, 058 (2015), arXiv:1505.04794
    [16] F. Caola, K. Melnikov, and R. Röntsch, Eur. Phys. J. C 77, 248 (2017), arXiv:1702.01352 doi: 10.1140/epjc/s10052-017-4774-0
    [17] L. Magnea, E. Maina, G. Pelliccioli et al., 12, 107 (2018), arXiv: 1806.09570 [InSPIRE], [Erratum: JHEP 06, 013 (2019)]
    [18] L. Magnea, E. Maina, G. Pelliccioli et al., JHEP 12, 062 (2018), arXiv:1809.05444
    [19] C. Anastasiou and K. Melnikov, Nucl. Phys. B 646, 220-256 (2002), arXiv:hep-ph/0207004
    [20] C. Anastasiou, L. J. Dixon, and K. Melnikov, Nucl. Phys. Proc. Suppl. 116, 193-197 (2003), arXiv:hep-ph/0211141 doi: 10.1016/S0920-5632(03)80168-8
    [21] C. Anastasiou, L. J. Dixon, K. Melnikov et al., Phys. Rev. Lett. 91, 182002 (2003), arXiv:hep-ph/0306192 doi: 10.1103/PhysRevLett.91.182002
    [22] K. Chetyrkin and F. Tkachov, Nucl. Phys. B 192, 159-204 (1981) doi: 10.1016/0550-3213(81)90199-1
    [23] A. Kotikov, Phys. Lett. B 254, 158-164 (1991) doi: 10.1016/0370-2693(91)90413-K
    [24] T. Gehrmann and E. Remiddi, Nucl. Phys. B 580, 485-518 (2000), arXiv:hep-ph/9912329 doi: 10.1016/S0550-3213(00)00223-6
    [25] O. Tarasov, Phys. Rev. D 54, 6479-6490 (1996), arXiv:hep-th/9606018 doi: 10.1103/PhysRevD.54.6479
    [26] R. Lee, Nucl. Phys. B 830, 474-492 (2010), arXiv:0911.0252 doi: 10.1016/j.nuclphysb.2009.12.025
    [27] R. N. Lee and K. T. Mingulov, DREAM, a program for arbitrary-precision computation of dimensional recurrence relations solutions, and its applications, arXiv: 1712.05173 [InSPIRE]
    [28] X. Liu, Y.-Q. Ma, and C.-Y. Wang, Phys. Lett. B 779, 353-357 (2018), arXiv:1711.09572
    [29] X. Liu and Y.-Q. Ma, Phys. Rev. D 99, 071501 (2019), arXiv:1801.10523
    [30] P. Zhang, C.-Y. Wang, X. Liu et al., JHEP 04, 116 (2019), arXiv:1810.07656
    [31] Y. Wang, Z. Li, and N. Ul Basat, Phys. Rev. D 101, 076023 (2020), arXiv:1901.09390 doi: 10.1103/PhysRevD.101.076023
    [32] X. Guan, X. Liu, and Y.-Q. Ma, Complete reduction of two-loop five-light-parton scattering amplitudes, arXiv: 1912.09294 [InSPIRE]
    [33] Q. Yang, R.-Y. Zhang, M.-M. Long et al., Chin. Phys. C 44, 093101 (2020), arXiv:2005.11010 doi: 10.1088/1674-1137/44/9/093101
    [34] C. Brønnum-Hansen and C.-Y. Wang, Contribution of third generation quarks to two-loop helicity amplitudes for W boson pair production in gluon fusion, arXiv: 2009.03742 [InSPIRE]
    [35] V. A. Smirnov, Commun. Math. Phys. 134, 109-137 (1990) doi: 10.1007/BF02102092
    [36] V. A. Smirnov, Mod. Phys. Lett. A 10, 1485-1500 (1995), arXiv:hep-th/9412063 doi: 10.1142/S0217732395001617
    [37] A. I. Davydychev and J. Tausk, Nucl. Phys. B 397, 123-142 (1993) doi: 10.1016/0550-3213(93)90338-P
    [38] D. J. Broadhurst, Eur. Phys. J. C 8, 311-333 (1999), arXiv:hep-th/9803091 doi: 10.1007/s100529900935
    [39] http://theor.jinr.ru/~kalmykov/three-loop/three.html#Master
    [40] B. Kniehl, A. Pikelner, and O. Veretin, JHEP 08, 024 (2017), arXiv:1705.05136
    [41] Y. Schroder and A. Vuorinen, JHEP 06, 051 (2005), arXiv:hep-ph/0503209
    [42] T. Luthe, Fully massive vacuum integrals at 5 loops. PhD thesis, Bielefeld U., 2015 [InSPIRE]
    [43] T. Luthe, A. Maier, P. Marquard et al., JHEP 03, 020 (2017), arXiv:1701.07068
    [44] A. Smirnov and A. Petukhov, Lett. Math. Phys. 97, 37-44 (2011), arXiv:1004.4199 doi: 10.1007/s11005-010-0450-0
    [45] S. Laporta, Int. J. Mod. Phys. A 15, 5087-5159 (2000), arXiv:hep-ph/0102033
    [46] C. Anastasiou and A. Lazopoulos, JHEP 07, 046 (2004), arXiv:hep-ph/0404258
    [47] A. V. Smirnov and F. S. Chuharev, FIRE6: Feynman Integral REduction with Modular Arithmetic, arXiv: 1901.07808 [InSPIRE]
    [48] P. Maierhöfer and J. Usovitsch, Kira 1.2 Release Notes, arXiv: 1812.01491 [InSPIRE]
    [49] A. von Manteuffel and C. Studerus, Reduze 2 - Distributed Feynman Integral Reduction, arXiv: 1201.4330 [InSPIRE]
    [50] R. N. Lee, J. Phys. Conf. Ser. 523, 012059 (2014), arXiv:1310.1145 doi: 10.1088/1742-6596/523/1/012059
    [51] A. von Manteuffel and R. M. Schabinger, Phys. Lett. B 744, 101-104 (2015), arXiv:1406.4513 doi: 10.1016/j.physletb.2015.03.029
    [52] T. Peraro, JHEP 12, 030 (2016), arXiv:1608.01902
    [53] J. Klappert and F. Lange, Comput. Phys. Commun. 247, 106951 (2020), arXiv:1904.00009 doi: 10.1016/j.cpc.2019.106951
    [54] J. Klappert, S. Y. Klein, and F. Lange, Interpolation of Dense and Sparse Rational Functions and other Improvements in FireFly, arXiv: 2004.01463 [InSPIRE]
    [55] J. Klappert, F. Lange, P. Maierhöfer et al., Integral Reduction with Kira 2.0 and Finite Field Methods, arXiv: 2008.06494 [InSPIRE]
    [56] T. Peraro, JHEP 07, 031 (2019), arXiv:1905.08019
    [57] J. Böhm, A. Georgoudis, K. J. Larsen et al., JHEP 09, 024 (2018), arXiv:1805.01873
    [58] P. Mastrolia and S. Mizera, JHEP 02, 139 (2019), arXiv:1810.03818
    [59] Z. Bern, L. J. Dixon, and D. A. Kosower, Dimensionally regulated one loop integrals, Phys. Lett. B 302, 299-308 (1993) [hep-ph/9212308] [InSPIRE], [Erratum: Phys. Lett. B 318, 649 (1993)]
    [60] E. Remiddi, Nuovo Cim. A 110, 1435-1452 (1997), arXiv:hep-th/9711188
    [61] J. M. Henn, Phys. Rev. Lett. 110, 251601 (2013), arXiv:1304.1806 doi: 10.1103/PhysRevLett.110.251601
    [62] R. N. Lee, JHEP 04, 108 (2015), arXiv:1411.0911
    [63] W. Bernreuther, C. Bogner, and O. Dekkers, JHEP 06, 032 (2011), arXiv:1105.0530
    [64] W. Bernreuther, C. Bogner, and O. Dekkers, JHEP 10, 161 (2013), arXiv:1309.6887
    [65] O. Dekkers and W. Bernreuther, Phys. Lett. B 738, 325-333 (2014), arXiv:1409.3124 doi: 10.1016/j.physletb.2014.09.060
    [66] V. Magerya and A. Pikelner, JHEP 12, 026 (2019), arXiv:1910.07522
    [67] A. Gehrmann-De Ridder, T. Gehrmann, and G. Heinrich, Nucl. Phys. B 682, 265-288 (2004), arXiv:hep-ph/0311276
    [68] A. Gehrmann-De Ridder, T. Gehrmann, and E. Glover, Nucl. Phys. B 691, 195-222 (2004), arXiv:hep-ph/0403057 doi: 10.1016/j.nuclphysb.2004.05.017
    [69] S. Groote, J. Korner, and A. Pivovarov, Annals Phys. 322, 2374-2445 (2007), arXiv:hep-ph/0506286 doi: 10.1016/j.aop.2006.11.001
    [70] T. Huber and D. Maitre, Comput. Phys. Commun. 178, 755-776 (2008), arXiv:0708.2443 doi: 10.1016/j.cpc.2007.12.008
    [71] M. Y. Kalmykov and B. A. Kniehl, JHEP 07, 031 (2017), arXiv:1612.06637
    [72] C. Anastasiou, C. Duhr, F. Dulat et al., JHEP 07, 003 (2013), arXiv:1302.4379
    [73] A. P. Prudnikov, Y. A. Brychkov, and O. I. Marichev, Integrals and Series: More Special Functions, Taylor and Francis Ltd. vol.3 (1989) [InSPIRE]
  • 加载中

Figures(4)

Get Citation
Xiao Liu, Yan-Qing Ma, Wei Tao and Peng Zhang. Calculation of Feynman loop integration and phase-space integration via auxiliary mass flow[J]. Chinese Physics C. doi: 10.1088/1674-1137/abc538
Xiao Liu, Yan-Qing Ma, Wei Tao and Peng Zhang. Calculation of Feynman loop integration and phase-space integration via auxiliary mass flow[J]. Chinese Physics C.  doi: 10.1088/1674-1137/abc538 shu
Milestone
Received: 2020-09-03
Article Metric

Article Views(1134)
PDF Downloads(58)
Cited by(0)
Policy on re-use
To reuse of Open Access content published by CPC, for content published under the terms of the Creative Commons Attribution 3.0 license (“CC CY”), the users don’t need to request permission to copy, distribute and display the final published version of the article and to create derivative works, subject to appropriate attribution.
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Email This Article

Title:
Email:

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

    Corresponding author: Xiao Liu, xiao6@pku.edu.cn
    Corresponding author: Yan-Qing Ma, yqma@pku.edu.cn
    Corresponding author: Wei Tao, ygtw@pku.edu.cn
    Corresponding author: Peng Zhang, p.zhang@pku.edu.cn
  • 1. School of Physics and State Key Laboratory of Nuclear Physics and Technology, Peking University, Beijing 100871, China
  • 2. Center for High Energy Physics, Peking University, Beijing 100871, China
  • 3. Collaborative Innovation Center of Quantum Matter, Beijing 100871, China

Abstract: We extend the auxiliary-mass-flow (AMF) method originally developed for Feynman loop integration to calculate integrals which also involve phase-space integration. The 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 the traditional differential-equation method, providing systematic boundary conditions and a highly nontrivial self-consistency check. The method is described in detail using 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.

    HTML

    I.   INTRODUCTION
    • With the good performance of the LHC, particle physics has entered an era of precision measurement. To further test the Standard Model of particle physics and to probe new physics, theoretical calculations at high order in the framework of perturbative quantum field theory are needed to match the precision of experimental data. One of the main difficulties for higher-order calculation is the phase-space integration. On the one hand, there are usually soft and collinear divergences under integration which make it impossible to calculate phase-space integration directly using a Monte Carlo numerical method. On the other hand, in general it is hard to express the results in terms of known analytical special functions. Significant progress has been made in the last few decades, however.

      The mainstream strategy to calculate divergent phase-space integration is to divide the integrals into a singular part and a 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 a Monte Carlo method [1-18]. If the process under consideration is sufficient inclusive, one can map phase-space integrals onto corresponding loop integrals by using the reverse unitarity relation [19-21]

      $ (2\pi) \delta({\cal{D}}^{\rm{c}}_i) = \frac{ {\rm{i}}}{{\cal{D}}^{\rm{c}}_i+ {\rm{i}}0^+}+\frac{- {\rm{i}}}{{\cal{D}}^{\rm{c}}_i- {\rm{i}}0^+}, $

      (1)

      where $ {\cal{D}}^{\rm{c}}_i = k_i^2-m_i^2 $ can be interpreted either as a mass shell condition or as an inverse propagator on the cut. In this way, techniques developed for calculating loop integration can be used, such as integration-by-parts (IBP) relations [22], differential equations (DEs) [23, 24], dimensional recurrence relations [25-27], and also methods developed by introducing auxiliary mass (AM) [28-34]. Furthermore, 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 (denoted as $ \{l_i^+\} $) on the l.h.s. of the cut, $ L^- $ is the number of loop momenta (denoted 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 (denoted 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 (denoted as $ \{k_i\} $) which are on mass shell, with $ m_i $ being the corresponding particle masses. $ L^+\geqslant 0 $, $ L^-\geqslant 0 $, $ M\geqslant 1 $ and $ N\geqslant 1 $ are reasonable. We denote Q as the total cut momentumwhich satisfies $ Q = \displaystyle\sum\limits_{i = 1}^M q_i = \displaystyle\sum\limits_{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.

      Figure 1.  Schematic diagram for a process with $ L = L^++L^- $ loops, M unintegrated external legs, and N cut legs.

      A general phase-space and loop integration with AM to be studied in this work is

      $ {F}(\vec{\nu};\vec s, \vec \eta)\equiv\int {\rm{dPS}}_N \prod\limits_{\alpha}\frac{1}{({\cal{D}}^{ \rm{t}}_{\alpha}+\eta^{ \rm{t}}_\alpha)^{\nu^{ \rm{t}}_{\alpha}}} \int \prod\limits_{i = 1}^{L^+} \frac{ {\rm{d}}^Dl_i^+}{(2\pi)^D}\prod\limits_{\beta}\frac{1}{({\cal{D}}^+_{\beta}+ {\rm{i}}\eta^+_\beta)^{\nu^+_{\beta}}} \int \prod\limits_{j = 1}^{L^-}\frac{ {\rm{d}}^Dl_j^-}{(2\pi)^D}\prod\limits_{\gamma}\frac{1}{({\cal{D}}^-_{\gamma}- {\rm{i}}\eta^-_\gamma)^{\nu^-_{\gamma}}} (l_i^+\cdot l_j^-)^{-\nu^{\pm}_{ij}} , $

      (2)

      where the non-integer spacetime dimension $ D = 4-2\epsilon $ is introduced to regularize all possible divergences, $ {\cal{D}}^{ \rm{t}}_{\alpha} $ are the inverse of tree propagators which do not depend on loop momenta, $ {\cal{D}}^+_{\alpha} $ are the 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, $ {\cal{D}}^-_{\alpha} $ are the 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}^{ \rm{t}}_1, {\nu}^{ \rm{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\leqslant 0 $, the vector $ \vec{\eta}\equiv({\eta}^{ \rm{t}}_1,{\eta}^{ \rm{t}}_2,\cdots,{\eta}^+_1,{\eta}^+_2,\cdots,{\eta}^-_1,{\eta}^-_2,\cdots) $ denotes the introduced AM terms, and $ {\rm{ dPS}}_N $ is the measure of N-particle-cut phase-space integration. For the total cross section, we have

      $ {\rm{ dPS}}_N\equiv(2\pi)^{D}\delta^D\left(Q-\sum\limits_{i = 1}^N k_i\right)\prod\limits_{i = 1}^{N}\frac{ {\rm{d}}^Dk_i}{(2\pi)^{D}}(2\pi)\delta({\cal{D}}^{\rm{c}}_i)\Theta\!(k_{i}^{0}-m_i), $

      (3)

      where $ \Theta $ is the Heaviside function. Differential cross sections can be obtained by introducing constraints into $ {\rm{ dPS}}_N $, e.g. an introduction of $ \delta(y-k_i\cdot p_1/k_i\cdot p_2) $ can give rise to the rapidity distribution of the i-th particle [20, 21], where $ p_1 $ and $ p_2 $ are the momenta of initial-state particles.

      The corresponding physical integral can be obtained from the above modified integral by taking all AMs to zero,

      $ {F}(\vec{\nu};\vec s,0)\equiv \mathop {\lim }\limits_{\vec \eta \to {0^ + }} {F}(\vec{\nu};\vec s,\vec \eta). $

      (4)

      It is reasonable 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. For $ \eta^{ \rm{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^{ \rm{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.

      For our purpose, we can choose components of $ \vec \eta $ to be either fully related to each others or completely independent. One extreme is to choose all the components of $ \vec \eta $ to be the same, and the 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. [28] for pure loop integration, where flow of $ \eta $ from $ \infty $ to $ 0^+ $ is obtained by solving DEs w.r.t. $ \eta $ . We will see that this method is not only systematic and efficient, but can also give high-precision results. 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 using the 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. IV. The calculation of basal phase-space integration without a denominator in the integrand is given in Appendix A.

    II.   AUXILIARY MASS EXPANSION AND FLOW
    • 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

      $ \frac{1}{{\cal{D}}^{ \rm{t}}_\alpha+\eta}\; {\xlongequal[\;] {\eta \to \infty }} \frac{1}{\eta}\sum\limits_{j = 0}^{+\infty}\left(\frac{-{\cal{D}}^{ \rm{t}}_\alpha}{\eta}\right)^j, $

      (5)

      $ \frac{1}{{\cal{D}}^{ \rm{t}}_\alpha}\; {\xlongequal[\;] { \eta \to \infty }} \; \frac{1}{{\cal{D}}^{ \rm{t}}_\alpha}, $

      (6)

      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 [35, 36] 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 perform the following AME,

      $ \dfrac{1}{{\cal{D}}^+_\alpha+ {\rm{i}} \eta}\; \xlongequal[\;]{\eta \to \infty} \dfrac{1}{\widetilde{{\cal{D}}}^+_{\alpha}+ {\rm{i}}\eta}\sum\limits_{j = 0}^{+\infty}\left(\frac{-K_\alpha}{\widetilde{{\cal{D}}}^+_{\alpha}+ {\rm{i}}\eta}\right)^j, $

      (7)

      $ \dfrac{1}{{\cal{D}}^+_\alpha\!+\! {\rm{i}} 0^+}\xlongequal[\;]{\eta \to \infty} \left\{\!\!\!\!\!{\begin{array}{*{20}{l}} {\dfrac{1}{\widetilde{{\cal{D}}}^+_{\alpha}+ {\rm{i}}0^+}\displaystyle\sum\limits_{j = 0}^{+\infty}\left(\frac{-K_\alpha}{\widetilde{{\cal{D}}}^+_{\alpha}+ {\rm{i}}0^+}\right)^j}\!\!&\!\!{\rm{if }}\;\widetilde{{\cal{D}}}^+_{\alpha}\neq0,\\ {\dfrac{1}{{{\cal{D}}}^+_{\alpha}+ {\rm{i}}0^+}}\!\!&\!\!{{\rm{if}} \;\widetilde{{\cal{D}}}^+_{\alpha} = 0,} \end{array}} \right. $

      (8)

      where we decompose $ {\cal{D}}^+_\alpha = \widetilde{{\cal{D}}}^+_{\alpha}+K_\alpha $ with $ \widetilde{{\cal{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{{\cal{D}}}^+_{\alpha} = 0 $) or decouples some loop momenta at the order of $ {|\eta|}^{1/2} $ from kinematical invariants. The latter effect results in some single-scale vacuum integrals being multiplied by integrals with fewer loop momenta.

      We find that, as $ \eta\to \infty $, $ {F}(\vec{\nu};\vec s, \eta) $ is simplified to a linear combination of integrals with fewer 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

      $ {F}(\vec{\nu}; \vec s, \eta)\longrightarrow\sum c\times F^{\rm{cut}} \times F^{\rm{bub}}, $

      (9)

      where c are rational functions of $ \vec {s\,} $ and $ \eta $, $ F^{\rm{bub}} $ includes only single-scale vacuum bubble integrals, and $ F^{\rm{cut}} $ denotes basal phase-space integrations with the integrands being polynomials of scalar products between cut momenta. $ F^{\rm{bub}} $ have been well studied up to five-loop order (see Refs. [37-43] and references therein). $ F^{\rm{cut}} $ can easily be dealt with because the only nontrivial information is the cut propagators, which will be explicitly studied in Appendix A. With this information in hand, the next question is how to obtain physical integrals.

      It was shown in Ref. [44] that, for any given problem, Feynman loop integrals form a finite-dimensional vector space, the basis of which is master integrals (MIs). The step to express all loop integrals as linear combinations of MIs is called reduction. With the reverse unitarity relation in Eq. (1), one can map phase-space integrations onto corresponding loop integrations. Therefore, integrals over the phase-space and loop momenta defined in Eq. (2) can also be reduced to corresponding MIs.

      Reduction of a general integral to MIs can traditionally be achieved by using IBP relations based on Laporta's algorithm [45-50]. Alternatively, one can achieve IBP reduction using finite-field interpolation [51-56], module intersection [57], intersection theory [58], or AME [29]. 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 DEs among MIs [23, 24, 59-62]. Especially, DEs w.r.t. the AM $ \eta $ are given by

      $ \frac{\partial}{\partial \eta} \vec{J}(\vec s, \eta) = M(\vec s, \eta) \vec{J}(\vec s, \eta), $

      (10)

      where $ \vec{J}(\vec s, \eta)\equiv \left\{ {F}(\vec{\nu}', \vec s, \eta), {F}(\vec{\nu}'', \vec s, \eta), \cdots\right\} $ is a complete set of MIs and $ M(\vec s, \eta) $ is the coefficient matrix as a rational function of $ \vec s $ and $ \eta $. The boundary condition of the DEs can be chosen at $ \eta\to\infty $, which can easily be 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} $.

      We emphasize that, in practice, one should first reduce the scattering amplitude to MIs without AM, and then calculate these MIs by introducing AM. Therefore, the AMF method does not change the step of reducing scattering amplitude. It does make the step of setting up DEs w.r.t. AM $ \eta $ more complicated than the traditional DE method because one more scale is involved. However, the cost is always tolerable as it is still simpler than the step of reducing full scattering amplitude.

    III.   EXAMPLES: MASTER INTEGRALS FOR $ e^+e^-\rightarrow \gamma^* \rightarrow t\bar{t}+X $ AT NNLO
    • As a simple but nontrivial example, we calculate the MIs encountered in the NNLO correction for $ t\bar{t} $ production in an $ e^+ e^- $ collision mediated by a virtual photon to demonstrate the validity of the 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

      $ \hat{F}(\vec{\nu};x,y)\equiv s^{N-\frac{N-1+L}{2}D+\nu} {F}(\vec{\nu};\vec{s},\eta), $

      (11)

      where $ s = Q^2 $, $x = \dfrac{4m_t^2}{s}$, $y = \dfrac{\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:

      $ \begin{aligned}[b] &{\eta}^t_1 = {\eta}^{ \rm{t}}_2 = \cdots, \\ &{\eta}^+_1 = {\eta}^+_2 = \cdots = {\eta}^-_1 = {\eta}^-_2 = \cdots, \end{aligned} $

      (12)

      with $ \eta_1^{ \rm{t}}\ll \eta_1^+ $. More precisely, if $ \{{\cal{D}}^-_{\alpha}\} $ or $ \{{\cal{D}}^+_{\alpha}\} $ depends on $ \vec {s\,} $, we choose $ \eta^{ \rm{t}}_1 = 0^+ $ and $ \eta^+_1 = \eta $; otherwise, we choose $ \eta^{ \rm{t}}_1 = \eta $ (the introduction of $ \eta_1^+ $ is unnecessary in this case). The publicly available systematic package FIRE6 [47] is sufficient to do all needed reductions in this simple exercise.

      There are many subprocesses at the partonic level. We use $ {\rm{V}}^{L^+}{\rm{R}}^{N-1}{\rm{V}}^{L^-} $ to distinguish processes with different numbers of independent loop momenta and cut momenta. The presence of $ N-1 $ instead 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. The calculation of MIs for 4-particle cuts at NNLO (RRR) is presented in Sec. III B. The calculation of MIs for 3-particle cuts at NNLO (VRR) is presented in Sec. III C, and the calculation of MIs for 2-particle cuts at NNLO (VVR or VRV) is presented in Sec. III D. Verification of these results is given in Sec. III E. We note that all these results have already been calculated in the literature using other methods (see, e.g., Refs. [63-66]) although they are not fully publicly available. We provide our results as ancillary file with the arXiv version of this paper.

    • A.   RR and VR

    • For RR, i.e. the $ \gamma^*\to t\bar{t}g $ process, the inverse propagators can be written as

      $ \begin{aligned}[b] &{\cal{D}}^{\rm{c}}_1 = k_1^2-m_t^2,~ {\cal{D}}^{\rm{c}}_2 = k_2^2-m_t^2,~ {\cal{D}}^{\rm{c}}_3 = (Q-k_1-k_2)^2; \\ &{\cal{D}}^ {\rm{t}}_1 = (Q-k_1)^2-m_t^2,~ {\cal{D}}^ {\rm{t}}_2 = (Q-k_2)^2-m_t^2. \end{aligned}$

      (13)

      The 2 physical MIs obtained are $ F^{\rm{cut}}_{2,3,1} $ and $ F^{\rm{cut}}_{2,3,2} $, which are calculated in Appendix A.

      For VR, the inverse propagators can be written as

      $ \begin{aligned}[b] &{\cal{D}}^{\rm{c}}_1 = k_1^2-m_t^2,\quad {\cal{D}}^{\rm{c}}_2 = (Q-k_1)^2-m_t^2; \\ &{\cal{D}}^+_1 = (k_1+l^+_1)^2-m_t^2,\\& {\cal{D}}^+_2 = (Q-k_1-l^+_1)^2-m_t^2, {\cal{D}}^+_3 = {l^+_1}^2. \end{aligned}$

      (14)

      We obtain 1 family with $ \vec \nu $ of the characterized MIs being $ \{1,1,0\} $. Physical MIs are given by

      $ \begin{aligned}[b] &\int{\rm{dPS}}_2\int\frac{ {\rm{d}}^Dl^+_1}{(2\pi)^D}\frac{1}{{\cal{D}}^+_1+ {\rm{i}}0^+} = \Phi_2\int\frac{ {\rm{d}}^Dl^+_1}{(2\pi)^D}\frac{1}{{l^+_1}^2-m_t^2+ {\rm{i}}0^+},\,\\ &\int{\rm{dPS}}_2\int\frac{ {\rm{d}}^Dl^+_1}{(2\pi)^D}\frac{1}{({\cal{D}}^+_1+ {\rm{i}}0^+)({\cal{D}}^+_2+ {\rm{i}}0^+)} =\Phi_2\int\frac{ {\rm{d}}^Dl^+_1}{(2\pi)^D}\frac{1}{({l^+_1}^2-m_t^2+ {\rm{i}}0^+)((l^+_1-Q)^2-m_t^2+ {\rm{i}}0^+)}, \end{aligned} $

      (15)

      where $ \Phi_2 = \int{\rm{dPS}}_2 $ is defined in Eq. (A5) and the remaining 1-loop integrals are easy to calculate analytically.

    • B.   RRR

    • For RRR, i.e., $ \gamma^*\to t\bar{t}gg $ or $ t\bar{t} q\bar{q} $ processes, the inverse propagators can be written as

      $ \begin{aligned}[b] {\cal{D}}^{\rm{c}}_1 = &k_1^2-m_t^2, \quad{\cal{D}}^{\rm{c}}_2 = {k_2}^2-m_t^2,\\ {\cal{D}}^{\rm{c}}_3 =&k_3^2,\quad {\cal{D}}^{\rm{c}}_4 = (Q-k_1-k_2-k_3)^2;\,\\ {\cal{D}}^{{\rm{t}}}_1 = &(Q-k_1-k_2)^2,\quad {\cal{D}}^{{\rm{t}}}_2 = (k_1+k_3)^2-m_t^2,\\ {\cal{D}}^{{\rm{t}}}_3 =&(Q-k_2-k_3)^2-m_t^2, \quad{\cal{D}}^{{\rm{t}}}_4 = (Q-k_2)^2-m_t^2,\,\\ {\cal{D}}^{{\rm{t}}}_5 = &(k_2+k_3)^2-m_t^2, \quad{\cal{D}}^{{\rm{t}}}_6 = (Q-k_1-k_3)^2-m_t^2, \\ {\cal{D}}^{{\rm{t}}}_7 =& (Q-k_1)^2-m_t^2, \end{aligned} $

      (16)

      where $ Q^2 = s $. The phase-space integrals can then be expressed as

      $ \hat{F}(\vec{\nu};x,y) = s^{4-\frac{3}{2}D+\sum \nu^{{\rm{t}}}_{\alpha}} \int{\rm{dPS}}_4\prod\limits_{\alpha = 1}^{7}\frac{1}{({\cal{D}}^{{\rm{t}}}_\alpha+\eta)^{\nu^{{\rm{t}}}_{\alpha}}}. $

      (17)

      Using FIRE6, we find that there are 37 MIs for finite $ \eta $ and the number is reduced to 15 as $ \eta $ vanishing. These MIs can be classified into 2 families, with $ \vec \nu $ of the characterized MIs being

      $ \{0,1,1,0,1,1,0\} \; \; \; {\rm{and}}\; \; \; \{0,1,0,1,1,0,1\}, $

      (18)

      which can also characterized by Feynman diagrams Fig. 2(a) and (b), respectively. Using our method, we can calculate all physical MIs with any fixed x. The result of the most complicated MI with e.g. $ x = 1/2 $ is

      Figure 2.  Representative Feynman diagrams for $ \gamma^*\rightarrow t\bar{t}gg $ or $ t\bar{t} q\bar{q} $ processes, where (a) and (b) define the two most complicated families and (c) defines a sub-family of (b). The thick curves represent top quarks, the thin curves represent massless particles, and the vertical dashed lines represent the final state cut.

      $ \begin{aligned}[b] \hat F\left( \left\{-1, 1, 0, 1, 1, 0, 1\right\};\dfrac{1}{2},0\right) =&4.54087957883195468901389004370\times 10^{-7} +0.0000105911293014536670979999400899 \epsilon \\ & + 0.000124406630344529071923178953167 \epsilon^2 + 0.000983063887963543479220271851806 \epsilon^3 \\ &+ 0.00589205324016960475844728350032 \epsilon^4 + 0.0286458793127880349046701435743 \epsilon^5 \\ &+ 0.118038608602644851031978155328 \epsilon^6 + 0.425508191586298756000765241113 \epsilon^7 \\ &+ 1.37520023939856884048640232783 \epsilon^8 +\cdots, \end{aligned} $

      (19)

      where we truncate the expansion to order $ \epsilon^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 for MIs. For brevity, we define the MIs for this family as:

      $ \hat F(\{\nu^{ \rm{t}}_1,\nu^{ \rm{t}}_2,\nu^{ \rm{t}}_7\};x,y) \!=\! s^{4-\frac{3}{2}D+\nu^{ \rm{t}}_1+\nu^{ \rm{t}}_2+\nu^{ \rm{t}}_7} \!\!\!\!\int\!\!\!\!{\rm{dPS}}_4\frac{({\cal{D}}^{{\rm{t}}}_1\!+\!\eta)^{-\nu^{ \rm{t}}_1}({\cal{D}}^{{\rm{t}}}_2\!+\!\eta)^{-\nu^{ \rm{t}}_2}}{({\cal{D}}^{{\rm{t}}}_7+\eta)^{\nu^{ \rm{t}}_7}}, $

      (20)

      with $ \nu^{ \rm{t}}_1, \nu^{ \rm{t}}_2\leqslant 0 $ and $ \nu^{ \rm{t}}_7\geqslant 0 $. This family contains 5 MIs for finite $ \eta $ (or y)

      $ \begin{aligned}[b] & \{\hat F(\{0,0,0\};x,y),\hat F(\{-1,0,0\};x,y),\hat F(\{0,-1,0\};x,y),\\ &\hat F(\{0,0,1\};x,y),\hat F(\{-1,0,1\};x,y)\}, \end{aligned} $

      (21)

      and 4 MIs as $ \eta\to0^+ $

      $ \begin{aligned}[b] &\Big\{\hat F(\{0,0,0\};x,0),\hat F(\{-1,0,0\};x,0),\\ &\hat F(\{0,-1,0\};x,0),\hat F(\{0,0,1\};x,0) \Big\}. \end{aligned}$

      (22)

      It is clear that the first three MIs of the two sets are just linear combinations of $ F^{\rm{cut}} $ that have been calculated in Appendix A.

      To calculate the last physical MI in Eq. (22), we set up DEs for corresponding MIs with $ x = 1/2 $ and finite $ \eta $, which gives

      $\begin{aligned}[b] \frac{\partial}{\partial{y}} \left( \begin{array}{c} \hat F\left(\{0,0,1\};\dfrac{1}{2},y\right)\\ \hat F\left(\{-1,0,1\};\dfrac{1}{2},y\right) \end{array} \right)=&\left( \begin{array}{cc} \dfrac{ \binom{3 - 4 \epsilon + 14y - 12\epsilon y -8 y^2 }{+52 \epsilon y^2 - 16 y^3 + 48 \epsilon y^3}}{y (1- 8 y) (1 + 4 y + 2 y^2)}& \dfrac{ 12 (-1 + \epsilon)}{y (1- 8 y)} \\ \dfrac{-4(-1 + \epsilon)}{ 1 - 8 y }& \dfrac{ 2 \binom{-7 + 6\epsilon - 39 y +46\epsilon y}{ - 24 y^2 + 32 \epsilon y^2}}{(1- 8 y) (1 + 4y + 2 y^2) } \end{array} \right)\left(\begin{array}{c} \hat F\left(\{0,0,1\};\dfrac{1}{2},y\right)\\ \hat F\left(\{-1,0,1\};\dfrac{1}{2},y\right) \end{array} \right)\\ &+\left( \begin{array}{c} \dfrac{-2 \binom{15- 16\epsilon + 53 y -46 \epsilon y}{ - 124 y^2 + 132 \epsilon y^2}}{y (1 - 8 y) (1 + 4y + 2 y^2)}\quad \dfrac{ - 48 (-1 + \epsilon) (1 - y)}{y (1 - 8 y) (1 + 4 y + 2 y^2) }\quad \dfrac{ -144 (-1 + \epsilon) (1 - y)}{y (1 - 8 y) (1 + 4 y + 2 y^2) } \\ \dfrac{- 2 \binom{17 - 18 \epsilon + 85y }{ -78\epsilon y - 88 y^2 + 96 \epsilon y^2}}{(1- 8 y) (1 + 4 y + 2 y^2) }\quad \dfrac{-56(-1 + \epsilon)}{(1 - 8 y) (1 + 4 y + 2 y^2) }\quad \dfrac{ -8 (-1 + \epsilon) (23- 16 y)}{(1 - 8 y) (1 + 4y + 2 y^2) } \end{array} \right)\\ &\times \left( \begin{array}{c} \hat F\left(\{0,0,0\};\dfrac{1}{2},y\right)\\ \hat F\left(\{-1,0,0\};\dfrac{1}{2},y\right)\\ \hat F\left(\{0,-1,0\};\dfrac{1}{2},y\right) \end{array} \right). \end{aligned} $

      (23)

      The boundary condition of these DEs is given by

      $ \begin{aligned}[b] & \hat F\left(\{0,0,1\};\frac{1}{2},y\right) {\xlongequal[\;] {\eta \to \infty }}s^{5-\frac{3}{2}D}\left(\frac{ 1}{\eta}\int{\rm{dPS}}_4-\frac{1}{\eta^2} \int{\rm{dPS}}_4{\cal{D}}^{{\rm{t}}}_7\right)_{x = 1/2},\\ &\hat F\left(\{-1,0,1\};\frac{1}{2},y\right) {\xlongequal[\;] {\eta \to \infty }}s^{4-\frac{3}{2}D}\left( \int{\rm{dPS}}_4-\frac{1}{\eta} \int{\rm{dPS}}_4({\cal{D}}^{{\rm{t}}}_1-{\cal{D}}^{{\rm{t}}}_7)\right)_{x = 1/2}. \end{aligned} $

      (24)

      By solving the DEs (23) with the boundary condition we obtain, e.g.,

      $ \begin{aligned}[b] \hat F\left(\{0,0,1\};\frac{1}{2},0\right) =& 1.08703304446867684362962983900\times10^{-8} + 2.35523252532576916745398951290\times10^{-7} \epsilon \\ & + 2.54552584682491818967669600707\times10^{-6} \epsilon^2 + 0.0000183017607702920607720306203025 \epsilon^3 \\ & + 0.0000985014117985155673601315129686 \epsilon^4 + 0.000423451568312682839128853965245 \epsilon^5 \\ & + 0.00151532069472999876700566927592 \epsilon^6 + 0.00464561935910925265889327482252 \epsilon^7 \\ & + 0.0124659262732530489584879132787 \epsilon^8 +\cdots. \end{aligned} $

      (25)

        Finally, the DE w.r.t. x is given by

      $ \begin{aligned}[b] \frac{\partial}{\partial{x}} \hat F(\{0,0,1\};x,0) = & \frac{-2 + 3 x + 4 \epsilon - 6 x \epsilon}{2 (-1 + x) x} \hat F(\{0,0,1\};x,0)+ \frac{-11 + 2 x + 12 \epsilon - 2 x \epsilon}{(-1 + x)^2 x}\hat F(\{0,0,0\};x,0)\\ &+\frac{ - 12 (-1 + \epsilon)}{(-1 + x)^2 x}\hat F(\{-1,0,0\};x,0)+ \frac {- 36 (-1 + \epsilon)}{(-1 + x)^2 x}\hat F(\{0,-1,0\};x,0). \end{aligned} $

      (26)

      With the boundary condition at $ x = 1/2 $ given in Eq. (25), solving the above DE can also give $ \hat F(\{0,0,1\};x,0) $ at any value of x.

    • C.   VRR

    • For VRR, the inverse propagators can be written as

      $ \begin{aligned}[b] {\cal{D}}^{\rm{c}}_1 =& k_1^2-m_t^2, \\{\cal{D}}^{\rm{c}}_2 =& k_2^2-m_t^2,\\{\cal{D}}^{\rm{c}}_3 =& (Q-k_1-k_2)^2;\,\\ {\cal{D}}^{ \rm{t}}_1 =& (Q-k_2)^2 - m_t^2, \\{\cal{D}}^{ \rm{t}}_2 =& ( Q-k_1 )^2- m_t^2;\,\\ {\cal{D}}^+_1 =& (k_1 + l^+_1)^2 - m_t^2,\\{\cal{D}}^+_2 =& (k_2 - l^+_1 )^2-m_t^2 ,\\ {\cal{D}}^+_3 =& {l^+_1}^2,\\ {\cal{D}}^+_4 =& (Q-k_1-k_2 + l^+_1)^2,\,\\ {\cal{D}}^+_5 =& (Q-k_2 + l^+_1)^2 - m_t^2, \\{\cal{D}}^+_6 =& ( Q-k_1 - l^+_1 )^2-m_t^2 , \end{aligned} $

      (27)

      where $ Q^2 = s $. Then the 1-loop $ 1\to3 $ phase-space integrals can be expressed as

      $\begin{aligned}[b] \hat{F}(\vec{\nu};x,y) =& s^{3-\frac{3}{2}D+\sum \nu^{{\rm{t}}}_{\alpha}+\sum \nu^+_{\alpha}} \int{\rm{dPS}}_3\prod\limits_{\alpha = 1}^{2}\frac{1}{({\cal{D}}^{{\rm{t}}}_{\alpha})^{\nu^{{\rm{t}}}_{\alpha}}}\\&\times \int\frac{ {\rm{d}}^Dl^+_1}{(2\pi)^D}\prod\limits_{\alpha = 1}^{6}\frac{1}{({\cal{D}}^+_{\alpha}+ {\rm{i}}\eta)^{\nu^+_{\alpha}}}. \end{aligned}$

      (28)

      We find that there are 71 MIs for finite $ \eta $ (or y) and the number is reduced to 27 when $ \eta\to 0^+ $. These MIs can be classified into 2 families,

      $ \{ 1, 0, 1, 1, 1, 0, 0, 1\}\; \; \; {\rm{and}}\; \; \; \{ 1, 0, 0, 1, 0, 1, 1, 1\}. $

      (29)

      The first family is the most complicated. It can also be characterized by the Feynman diagram in Fig. 3(a). Using our method, the result of the most complicated MI with, e.g., $ x = 1/2 $ is

      Figure 3.  Representative Feynman diagrams in VRR, where (a) defines the most complicated family and (b) defines a sub-family of (a). The thick curves represent top quarks, the thin curves represent massless particles, and the vertical dashed lines represent the final state cut.

      $ \begin{aligned}[b] \hat{F}\left(\{ 2, 0, 1, 1, 1, 0, 0, 1\};\frac{1}{2},0\right) =& ( 4.35941187166437229714484148598\times10^{-6} - 4.87955595721663859057448350469\times10^{-6} {\rm{i}})\epsilon^{-2} \\ & +( 0.000290878052291807102955726741096 + 0.000037713113716691251718223752983 {\rm{i}})\epsilon^{-1} \\ &+(0.00232637225490549068317799097260 + 0.00161913500108049877728334443544 {\rm{i}}) \\ & +(0.0064623807207395294287699841138 + 0.0143681169585383071360453670293 {\rm{i}}) \epsilon \\ & -(0.0244064366687345660260807481505 - 0.0681410637818588674321247511745 {\rm{i}}) \epsilon^2 \\ & -(0.324062211265101500673336745067 - 0.176280502471673558710833935083 {\rm{i}}) \epsilon^3 \\ & -(1.84537524457279855776436524767 - 0.05103632828295095891823713650 {\rm{i}}) \epsilon^4 \\ & -(7.60322033887241114962526507189 + 2.00523602305209729225222600674 {\rm{i}}) \epsilon^5 \\ & -(27.2739377963678526029153832666 + 12.0357468364668313721542998217 {\rm{i}}) \epsilon^6+\cdots. \end{aligned} $

      (30)

      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 for MIs. For brevity, we define the MIs for this family as:

      $ \hat F(\{\nu^{ \rm{t}}_1,\nu^+_1,\nu^+_2\};x,y) =s^{3-\frac{3}{2}D+\nu^{ \rm{t}}_1+\nu^+_1+\nu^+_2} \int{\rm{dPS}}_3\frac{1}{({\cal{D}}^{{\rm{t}}}_1)^{\nu^{{\rm{t}}}_1}} \\ \int\frac{ {\rm{d}}^Dl^+_1}{(2\pi)^D}\frac{1}{({\cal{D}}^+_{1}+ {\rm{i}}\eta)^{\nu^+_1}({\cal{D}}^+_{2}+ {\rm{i}}\eta)^{\nu^+_2}}. $

      (31)

      This family contains 7 MIs for finite $ \eta $

      $ \{ \hat F(\{0,0,1\};x,y),\hat F(\{-1,0,1\};x,y),\hat F(\{0,1,1\};x,y),\hat F(\{-1,1,1\};x,y),\hat F(\{0,1,2\};x,y),\hat F(\{1,1,1\};x,y),\hat F(\{1,1,2\};x,y) \}, $

      (32)

      and 6 MIs as $ \eta\to0^+ $

      $ \{\hat F(\{0,0,1\};x,0),\hat F(\{-1,0,1\};x,0),\hat F(\{0,1,1\};x,0),\hat F(\{0,1,2\};x,0), \hat F(\{1,1,1\};x,0),\hat F(\{1,1,2\};x,0)\}. $

      (33)

        To calculate the last 2 physical MIs, we set up DEs for corresponding MIs with $ x = 1/2 $ and finite $ \eta $, which gives

      $ \begin{aligned}[b] \frac{\partial}{\partial{y}}\hat F\left(\{1,1,1\};\frac{1}{2},y\right) = & -2 {\rm{i}}\hat F\left(\{1,1,2\};\frac{1}{2},y\right),\\ \frac{\partial}{\partial{y}}\hat F\left(\{1,1,2\};\frac{1}{2},y\right) =& \frac{ 16 {\rm{i}} (1 - 2 \epsilon) ( 3 \epsilon-1)}{(- {\rm{i}} + 8 y) ( {\rm{i}} + 8 y)}\hat F\left(\{1,1,1\};\frac{1}{2},y\right)+ \frac{ 4 ( {\rm{i}} + 8 y - 4 {\rm{i}} \epsilon - 64 y \epsilon)}{(- {\rm{i}} + 8 y) ( {\rm{i}} + 8 y)}\hat F\left(\{1,1,2\};\frac{1}{2},y\right)\\ &-\frac{64 (5 - 11 \epsilon + 6 \epsilon^2)}{ y (- {\rm{i}} + 8 y) ( {\rm{i}} + 8 y)}\hat F\left(\{0,0,1\};\frac{1}{2},y\right)+ \frac{768 (-1 + \epsilon)^2}{ y (- {\rm{i}} + 8 y) ( {\rm{i}} + 8 y)}\hat F\left(\{-1,0,1\};\frac{1}{2},y\right)\\ &+\frac{ 8 (1 - 2 \epsilon) (5 - 4 {\rm{i}} y -6 \epsilon + 8 {\rm{i}} y \epsilon)}{ y (- {\rm{i}} + 8 y) ( {\rm{i}} + 8 y)}\hat F\left(\{0,1,1\};\frac{1}{2},y\right) +\frac{ 32 (1 - 2 \epsilon) (-3 + 4 \epsilon)}{ y (- {\rm{i}} + 8 y) ( {\rm{i}} + 8 y)}\hat F\left(\{-1,1,1\};\frac{1}{2},y\right)\\ &-\frac{4 (- {\rm{i}} + 4 y) (-1 + 2 \epsilon)}{ y (- {\rm{i}} + 8 y)}\hat F\left(\{0,1,2\};\frac{1}{2},y\right). \end{aligned} $

      (34)

      The boundary condition for these DEs is given by

      $ \begin{aligned}[b] &\hat F\left(\{1,1,1\};\frac{1}{2},y\right) {\xlongequal[\;] {\eta \sim \infty }} s^{6-\frac{3}{2}D} \eta^{\frac{D}{2}-2}\frac{ {\rm{i}}(D-2)}{2}F^{\rm{bub}}_{1,1}(D)\left(\int{\rm{dPS}}_3\frac{1}{{\cal{D}}^ {\rm{t}}_1}\right)_{x = 1/2},\\ &\hat F\left(\{1,1,2\};\frac{1}{2},y\right) {\xlongequal[\;] {\eta \sim \infty }} s^{7-\frac{3}{2}D} \eta^{\frac{D}{2}-3}\frac{(4 -D) (D-2)}{8}F^{\rm{bub}}_{1,1}(D)\left(\int{\rm{dPS}}_3\frac{1}{{\cal{D}}^ {\rm{t}}_1}\right)_{x = 1/2}, \end{aligned}$

      (35)

      where

      $ F^{\rm{bub}}_{1,1}(D) \equiv\int\frac{ {\rm{d}}^D{l^+_1}}{(2\pi)^D}\frac{1}{{l^+_1}^2+ {\rm{i}}}, $

      (36)

      and $\int{\rm{dPS}}_3\dfrac{1}{{\cal{D}}^ {\rm{t}}_1}$ can again be calculated by the method inSec. III B or obtained from RR in Sec. III A. Knowing the first 5 MIs, by solving the DEs (34) with the boundary condition we obtain e.g.,

      $ \begin{aligned}[b] \hat F\left(\{1,1,2\};\frac{1}{2},0\right) =& (7.78790446721069262502850093774\times10^{-6} + 2.91319469772237394135356308348\times10^{-6} {\rm{i}}) \\ & +(0.000130430373015787655604488198861 + 0.000068404169458201184291920092123 {\rm{i}}) \epsilon \\ & +(0.001077434813828191909787186362432 + 0.000750926876250745472210277589430 {\rm{i}}) \epsilon^2 \\ & +(0.00584278150920839062615612508136 + 0.00527570101382158661589031061691 {\rm{i}}) \epsilon^3 \\ & +(0.0233461280012444372334494219123 + 0.0270859736951617524563966282868 {\rm{i}}) \epsilon^4 \\ & +(0.0730918539437148667076104654800 + 0.1095165249743204252589933869672 {\rm{i}}) \epsilon^5 \\ & +(0.185975373883125986488613881520 + 0.366393770042708443331564801509 {\rm{i}}) \epsilon^6 \\ & +(0.393093986188519076512424694564 + 1.052172170765638116257825410632 {\rm{i}}) \epsilon^7 \\ & +(0.69775277299606861048706250047 + 2.67282546122008383022615104289 {\rm{i}}) \epsilon^8+\cdots. \end{aligned} $

      (37)

        Finally, the DEs w.r.t. x are given by

      $ \begin{aligned}[b] \frac{\partial}{\partial{x}} \left(\begin{array}{c} \hat F(\{1,1,1\};x,0)\\ \hat F(\{1,1,2\};x,0) \end{array}\right)=&\left(\begin{array}{c} \dfrac{ -\epsilon}{ x}\qquad\qquad\qquad \dfrac{ 1}{2}\\ \dfrac{(1 - 2 \epsilon) ( 3 \epsilon-1)}{(-1 + x) x}\quad \dfrac{x + 4 \epsilon - 10 x \epsilon}{2 (-1 + x) x} \end{array}\right) \left(\begin{array}{c} \hat F(\{1,1,1\};x,0)\\ \hat F(\{1,1,2\};x,0) \end{array}\right)\\ &+\left( \begin{array}{c} \dfrac{ 6 (-1 + \epsilon) (-4 - x + 4 \epsilon + 2 x \epsilon)}{(-1 + x)^2 x^2 (-1 + 2 \epsilon)}\quad \dfrac{ - 72 (-1 + \epsilon)^2}{(-1 + x)^2 x^2 (-1 + 2 \epsilon)}\quad \dfrac{ - 2 (-1 + 2 \epsilon)}{(-1 + x) x}\quad\dfrac{ -2}{-1 + x} \\ \dfrac{- 4 (4 + 3 x - 8 \epsilon - 8 x \epsilon + 4 \epsilon^2 + 5 x \epsilon^2)}{(-1 + x)^2 x^3}\quad \dfrac{ 48 (-1 + \epsilon)^2}{(-1 + x)^2 x^3}\quad \dfrac{ 2 (1-2\epsilon)^2}{(-1 + x)^2 x}\quad \dfrac{ 2 (-1 + 2 \epsilon)}{(-1 + x)^2} \end{array} \right)\\ &\times \left( \begin{array}{c} \hat F(\{0,0,1\};x,0)\\ \hat F(\{-1,0,1\};x,0)\\ \hat F(\{0,1,1\};x,0)\\ \hat F(\{0,1,2\};x,0) \end{array} \right). \end{aligned} $

      (38)

      With the boundary condition at $ x = 1/2 $ given in Eq. (37), by solving the above DEs we can also evaluate $ \hat F(\{1,1,1\};x,0) $ and $ \hat F(\{1,1,2\};x,0) $ at any value of x.

    • D.   VVR and VRV

    • For VVR, the inverse propagators can be written as

      $ \begin{aligned}[b] &{\cal{D}}^{\rm{c}}_1 = k_1^2-m_t^2, {\cal{D}}^{\rm{c}}_2 = (Q-k_1)^2-m_t^2;\, {\cal{D}}^+_1 = (k_1 + l^+_1 + l^+_2)^2 - m_t^2, {\cal{D}}^+_2 = (k_1 + l^+_2)^2 - m_t^2, {\cal{D}}^+_3 = (k_1 + l^+_1)^2 - m_t^2, \\ &{\cal{D}}^+_4 = (-k_1 - l^+_1 + Q)^2-m_t^2 ,\,{\cal{D}}^+_5 = {l^+_1}^2, {\cal{D}}^+_6 = {l^+_2}^2, {\cal{D}}^+_7 = (-k_1 - l^+_2 + Q)^2-m_t^2 , {\cal{D}}^+_8 = (-k_1 - l^+_1 - l^+_2 + Q)^2-m_t^2 , {\cal{D}}^+_9 = (l^+_1 + l^+_2)^2 . \end{aligned} $

      (39)

      Then the phase-space integrals can be expressed as

      $ \hat F(\vec{\nu};x,y) = s^{2-\frac{3}{2}D+ \sum \nu^+_{\alpha}} \int{\rm{dPS}}_2 \int\frac{ {\rm{d}}^Dl^+_1 {\rm{d}}^Dl^+_2}{(2\pi)^{2D}}\prod\limits_{\alpha = 1}^{9}\frac{1}{({\cal{D}}^+_{\alpha}+ {\rm{i}}\eta)^{\nu^+_{\alpha}}}. $

      (40)

      We find that there are 53 MIs for finite $ \eta $ and the number is reduced to 21 when $ \eta\to 0^+ $. These MIs can be classified into 3 families

      $ \begin{aligned}[b] &\{ 1, 1, 0, 1, 1, 1, 0, 1, 0\} \; \; \; {\rm{and}}\; \; \; \{ 1, 1, 1, 1, 0, 1, 0, 0,1\}\; \; \; {\rm{and}}\; \; \; \\ &\{1, 1, 1, 0, 0, 1, 1, 1, 0\}. \end{aligned} $

      (41)

      The first family is the most complicated. It can be characterized by the Feynman diagram in Fig. 4(a). We use the same method (note, $ F^{\rm{bub}} $ up to 2-loop is needed to calculate boundary conditions) in VRR to calculate MIs for these families.

      Figure 4.  Representative Feynman diagrams in VVR and VRV, where (a) defines the most complicated family for VVR and (b) defines the family for VRV. The thick curves represent top quarks, the thin curves represent massless particles, and the vertical dashed lines represent the final state cut.

      For VRV, the inverse propagators can be written as

      $ \begin{aligned}[b] &{\cal{D}}^{\rm{c}}_1 = k_1^2-m_t^2, {\cal{D}}^{\rm{c}}_2 = (Q-k_1)^2-m_t^2;\,\\ &{\cal{D}}^+_1 = (k_1 + l^+_1)^2 - m_t^2, {\cal{D}}^+_2 = (k_1 + l^+_1 - Q)^2 -m_t^2 , {\cal{D}}^+_3 = {l^+_1}^2;\,\\ &{\cal{D}}^-_1 = (k_1 + l^-_1)^2-m_t^2 , {\cal{D}}^-_2 = (k_1 + l^-_1- Q)^2-m_t^2 , {\cal{D}}^-_3 = {l^-_1}^2, \end{aligned} $

      (42)

      in addition to which there is a scalar product $ l_1^+\cdot l_1^- $. The obtained family is:

      $ \{1,1,0,1,1,0,0\}, $

      (43)

      which can also be characterized by the Feynman diagram Fig. 4(b). Because loop integrations in MIs of this case are factorized, their calculation is as simple as the 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 the results of MIs obtained by the two strategies for different values of x, e.g., $ x = 1/2 $ and $ x = 2/5 $, are consistent with each other. This provides a highly nontrivial self-consistency check because DEs w.r.t. $ \eta $ and w.r.t. x are significantly different.

      Our numerical values for $ \hat F(\vec{\nu};0,0) $, i.e. MIs for massless QCD, are in full agreement with the known analytical results of massless MIs in the literature [67, 68].

      For the RRR sub-process (Sec. III B), our numerical results for $ \hat F(\{0,0,0,0,0,0,1\};x,0) $ and $ \hat F(\{0,0,0,1,0,0,1\}; x,0) $ show excellent agreement with the corresponding analytical results in the literature [63] ($ T_4 $ and $ T_5 $ in this reference).

    IV.   SUMMARY AND OUTLOOK
    • In this paper, we have extended the AMF method originally developed for Feynman loop integrals [28] to calculate MIs which also involve phase-space integration. As a pedagogical example, we used this method to calculate the MIs encountered in $ e^+e^-\rightarrow \gamma^* \rightarrow t\bar{t}+X $ at NNLO. Our results agree with results obtained by using other methods (our full results are available as an ancillary file with the arXiv version of this paper). Although the AMF method depends on a reduction procedure to decompose all integrals to MIs and to set up the DEs of MIs w.r.t. $ \eta $, 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 systematically as the sector decomposition method. However, compared with the latter method, the AMF method is much more efficient and can provide very high precision. Its high-precision nature makes it possible to obtain analytical results with a proper ansatz. For problems with two or more kinematical invariants, where the DE method works, the AMF method can not only systematically provide as many boundary conditions as needed for DEs w.r.t. kinematical invariants, but also provides a highly nontrivial self-consistency check for the obtained results. All the above advantages make the AMF method very useful for perturbative calculation at high orders.

    ACKNOWLEDGMENTS
    • We thank Chen-Yu Wang, K.T. Chao and C. Meng, Xin Guan, Xiao-Bo Jin, Zhi-Feng Liu for many helpful communications and discussions.

    APPENDIX A: CALCULATION OF $ F^{\rm{cut}} $
    • In this appendix, we calculate the MIs of basal phase-space integration with no denominator. It is important to note that, besides the masses presented in cut lines, these MIs only depend on the square of the center of mass energy $ s\equiv Q^2 $, regardless of the configuration of external unintegrated momenta. We use $ F_{r,N,n}^{\rm{cut}} $ to denote the n-th MI for N-particle-cut integrals with $ m_1 = \cdots = m_r = m $ and $ m_{r+1} = \cdots = m_N = 0 $, and we provide explicit results for $ r = 0,\; 1 $, and $ 2 $. MIs for general cases can be studied as follows. By using the optical theorem (see e.g. Ref. [67]), the calculation of the basal phase-space integral is translated to the calculation of the imaginary part of the corresponding sunrise pure loop integral. The latter 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. [69].

      For $ r = 0 $ and $N\geqslant 2$, there is only one MI, which is given by

      $ F_{0,N,1}^{\rm{cut}}\!\!\equiv\!\!\int\!\!{\rm{dPS}}_N \!=\! \frac{2^{5-4N-2\epsilon+2N\epsilon}\pi^{3-2N-\epsilon+N\epsilon} \Gamma(1-\epsilon)^N}{\Gamma((N-1)(1-\epsilon))\Gamma(N(1-\epsilon))}s^{N-2+\epsilon-N\epsilon}. \tag{A1} $

      (A1)

      For $ r = 1 $, in general there are two MIs, which are given by ($ n = 1,2 $)

      $ \begin{aligned}[b] F^{\rm{cut}}_{1,N,n}\equiv & \int{\rm{dPS}}_N\left((Q-k_1)^2\right)^{n-1} = \frac{2^{2(5-n-3N-2\epsilon+2N\epsilon)}\pi^{\frac{7}{2}+N(\epsilon-2)-\epsilon}\Gamma(1-\epsilon)^{N-1}\Gamma(n-3+N+2\epsilon-N\epsilon)}{ \Gamma((N-1)(1-\epsilon))\Gamma\left(n-\dfrac{3}{2}+N+\epsilon-N\epsilon\right)\Gamma((N-2)(1-\epsilon))}\\ &\times s^{n-3+N+\epsilon-N\epsilon}\left(1-\frac{m^2}{s}\right)^{2(n+N+\epsilon-N\epsilon)-5}\\ &\times {}_2F_1\left(n-2+N+\epsilon-N\epsilon,n-3+N+2\epsilon-N\epsilon;2(n-2+N+\epsilon-N\epsilon);1-\frac{m^2}{s}\right), \end{aligned} \tag{A2}$

      (A2)

      where $ {}_pF_q $ are generalized hypergeometric functions which can be evaluated using the publicly available program HypExp [70]. In the special case of $ N = 2 $, there is only one MI $ F^{\rm{cut}}_{1,2,1} $. For $ r = 2 $, in general there are three MIs, which are given by ($ n = 1,2,3 $)

      $ \begin{aligned}[b] F^{\rm{cut}}_{2,N,n}\equiv & \int{\rm{dPS}}_N\left((k_1+k_2)^2\right)^{n-1} = \frac{2^{5+2N(\epsilon-2)-2\epsilon}\pi^{3+N(\epsilon-2)-\epsilon}\Gamma(1+n-2\epsilon)\Gamma(1-\epsilon)^{N-1}\Gamma(n-\epsilon)}{\Gamma(2-2\epsilon) \Gamma(n-1+N-N\epsilon)\Gamma(n-2+N+\epsilon-N\epsilon)} s^{n-3+N+\epsilon-N\epsilon}\\ & \times {}_3F_2\left(\epsilon-\dfrac{1}{2},2-n-N+N\epsilon,3-n-N-\epsilon+N\epsilon;1-n+\epsilon,2\epsilon-n;\frac{4m^2}{s}\right)\\ &+\frac{2^{4+2N(\epsilon-2)}\pi^{\frac{7}{2}+N(\epsilon-2)-\epsilon}\Gamma(1-\epsilon)^{N-1}\Gamma(\epsilon-n)}{\Gamma\left(\dfrac{3}{2}-n\right) \Gamma((N-1)(1-\epsilon))\Gamma((N-2)(1-\epsilon))} s^{n-3+N+\epsilon-N\epsilon}\left(\frac{4m^2}{s}\right)^{n-\epsilon}\\ &\times {}_3F_2\left(n-\dfrac{1}{2},3-N-2\epsilon+N\epsilon,2-N-\epsilon+N\epsilon;1+n-\epsilon,\epsilon;\frac{4m^2}{s}\right)\\ &+\frac{2^{4+2N(\epsilon-2)}\pi^{\frac{7}{2}+N(\epsilon-2)-\epsilon}\Gamma(1-\epsilon)^{N-2}\Gamma(\epsilon-1)\Gamma(2\epsilon-1-n)}{\Gamma\left(\dfrac{1}{2}-n+\epsilon\right) \Gamma((N-2)(1-\epsilon))\Gamma((N-3)(1-\epsilon))} s^{n-3+N+\epsilon-N\epsilon}\left(\frac{4m^2}{s}\right)^{1+n-2\epsilon}\\ &\times {}_3F_2\left(\dfrac{1}{2}+n-\epsilon,4-N-3\epsilon+N\epsilon,3-N-2\epsilon+N\epsilon;2+n-2\epsilon,2-\epsilon;\frac{4m^2}{s}\right). \end{aligned}\tag{A3}$

      (46)

      In the special case of $ N = 2 $, there is only one MI, which means $ F^{\rm{cut}}_{2,2,2} $ and $ F^{\rm{cut}}_{2,2,3} $ can be reduced to $ F^{\rm{cut}}_{2,2,1} $. In the special case of $ N = 3 $, there are only two MIs, which means $ F^{\rm{cut}}_{2,3,3} $ can be reduced to $ F^{\rm{cut}}_{2,3,1} $ and $ F^{\rm{cut}}_{2,3,2} $. The number of MIs obtained here is consistent with the expectation in Ref. [71]. The results (A2) and (A3) agree with the partial results in the literature (see Eq. (B.7) in Ref. [72] and Eqs. (3.14-3.16) in Ref. [63]).

      To obtain Eqs. (A2) and (A3), we decompose the measure of phase-space integration into two parts:

      $ \begin{aligned}[b]&{\rm{dPS}}_N(Q;\{k_i\})_{\overline {\;s = {Q^2}} }^{\underline {{s_A} = k_A^2} }\int_{(\sum\nolimits_{i\in\Omega}m_i)^2}^{(\sqrt{s}-\sum\nolimits_{i\notin\Omega}m_i)^2}\\&\frac{ {\rm{d}} s_A}{2\pi}{\rm{dPS}}_{N-\#_{\Omega}+1}(Q; \left\{k_{i\notin\Omega}\right\}, k_A){\rm{dPS}}_{\#_{\Omega}}(k_A; \left\{k_{i\in\Omega}\right\}), \end{aligned} \tag{A4} $

      (A4)

      where $ \Omega $ is a subset of $ \{1,2,\cdots,N\} $ and $ \#_{\Omega} $ is the number of elements of $ \Omega $. Specifically, $ \Omega $ 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 the two-particle cut is used:

      $\begin{aligned}[b] \Phi_2 =& \int{\rm{dPS}}_2 = \frac{2^{-3 + 2\epsilon}\pi^{-1 + \epsilon}\Gamma(1 - \epsilon)}{\Gamma(2 -2\epsilon)}s^{-\epsilon}\\&\times\left( 1 - \frac {2 (m_ 1^2 + m_ 2^2)} {s} + \left(\frac {m_ 1^2 - m_2^2}{s}\right)^2\right)^{\frac{1}{2}-\epsilon}.\end{aligned} \tag{A5} $

      (A5)

      The following relations (which can be found from e.g. Ref. [73]) are also useful:

      $ {}_2F_1\left(a,b;c;z\right) = (1-z)^{-a}{}_2F_1\left(a,c-b;c;\dfrac{z}{z-1}\right), \tag{A6} $

      (A6)

      $ {}_2F_1\left(a,1-a;c;z\right) = (1-z)^{c-1}\left(\sqrt{1-z}+\sqrt{-z}\right)^{2-2a-2c} {}_2F_1\left(a+c-1,c-\frac{1}{2};2c-1;\frac{4\sqrt{z^2-z}}{(\sqrt{1-z}+\sqrt{-z})^2}\right),\tag{A7} $

      (A7)

      $ \begin{aligned}[b] \int_{z}^{y} {\rm{d}}{x}x^{\alpha-1}(y-x)^{c-1}(x-z)^{\beta-1}{}_2F_1\left(a,b;c;1-\frac{x}{y}\right) = & \frac{\Gamma(c)\Gamma(c-a-b)\Gamma(\beta)\Gamma(1-\alpha-\beta)}{\Gamma(c-a)\Gamma(c-b)\Gamma(1-\alpha)}y^{c-1}z^{\alpha+\beta-1}\\ &\times {}_3F_2\left(a-c+1,b-c+1,\alpha;a+b-c+1,\alpha+\beta;\frac{z}{y}\right) \end{aligned} $

      $ \begin{aligned}[b] &+ \frac{\Gamma(c)\Gamma(a+b-c)\Gamma(\beta)\Gamma(a+b-c-\alpha-\beta+1)}{\Gamma(a)\Gamma(b)\Gamma(a+b-c-\alpha+1)}y^{a+b-1}z^{c-a-b+\alpha+\beta-1} \times {}_3F_2\bigg(1-a,1-b,c-a-b+\alpha;c-a-b+1,\\ &\quad c-a-b+\alpha+\beta;\frac{z}{y}\bigg)+ \frac{\Gamma(c)\Gamma(\alpha+\beta-1)\Gamma(c-a-b+\alpha+\beta-1)}{\Gamma(c-a+\alpha+\beta-1)\Gamma(c-b+\alpha+\beta-1)}y^{c+\alpha+\beta-2}\\ &\times {}_3F_2\left(1-\beta,a-c-\alpha-\beta+2,b-c-\alpha-\beta+2; 2-\alpha-\beta,a+b-c-\alpha-\beta+2;\frac{z}{y}\right). \end{aligned} \tag{A8} $

      (A8)
Reference (73)

目录

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return