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.
  • 加载中
  • [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. B485 (1997) 291–419[hep-ph/9605323][InSPIRE]. [Erratum: Nucl. Phys.B510, 503(1998)].
    [7] V. Del Duca, C. Duhr, G. Somogyi, F. Tramontano, and Z. Trócsányi, JHEP 04, 036 (2015), arXiv:1501.07226
    [8] B. W. Harris and J. F. Owens, Phys. Rev. D65, 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, S. Jones, M. Kerner, J. Schlenk, and T. Zirke, 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, G. P. Salam, and G. Zanderighi, Fully Differential Vector-Boson-Fusion Higgs Production at Next-to-Next-to-Leading Order, Phys. Rev. Lett. 115 (2015) 082002[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, and F. Petriello, Phys. Rev. Lett. 115, 062002 (2015), arXiv:1504.02131 doi: 10.1103/PhysRevLett.115.062002
    [15] J. Gaunt, M. Stahlhofen, F. J. Tackmann, and J. R. Walsh, 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, C. Signorile-Signorile, P. Torrielli, and S. Uccirati, Local analytic sector subtraction at NNLO, JHEP 12 (2018) 107[arXiv: 1806.09570][InSPIRE]. [Erratum: JHEP 06, 013 (2019)].
    [18] L. Magnea, E. Maina, G. Pelliccioli, C. Signorile-Signorile, P. Torrielli, and S. Uccirati, JHEP 12, 062 (2018), arXiv:1809.05444
    [19] C. Anastasiou and K. Melnikov, Nucl. Phys. B646, 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, and F. Petriello, 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. B779, 353-357 (2018), arXiv:1711.09572
    [29] X. Liu and Y.-Q. Ma, Phys. Rev. D99, 071501 (2019), arXiv:1801.10523
    [30] P. Zhang, C.-Y. Wang, X. Liu, Y.-Q. Ma, C. Meng, and K.-T. Chao, 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, S.-M. Wang, W.-G. Ma, J.-W. Zhu, and Y. Jiang, 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, and Y. Schroder, 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, and J. Usovitsch, 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, H. Schönemann, and Y. Zhang, 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 (1993) 299–308 [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. B682, 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, and B. Mistlberger, 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.
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. shu
Milestone
Received: 2020-09-24
Article Metric

Article Views(52)
PDF Downloads(18)
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 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.

    HTML

    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.

      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 [1-18]. If the process in 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 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 [22], differential equations (DEs) [23, 24], dimensional recurrence relations [25-27], and also methods developed by introducing auxiliary mass (AM) [28-34]. 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^+\geqslant 0 $, $ L^-\geqslant 0 $, $ M\geqslant 1 $ and $ N\geqslant 1 $ are understandable. We denote Q as the total cut momentum which 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.  A 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 inverse of tree propagators which do not depend on loop momenta, $ {\cal{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, $ {\cal{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}^{ \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 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 rapidity distribution of the i-th particle [20, 21], where $ p_1 $ and $ p_2 $ are 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 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^{ \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. 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. [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 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.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 $ \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} \underline{\underline {\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}\underline{\underline {\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 do the following AME,

      $ \dfrac{1}{{\cal{D}}^+_\alpha+ {\rm{i}} \eta}\underline{\underline {\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^+}\underline{\underline {\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 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

      $ {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}} $ include only single-scale vacuum bubble integrals, and $ F^{\rm{cut}} $ denote basal phase-space integrations with integrand being polynomials of scalar products between cut momenta. $ F^{\rm{bub}} $ have been well studied up to five-loop order (see [37-43] and references therein). $ F^{\rm{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. [44] 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 [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 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} $.

      We emphasize that, in practice, one should first reduce 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. But it makes the step of setting up DEs w.r.t. AM $ \eta $ more complicated than 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 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

      $ \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}\} $ depend on $ \vec s $, we choose $ \eta^{ \rm{t}}_1 = 0^+ $ and $ \eta^+_1 = \eta $; Otherwise, we choose $ \eta^{ \rm{t}}_1 = \eta $ (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 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. [63-66]) 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., $ \gamma^*\to t\bar{t}g $ process, 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 obtained 2 physical MIs are $ F^{\rm{cut}}_{2,3,1} $ and $ F^{\rm{cut}}_{2,3,2} $, which are calculated in Appendix.A.

      For VR, 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 = (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 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, 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 =k_3^2, {\cal{D}}^{\rm{c}}_4 = (Q-k_1-k_2-k_3)^2;\,\\ {\cal{D}}^{{\rm{t}}}_1 = &(Q-k_1-k_2)^2, {\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, {\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, {\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 $. Then phase-space integrals can 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)

      By 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 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} $ process, where (a) and (b) define the two most complicated family and (c) defines a sub-family of (b). Here thick curves represent top quark, thin curves represent massless particle, and vertical dashed lines represent 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 of MIs. For brevity, we define 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] &\{\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)\}. \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) \underline{\underline {\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) \underline{\underline {\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 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 obtain $ \hat F(\{0,0,1\};x,0) $ at any value of x.

    • C.   VRR

    • For VRR, 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

      $ \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}}} \int\frac{ {\rm{d}}^Dl^+_1}{(2\pi)^D}\prod\limits_{\alpha = 1}^{6}\frac{1}{({\cal{D}}^+_{\alpha}+ {\rm{i}}\eta)^{\nu^+_{\alpha}}}. $

      (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 most complicated family is the first family, which can also be characterized by Feynman diagram 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). Here thick curves represent top quark, thin curves represent massless particle, and vertical dashed lines represent 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 of MIs. For brevity, we define 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 of these DEs is given by

      $ \begin{aligned}[b] &\hat F\left(\{1,1,1\};\frac{1}{2},y\right) \underline{\underline {\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) \underline{\underline {\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 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.,

      $ \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 is 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), solving the above DEs 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, 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 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 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^{\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. Here thick curves represent top quark, thin curves represent massless particle, and vertical dashed lines represent final state cut.

      For VRV, 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 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. $ \eta $ and that 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 literature [67, 68].

      For 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 literature [63] ($ T_4 $ and $ T_5 $ in this reference).

    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^-\rightarrow \gamma^* \rightarrow t\bar{t}+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. $ \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 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 DE 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.

    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 MIs of basal phase-space integration with no denominator. It is important to note that, besides masses presented in cut lines, these MIs only depend on the square of 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 will provide explicit results for $ r = 0,\; 1 $, and $ 2 $. MIs for general cases can be studied as following. By using optical theorem (see, e.g., Ref. [67]), 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. [69].

      For $ r = 0 $ and $N\geqslant 2$, theres 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 by using the public 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 are 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 partial results in literature [see Eq.(B.7) in Ref. [72] and Eqs. (3.14-3.16) in Ref.[63]].

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

      $ \begin{aligned} {\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, the $ \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 two-particle cut is used:

      $ \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}. \tag{A5} $

      (A5)

      The following relations (can be found, e.g., from 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)\\ &+ \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\left(1-a,1-b,c-a-b+\alpha;c-a-b+1,c-a-b+\alpha+\beta;\frac{z}{y}\right)\\ &+ \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;\right.\\ &\left.\quad 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