A Systematic and Efficient Method to Compute Multi-loop Master Integrals

We propose a novel method to compute multi-loop master integrals by constructing and numerically solving a system of ordinary differential equations, with almost trivial boundary conditions. Thus it can be systematically applied to problems with arbitrary kinematic configurations. Numerical tests show that our method can not only achieve results with high precision, but also be much faster than the only existing systematic method sector decomposition. As a by product, we find a new strategy to compute scalar one-loop integrals without reducing them to master integrals.


Introduction
With the continuous improvement of statistics and experimental systematics at the Large Hadron Collider, the aim of testing the particle physics Standard Model and discovering new physics strongly demands theoretical predictions to also improve uncertainty to the same level. For many important processes, high order perturbative calculations are needed to this end. At the one-loop order, thanks to the improvement of traditional tensor reduction [1] and the development of unitarity-based reduction [2][3][4], one can efficiently express scattering amplitudes in terms of linear combinations of master integrals (MIs). As the computation of oneloop MIs is a solved problem [5][6][7], one-loop calculations can now be done automatically. Expressing multi-loop scattering amplitudes in terms of MIs is also possible using such as the integrationby-parts (IBP) reduction [8][9][10][11][12][13] or the unitarity-based multi-loop reduction [14][15][16][17][18][19][20][21][22][23][24][25][26]. Then, one of the main obstacles for multi-loop calculation is the computation of multi-loop MIs.
We take two recent studies in literature as examples to demonstrate how hard the multi-loop MIs computation is. One example is a two-loop calculation of pseudoscalar quarkonium inclusive decay [27], where the computational expense of MIs is about O(10 5 ) CPU core-hour. Another example is a calculation of four-loop nonplanar cusp anomalous dimension [28]. The reduction of amplitudes to MIs in this problem has been done much earlier in Ref. [29], yet the computation of these MIs is very challenging. The final numerical result obtained in Ref. [28] has uncertainty as large as 10%, which we believe is already the best precision that one can get with a tolerable computational expense.
Currently, the only method that can systematically compute any MI is the sector decomposition [30]. Unfortunately, this method is extremely time-consuming, besides that it is hard to achieve high precision. Mellin-Barnes representation [31] is another widely used method, yet it has difficulty to deal with non-planar diagrams, at least not in a systematic way (see Ref. [32] and references therein for recent progress). The differential equation (DE) method [33][34][35][36] is a powerful tool to compute multi-loop MIs, which bases on the fact that derivation of a MI with respect to its kinematic variables (including Mandelstam variables and internal masses) can be re-expressed as a linear combination of MIs using aforementioned reductions. For simple problems, DE method can give analytical results thanks to the introduction of canonical form [37][38][39]; while for complicated problems, one can solve DEs numerically to achieve results with high precision (see [40][41][42][43] and references therein). However, it needs input of boundary conditions (BCs) of MIs evaluated at another set of kinematic configurations. As there is no general rule to obtain BCs for arbitrary problems at present, one needs to find good BCs case by case, which makes it hard for DE method to be systematical. In practice, sector decomposition method is employed in Ref. [27], and both sector decomposition method and Mellin-Barnes representation method are employed in Ref. [28]. In this Letter, we develop a novel method to compute multiloop MIs by constructing and solving a system of ordinary differential equations (ODEs). Advantages of our method are as follows: 1) Our BCs are fully massive vacuum integrals with a single mass scale, which are much simpler to compute and have been well studied in literature [44]. As a result, our method can be systematically applied to any complicated problem; 2) ODEs can be numerically solved efficiently to high precision, no matter how many mass scales are involved in the problem; 3) Computing MIs with complex kinematic variables is very easy in our method, while it could be hard for other methods (note that introducing imaginary part to kinematic variables is needed for many purposes, e.g., to describe particle decay or to study the S-matrix theory). Numerical tests show that our method can be much faster than the only existing systematic method sector decomposition. As a by product, we find a new strategy to compute scalar one-loop integrals in arbitrary spacetime dimensions without reducing them to MIs.

The method
Let us introduce a dimensionally regularized L-loop MI, where D is the spacetime dimension, D α ≡ q 2 α − m 2 α are usual Feynman propagators, and q α are linear combinations of loop momenta i and external momenta p i . The actual integral that we want to get is with 0 + defining the causality of Feynman amplitudes. In the following, we will suppress the dependence on D and {ν α } whenever it does not introduce any confusion.
We set up ODEs by differentiating MIs with respect to η and then re-expressing them in terms of MIs, which results in where I(η) is the vector of a complete set of m MIs and A(η) is the m × m coefficient matrix. To obtain MIs at η = 0 + , we solve the ODEs with BCs chosen at η = ∞. As we will show, BCs are simply vacuum integrals with equal masses, which can be computed rather easily. Considering also that numerically solving these ODEs is well-studied mathematical problem, our method provide a systematic and efficient way to compute multi-loop MIs to high precision.

