Automatic computation of Feynman integrals containing linear propagators via auxiliary mass flow

We proposed a recipe to systematically calculate Feynman integrals containing linear propagators using the auxiliary mass flow method. The key of the recipe is to introduce a quadratic term for each linear propagator and then using differential equations to get rid of their effects. As an application, we calculated all master integrals of vacuum integrals containing a gauge link up to four loops, and we checked the results by nontrivial dimensional recurrence relations.


I. INTRODUCTION
As particle physics enters the era of precision, it is necessary to compute perturbative quantum field theory to high order in coupling constants. To this end, the evaluation of Feynman integrals is an unavoidable task. The main strategy at present is to decompose all Feynman integrals in a problem to a small set of bases, called master integrals (MIs), using integration-by-parts (IBP) reduction , and then calculate these MIs.
Many methods exist to calculate MIs [2,. Among them, the auxiliary mass flow (AMF) method [57][58][59][60], a special case of differential equations (DEs) method, calculates MIs by setting up and solving DEs of MIs with respect to an auxiliary mass term η (denoted as η-DEs). This is systematic and efficient, as far as η-DEs can be set up by using reduction strategy. Equipped with an iterative strategy [59] and the blocktriangular reduction form [11,12], the AMF method becomes so powerful that many Feynman integrals in cutting-edge problems, which are very challenging for other methods, can be calculated (see Fig. 4 in Ref. [59]). At the same time, the AMF method has already been used in many physical processes [61][62][63][64][65][66][67][68].
The AMF method is systematic only if boundary conditions of the η-DEs at η → ∞ can be systematically obtained. If all inverse propagators defining Feynman integrals are quadratic in loop momenta, where these kinds of Feynman integrals are called quadratic integrals, then regions at η → ∞ can be systematically identified [69,70], and the contribution of each region can be calculated within the AMF framework [71]. However, life becomes not that easy when some inverse propagators defining Feynman integrals are linear in loop momenta, where these kinds of Feynman integrals are called linear integrals. In this case, identification of regions at η → ∞ is not that systematic, and the calculation of each region becomes nontrivial [61,62]. As linear integrals show up frequently in region expansion [69] and in effective field theories, it will be good if they can be calculated systematically by the powerful AMF method.
In this paper, by introducing an auxiliary quadratic term for each linear propagator, we develop a recipe to calculate any linear integral systematically within the AMF framework. In the rest of the paper, we first present the general idea of our recipe and then give some examples to illustrate its details. As an application, we calculate all master integrals of four-loop vacuum integrals containing a gauge link for the first time, which are useful for studying parton distribution functions [72][73][74][75][76]. All results in this paper have been checked by using dimensional recurrence relations [36,77,78].

Consider a general Feynman integral
where L is the number of loops and N is the number of propagators and irreducible scalar products, and D a are inverse propagators and ν = (ν 1 , ν 2 , · · · , ν N ) integers. We introduce a corresponding auxiliary quadratic integral I ν (x) defined by If D a is an irreducible scalar product or is quadratic in loop momenta, then D a = D a . Otherwise, it can be generally expressed as 1 Taking advantage of the DEs (7), there are only n independent J s i j k , which can be fixed by matching the expansion (8) with values at a given point − → J (x 0 ), as far as x 0 is chosen sufficiently small so that the above expansion is convergent at x = x 0 . Therefore, once − → J (x 0 ) is known, the asymptotic expansion of − → J (x) around x = 0 can be achieved, and thus, I ν can be obtained.
For a given regular point x = x 0 , − → J (x 0 ) are quadratic integrals that can be calculated by using the AMF method [57][58][59][60]. A convenient choice can be replacing i0 + by −η + i0 + in all modified propagators, i.e., propagators which are linear in loop momenta when x → 0 + . The advantage of this choice is that the number of MIs does not increase after introducing η. As η → ∞, the newly defined integrals in the AMF method are essentially single-mass vacuum integrals [59], which can be easily calculated within the AMF framework [71]. With this boundary condition, we can obtain − → J (x 0 ) by flowing η from ∞ to i0 − using η-DEs.
As numerically solving DEs is very efficient, our method is thus systematic and efficient as far as DEs can be set up.

III. A PEDAGOGICAL EXAMPLE
To elaborate how to use the above proposed recipe to calculate linear integrals, we take the simplest linear integral with one linear propagator and one quadratic propagator as an example, where ℓ is the loop momentum, and q and n are external momenta with n 2 = −1. Clearly, this integral only depends on the scale q · n. Without loss of generality, we choose q · n = 1. This Feynman integral is relevant to the Feynman diagram shown in Fig. 1, where the double line stands for a gauge link which contributes linear propagators. We call this a one-loop bubble integral because there is no external legs of particles. The corresponding auxiliary quadratic integral is defined as with For finite x, there are two MIs which can be defined as with x . Using IBP, we get the system of x-DEs, Besides x = ∞, there are two additional singularities x = 0 and x = − 1 4 , which correspond to an infinity mass term and zero mass term in D 2 , respectively. Therefore, we can choose boundary conditions at, e.g., x 0 = 1 16 , which point can be well estimated by asymptotic expansions around x = 0. Based on DEs (17), we get the following asymptotic expansions: which gives As x → 0 + , we find J 1 (0 + ) = 0 and I (1,1) (0 + ) = c 0 . The DEs (17) also provide recurrence relations of coefficients, where i ≥ 1. Therefore, only a 0 and c 0 are free coefficients.
To fix a 0 and c 0 , we need to know where we only keep 10 orders in ǫ expansion with 16-digit precision, although higher order and more digits can be obtained easily. We keep the same accuracy in all examples in this paper. Calculation of − → J ( 1 16 ) in this example is certainly a very simple problem, but for general complicated problems we can still use the AMF method to calculate them. By matching the asymptotic expansions , we can determine a 0 and c 0 . Eventually we obtain

IV. VACUUM INTEGRALS CONTAINING A GAUGE LINK UP TO FOUR LOOPS
Now, we apply our recipe to calculate vacuum integrals containing a gauge link up to four loops. These results are useful for studying parton distribution functions [72][73][74][75][76]. On the one hand they can serve as boundary conditions for four-loop parton correlation functions calculation, and on the other hand, they are needed to renormalize three-loop parton correlation functions [79,80]. In general, for an L-loop vacuum integral with a gauge link, there must be at least L quadratic propagators and at most L linear propagators, or else, it is either reducible or scaleless. With this restriction, up to four loops we find that there are 64 nonzero and irreducible sectors in total which give 65 MIs. The 64 nonzero and irreducible sectors are displayed in the the supplementary material together with results of all the 65 MIs. We take the most complicated integral family shown in Fig. 2 as an example to explain our calculation details. All other families can be calculated similarly and final results are provided in the supplementary material.
A general integral in this family can be expressed as where D 1 , D 2 , · · · , D 8 are quadratic denominators, D 9 , D 10 , D 11 linear denominators, and D 12 , D 13 , D 14 irreducible numerators, with explicit form where Feynman prescription +i0 + is omitted. This family has 32 MIs in total. Similar to the one-loop case, the single scale q ·n can be extracted as an overall factor, and thus we can simply set q · n = 1.
Based on our recipe, we define a corresponding auxiliary family of integrals, where and D i = D i for i = 9, 10, 11. After changing linear propagators to quadratic propagators, the number of MIs in the auxiliary family increases to 120. We use − → J (x) to denote these MIs. Each J i (x) is defined by a specific list of exponents ν i via the following relation where, for convenience, the factor x ν i 9 +ν i 10 +ν i 11 is introduced to cancel the overall factor x in D 9 , D 10 , D 11 .
For a regular point, say x 0 = 1 16 , we calculate the corresponding MIs − → J ( 1 16 ) using the package AMFlow[60] with the mass mode at the first step of the iteration. The calculation is automatic and can be decomposed to the following steps. First, i0 + in D 9/10/11 is replaced by −η + i0 + . Then, η-DEs are set up, and boundary conditions at η → ∞ are calculated [59,71], where Kira [18] is employed to do IBP reduction. Finally, the η-DEs are solved to realize the flow of η from ∞ to i0 − . Note that higher precision and higher orders in ǫ can also be achieved easily.

V. INTEGRALS WITH EXTERNAL LEGS
Now let us consider more general linear integrals containing both gauge links and external legs. Take fourloop integrals containing one external leg shown in Fig.  3 as an example. The top-sector corner integral can be expressed as 1,1,1,1,1) with where p is a lightlike external momentum, and Feynman prescription +i0 + for each inverse propagator is omitted. To calculate integrals in this example, we define auxiliary quadratic integrals by using the following inverse propagators and D i = D i for i ≤ 4. After changing linear propagators to quadratic propagators, the number of MIs increases from 2 to 6. Using AMFlow [60], we calculate auxiliary quadratic integrals at a fixed point x = x 0 , and then solve x-DEs to obtain integrals at x = 0 + . Without loss of generality, for the phase-space point q ·n = 1, p·n = 1 2 , we obtain 1,1,1,1,1) When there are more gauge links or more external legs, the calculation process is similar.

VI. SELF-CONSISTENCY CHECK BY DIMENSIONAL RECURRENCE RELATIONS
For any integral I D ν , where we explicitly write the dependence of spacetime dimension D, there is a parametric representation (see, e.g., [81]), and thus, As U is a homogeneous polynomial of x i , the above expression can be translated into a linear combination of integrals at dimension D with different ν, which can be reduced to MIs by using IBP identities. Supposing J D is a set of MIs, the above procedure results in dimensional recurrence relations [77,78], where C D is a matrix rationally depending on D.
Up to three loops, vacuum integrals containing a gauge link have been calculated in Ref. [80] by solving dimensional recurrence relations using the method in Ref. [36]. The obtained results agree with ours. For four-loop integrals containing a gauge link, we find it very hard to solve the corresponding dimensional recurrence relations. However, the dimensional recurrence relations can provide a highly nontrivial self-consistency check for our results. We calculate MIs at two different space-time dimensions D 0 and D 0 − 2, and we find that the obtained results satisfy the relations Eq. (36) within precision.

VII. SUMMARY
In summary, we develop a recipe to calculate linear integrals using the AMF method [57][58][59]. For any given linear integral, our recipe is to introduce an auxiliary quadratic term for each linear propagator, and the obtained auxiliary quadratic integral can be calculated systematically using the AMF method. Taking the result of the auxiliary quadratic integral calculated at fixed auxiliary quadratic terms as the boundary condition and using differential equations to push the auxiliary quadratic terms to zero, effects of auxiliary quadratic terms will die out eventually, and we get the result of the target linear integral. This recipe of calculating linear integrals is very systematic and has been implemented in the package AMFlow [60].
As linear integrals show up frequently in region expansion and in effective field theories, our recipe will be useful in phenomenological studies. As the first application, we have calculated all MIs of vacuum integrals containing a gauge link up to four loops, which are useful to study parton distribution functions. Our results have been checked by nontrivial dimensional recurrence relations.

VIII. ACKNOWLEDGMENTS
We thank K.T. Chao, X. Liu For vacuum integrals containing a gauge link, x-DEs can be obtained from η-DEs. We take the family in Fig. 1 as an example to explain this, and the final relation is independent of the specific family. MIs after introducing where Feynman prescription +i0 + in each denominator is omitted. The corresponding x-DEs are denoted as When calculating − → J (x 0 ) using AMF method, we define modified integrals The corresponding η-DEs are denoted as We express A(x) by M (η). Let us define − → K (p 2 , m 2 ) = K 1 (p 2 , m 2 ), K 2 (p 2 , m 2 ) .
Because Feynman integrals have a definite mass dimension, there is a rescaling relation, where α is a diagonal matrix with diagonal elements denoting mass dimension. In this example, it is given by (A10) Expansion of Eq. (A9) to O(δ) gives For convenience, we denote which connects x-DEs to η-DEs for any vacuum integrals containing a gauge link. It is useful because M can be more easily obtained.