Boundary conditions
Before studying BCs rigorously, let us explain the idea of choosing BCs at η = ∞. With a sufficiently large imaginary part in all denominators, we expect all kinematic variables to be negligible because they are finite. Thus we should be able to set both internal masses m α and external momenta p i to zero at the boundary, which results in simple vacuum integrals with equal masses. The only loophole in this argument is that, as loop momenta i can be arbitrarily large, it is not obvious that i · p j are negligible comparing with η even if η → ∞. The loophole can be fixed by studying its Feynman parametric representation, and then our naïve expectation holds in general.
We assume ν α > 0 for all α in Eq. (1) to simplify our discussion, although our final conclusion is unchanged even without this where U and F are so-called graph polynomials that can be related to the spanning 1-tree and 2-tree of the original Feynman diagram, respectively (see e.g. Ref. [45]), and ν is short for α ν α .
Feynman parameters. By definition, a 2-tree can be generated by a 1-tree, i.e. there exists a pair of indexes j and k so that where I bub (1) can be interpreted as a vacuum integral with equal internal squared masses m 2 = −i. It is worth mentioning that the object J (η) ≡ η ν−LD/2 I(η) is analytic near η = ∞ based on the above discussion.
To compute I bub (1), we again reduce it to linear combination of corresponding vacuum MIs, diagrams of which up to 3 loops are shown in Fig. 1. Computation of these vacuum MIs is well studied, with analytical results available up to 3 loops [46][47][48] (see [49] and references therein for some pioneering works) and numerical results available up to 5 loops [44,50,51]. We therefore conclude that the computation of BCs in our method is a solved problem.

Solving ODEs numerically
Knowing BCs, solving the ODEs numerically to obtain MIs at η = 0 + is a well-studied mathematical problem. The solution can be obtained efficiently to high precision. Having located all singular points, an algorithm to solve the system of ODEs can be roughly divided into three steps: • Step 2: Using ODEs to take Taylor expansion of MIs around η = η j ( j = 0, 1, · · · ) and obtain their values at η = η j+1 . Do it repeatedly until η j+1 is inside of the small circle, denoting η = η last in Fig. 2. • Step 3: Using ODEs to take asymptotic expansion of MIs around η = 0 and matching with their values at η = η last to obtain dimensionally regularized MIs at η = 0 + . Modified MIs at the step 1 are combinations J k (η) = η ν k −LD k /2 I k (η), which are analytic around η = ∞ as pointed out above. From BCs at η = ∞, ODEs tell us their first-order derivatives, based on which ODEs then tell us their second-order derivatives, and so on. Eventually we get the Taylor series at η = ∞.
Taylor expansion at the step 2 can be obtained similarly.
Having obtained the values of I(η last ) with η last < η min from step 2, we then do an asymptotic expansion of MIs at η = 0. We In general, the matrix Ã can be divergent as η −p with positive integer p, which is called the Poincaré rank at η = 0. However, because MIs have integral representation, they can only have regular singularities on the whole complex plane of η. This enables us to reduce the Poincaré rank at η = 0 by a transformation (see Ref. [38] and references therein for different algorithms). So, without loss of generality, we can assume that Ã (0) is finite. Then the general solution of Eq. (6) near η = 0 has the form 1 where the matrix P (η) has been expanded as ∞ n=0 P n η n with P 0 being the identity matrix. In the above expansion, non-analytical behaviors of η, like η a+bD or ln(η), are contained in the matrix exponential exp Ã (0)ln(η) = η A(0) . Substitute Eq. (7) into Eq. (6), we get a system of equations, which can be solved recursively to obtain P n (n ≥ 1). It is clear that P (η) is determined completely by the matrix Ã (η), or equivalently by the ODEs. As a consequence, all boundary information are contained in the vector v 0 , which can be determined uniquely by matching the value of I(η last ) obtained from step 2, Then we have For dimensional regularized Feynman integrals, all terms that are non-analytic in η vanish in the limit η → 0 + , which results in well-defined I(0). All infrared divergences emerge in this limit. As a result, our method works regardless of whether propagators in the Feynman integral are massive or massless. In practice, instead of using Eq. (10) directly, a more efficient way is to solve nonhomogeneous ODEs sequentially.

An one-loop example
A fascinating feature of our method is that ODEs at one-loop level can be easily constructed for all scalar integrals, not restricted to MIs. We will only consider scalar integrals with nonvanishing Gram determinant, because scalar integrals with vanishing Gram determinant can be easily reduced to the former case. Using a raising dimensional recurrence relation and a lowering dimensional recurrence relation [see e.g., Eqs. (14) and (18) in Ref. [53]], we get ∂ ∂η  α=1 R βα z α = C (β = 1, · · · , N). It is worth mentioning that both C and z α are independent of η. Therefore, Eq. (11) has only a purely imaginary singularity at η = i 2 C , which is the leading Landau singularity of The DE (11) naturally connects a D dimensional scalar integral to some (D − 2) dimensional scalar integrals with one less propagator. Thus, we should better choose this set of scalar integrals as basis to solve the system of ODEs. Based on this, we can compute any one-loop scalar integral.
As an example, we compute an infrared-divergent multi-scale scalar integral shown in Fig. 3(a), with the cross symbol indicating that the power of the corresponding propagator is two.
We set D = 4 − 2 and express the final results as Laurent series up to O( ). Kinematic variables are chosen as p 2 s is close to the threshold (m 3 + m 4 ) 2 , and then we have η min ≈ 4.5|δ| and η max ≈ 10.2. We can run around 3 + log 2 (η max /η min ) ≈ 4.2 + 1.4 ln(1/|δ|) steps to go from η = ∞ to η = 0 + . Thus the computational complexity increases slowly as s close to threshold.
By taking Taylor expansion or asymptotic expansion to 30 orders at each step of the running process, we find that the obtained final numerical results for the above scalar integral have at least 10 correct significant digits for any choice of δ in the range 10 −7 ≤ |δ| ≤ 1. To compare with, we find that the sector decomposition approach FIESTA4 [54] can get result with tolerable uncertainty only for |δ| ≥ 10 −3 .

A two-loop example
We take the non-planar scalar integral shown in Fig. 3(b) as a two-loop example, where solid lines have masses m 2 = 1, dashed lines are massless, s = 4, and t = −1. This is the most timeconsuming MI in the study of Ref. [27].
For multi-loop integrals, similar relations as Eq. (11) are not available yet. One needs to rely on reduction methods. We use FIRE5 [12]  With η min ≈ 0.118 and η max = 4, we run 14 steps to go from η = ∞ to η = 0 + . By taking Taylor expansion or asymptotic expansion to 30 orders at each step, we get numerical result for the non-planar integral in Fig. 3(b), − (4.75677 − 12.0749i) + (23.9674 − 55.4214i) + · · · , (12) together with results of other 167 integrals in the same sector, all of which have at least 6 correct significant digits. We obtain these numerical results using our Mathematica code in a few minutes on a laptop. The efficiency can be further improved in a future Fortran code. To compute MIs in the same sector using FI-ESTA4 in Ref. [27], it cost about O(10 4 ) CPU core-hour to achieve results with 1‰ uncertainty. 2 Therefore, our method is faster by at least 10 5 times than FIESTA4 for this problem. Using finite integrals as bases [55] or using other sector decomposition codes, e.g. pySecDec [56], may speed up the calculation for this problem, but the efficiency should be still hard to compete with our method.

Summary and outlook
In this Letter, we presented a novel method to compute multiloop MIs, which works by constructing and solving ordinary differential equations of MIs with respect to Feynman prescription parameter η. On the one hand, the method retains all advantages of usual DE approach, such as efficient for numerical evaluation, stable around threshold region, and flexible to introduce imaginary kinematic variables; and on the other hand, the method overcomes the main difficulty of usual DE approach, with almost trivial boundary conditions. Therefore, our method can compute MIs systematically and efficiently to high precision, which makes it possible to define MIs as new special functions instead of expressing them in terms of other special functions.
Based on a two-loop example, it was found that our method can be much faster than the only existing systematic method, sector decomposition. With this significant improvement of efficiency, one can now use reasonable computer source to perform phase space integration in multi-loop calculations. To this purpose, one needs to construct and solve ODEs of MIs at each given phase space point, which is very natural in our method.
One thing which we want to clarify is that, to estimate the time consumed in our method, we only counted the time to solve ODEs, but ignored the time to construct ODEs. The reason is as follows. Whenever one wants values of MIs for a multi-loop calculation, one must have expressed scattering amplitudes as linear combination of them by reduction methods. The construction of ODEs of MIs is usually much faster than the full reduction of scattering amplitudes, and thus its time can be neglected. We further note that a new reduction method is proposed in Ref. [57], which may significantly reduce the time to construct ODEs.
Many variants of our method can be constructed easily. The basic idea of our method is to introduce an infinite large value in propagators at the boundary of ODEs, so that all kinematic variables can be dropped out. This aim can be equally met, e.g., by only differentiating the parameter η of some propagators in the construction of ODEs. Then to compute corresponding BCs, one again constructs ODEs but with respect to η of other propagators.