$\texttt{AlgRel.wl}$: Algebraic Relations for the Product of Propagators in Feynman integrals

Motivated by the foundational work of Tarasov, who pointed out that the algebraic relations of the type considered here can lead to functional reduction of Feynman integrals, we suitably modify the original method to be able to implement and automatize it and present a $\textit{Mathematica}$ package $\texttt{AlgRel.wl}$. The purpose of this package is to help derive the algebraic relations with arbitrary kinematic quantities, for the product of propagators. Under specific choices of the arbitrary parameters that appear in these relations, we can write the original integral with all massive propagators in general, as a sum of integrals which have fewer massive propagators. The resulting integrals are of reduced complexity for computational purposes. For the one-loop cases, with all different and non-zero masses, this would result in integrals with one massive propagator. We also devise a strategy so that the method can also be applied to higher-loop integrals. We demonstrate the procedure and the results obtained using the package for various one-loop and higher-loop examples. Due to the fact that the Feynman integrals are intimately related to the hypergeometric functions, a useful consequence of these algebraic relations is in deriving the sets of non-trivial reduction formulae. We present various such reduction formulae and further discuss how, more such formulae can be obtained than described here. The $\texttt{AlgRel.wl}$ package and an example notebook $\texttt{Examples.nb}$ can be found at https://github.com/TanayPathak-17/Algebraic-relation-for-the-product-of-propagators


Introduction
In this work, we consider the formalism first proposed by Tarasov to derive algebraic relations for the product of propagators for functional reduction [1]. We systematically develop an algorithm inspired by the original work and present a realization in Mathematica for the same, which is provided for the user as a package called AlgRel.wl . We have used the package to simplify and analyze many important and interesting Feynman Integrals that are amenable to treatment using this formalism.
Feynman integrals play an important role in precision calculations in quantum field theory. There are various methods to evaluate them [2,3]. Even with all these methods, it is at times still challenging to compute Feynman integrals. More often, other techniques are used to facilitate this computation. In [4] the method of functional reduction is introduced to derive functional relations between Feynman integrals. These relations reduce the original integral into a sum of integrals which are easier to evaluate. The focus of the present work is this new way to obtain functional relations by deriving the algebraic relations for the product of propagators. This method in turn then leaves some undetermined free parameters which can be chosen at will. Appropriate choices of these parameters result in various functional equations for Feynman integrals [5][6][7][8][9].
The method can be applied to any one-loop diagram, indeed as already pointed out in detail by Tarasov. Despite this, no working code has been provided in the past. In the present work, we provide an automated Mathematica package AlgRel.wl to derive the algebraic relation for the product of propagators. Our code here fills this gap in the possibility of finding widespread use of formalism. Since our goal is an efficient algorithmic implementation to find the algebraic relation, we introduce a recursive way of method. The free parameters in the resulting relation can then be chosen in an appropriate manner to derive the functional equations for the Feynman integrals. More specifically, for presentation purposes, we focus on the cases when all these free parameters are zero and the original Feynman integral with many massive propagators can be written as a sum of integrals with fewer massive propagators, which was also pointed out in [10] 1 . For the one-loop integrals, with all different and non-zero masses, this procedure can be used to reduce the original integral to a sum of integrals with one massive propagator. We apply the method for up to 6-point, one-loop integrals and show that the N −point one loop integral with all massive propagators and general external momenta can be written as a sum of 2 N −1 integrals with just one massive propagator. Though the method is not readily generalizable to the higher loop we yet extend the uses to cover certain cases of 2-and even 3-loop Feynman integrals. In a similar manner, this approach is also applicable to higher loops. Our findings show that we require at least 4 propagators in order for the formalism to be viable and to be of utility as far as the simplification is concerned. We explain this feature in some detail.
We, however, notice that such functional reduction is one of the many possibilities obtained after choosing the free parameters obtained from the algebraic relation [11]. In view of the proposed method of functional reduction of Feynman integrals, the package has been built in such a way that the final result still has arbitrary parameters which can be chosen suitably for the functional reduction procedure. Using a few of the analytical results available for the one-loop integrals, we explicitly show how the complexity in the evaluation of these integrals can be reduced. Whenever the Feynman integrals can be expressed in terms of hypergeometric functions [12][13][14][15] this reduction in complexity gives rise to reduction formulae for the hypergeometric functions. Thus it can be used to establish new relations between multi-variable hypergeometric functions. We discover many new reduction formulae for such hypergeometric functions, which, to the best of our knowledge, have not appeared anywhere in the literature. We also discuss in detail how further reduction formulae can be obtained from already available results for the one-loop cases. Such relations between hypergeometric functions are also obtained in [16], where explicit relations between hypergeometric functions are derived via the evaluation of Feynman integrals. In order to make the results accessible, we provide several examples in a single Mathematica notebook that allows the reader to appreciate the power of the formalism, based on the code that is provided along with it.
The article is organized as follows. In section 2 we discuss the method in detail with one loop bubble integral as an example and explicitly present how the reduction in complexity has been achieved for the integral. In section 3, we present the algorithm of the AlgRel.wl package and discuss its usage in detail. In section 4, various results obtained for one, two, and three-loop integrals are presented. In section 5, we discuss the various analytic results in terms of multi-variable hypergeometric functions already derived for the one-loop N −point integrals [17,18] and show how the present work helps in deriving the reduction formulae for the multi-variable hypergeometric functions using them. Finally, we conclude the paper with summary and discussions in section 6. In Appendix B, we provide a list of various reduction formulae that we derive, along with some details on how to further extend the list given there.
The package AlgRel.wl along with a Mathematica notebook Examples.nb , that contains all the examples discussed in the paper can be found in the GitHub repository.

The Method
We now explain the method to find the algebraic relation of the product of propagators with the help of the one-loop bubble integral. For this example and all the following examples we consider the integrals in d dimensions. Consider the one-loop bubble integral To find the algebraic relation for the product of propagators, we instead consider a more general propagator, depending on only one loop-momenta, of the following form where k is the loop-momentum, q i 's are dependent on external momenta and can be zero as well and m i is the mass of the propagator.
With the general propagators, we now have where substituting q 1 = 0 and q 2 = −p we recover Eq.(2.1). We seek the algebraic relation for the integrand, by introducing a new denominator D 1 along with coefficients x 1 and x 2 , of the following form 2). The unknowns that are introduced can be fixed using the above equation, while the remaining parameters are arbitrary and can be fixed at will in such a way that the resulting relationship gives rise to integrals which are easier to compute. Using Eq.(2.4) we get Comparing the coefficients of k 2 , k and the remaining k independent term we get Solving for x 1 , x 2 and P 1 we get following two sets of solutions and We remark that both of the above sets can be used for the purpose of finding the algebraic relation. However, for convenience, we would use the set given by Eq.(2.7). For the purpose of the algebraic relations, both of the choices are equivalent. One can use either of the solutions and can check that it satisfies the algebraic relation that we seek to find. In the above equation, M 1 is still an arbitrary variable that can be chosen at will. Choosing various values of M 1 results in different functional equations [11] for the bubble integral.
For the present work, we focus on one of the simple choices i.e. M 1 = 0. Integrating Eq.(2.4) and substituting q 1 = 0 and q 2 = −p we have Hence we see that the general two-point bubble integral with non-zero masses can be written in terms of two integrals with just one mass. Diagrammatically Eq.(2.9) can be represented as in Fig.2. To see how the complexity in the computation has been reduced in the Eq.(2.9), we refer to a few analytic results. The general result for the massive bubble diagram can be written in terms of the Appell F 4 function [19]. where, is the Appell F 4 hypergeometric series with region of convergence (ROC) given by |x| + |y| < 1. The analytic expression result for I 2 (p, m, 0) is readily available in [20,21].
p 2 m 2 (2.12) Using the above relation in Eq.(2.9), we get the following for the right-hand side The above relation can be viewed as a reduction formula without making reference to the underlying Feynman integral and the result is shown in Eq.(5.4) and Eq. (5.7). In a similar manner, evaluation of other Feynman integrals can be used to obtain the relationship between hypergeometric functions [16]. Such a reduction of hypergeometric functions with a higher number of variables to those with a lesser number of variables also helps when the analytic continuation has to be done to reach a certain kinematical region. For the case of Appell F 4 the elaborate analytic continuation has been performed explicitly in [22] or using automatized algorithms [23] for more general multi-variable hypergeometric functions. This whole process still does not guarantee convergence for all the values of the parameter space [24]. While for the case of 2 F 1 complete table of analytic continuations is available [25]. The procedure to find the analytic continuations also gets more complicated with the increase in the number of variables even with the use of automatized packages.

Algorithm
We now present a general algorithm for the case when we have N denominators to find algebraic relation recursively.
Consider the general situation with product of N denominators as 1 d 1 · · · d N .
1. We first find the algebraic relation by taking d 1 and d 2 2. We then multiply the above equation by 1 3. We then find the algebraic relation of each pair of d i s again using Eq.(3.1).
4. Then in the resulting relation, we repeat this process until all the denominators are exhausted.
The final result is a sum of 2 N −1 terms where N is the total number of denominators we started with. It is to be noted that the above procedure is a slight modification of the original method [1]. In [1], we start by seeking the following algebraic relation for the product of N propagators Comparing the coefficients of k 2 , k and using the constant term we get an over-determined set of equations. Such a system leave x 3 , x 4 · · · x N undetermined. Such procedure when used recursively with each term on the RHS of the above equation finally results in N ! total number of terms, unlike 2 N −1 terms using the procedure presented here. Also, the arbitrariness in the choice of coefficients x i s in the original algorithm is now present in the choice of parameters M i s.

Usage
The recursive algorithm presented previously has been automatized in the accompanying Mathematica package AlgRel.wl . Below we demonstrate the usage of the package AlgRel.wl . After downloading the package and putting it in the same directory as the notebook we can call the package as follows: The package has been made assuming the form d i = (k + p i ) 2 − m 2 i for the propagator, where k, p and m can be changed as per the convenience of the user. The only command of the package is AlgRel.wl , which can be called as follows • {k,q,m}: It is a list containing three variables corresponding to the general propa- k denotes the loop momenta, q denotes the combination of external momenta and can be zero too and m denotes the mass of the propagator.
• {P,M}: It is a list containing two variables. They are used to set the variables for the auxiliary propagator introduced for obtaining the algebraic relation, It automatically takes the k from the previous list.
• x: It is used to denote the variable for the coefficients in the algebraic relation, Eq.(3.3).
• Substitutions: It is a list of substitution for q i and M i .
The output of the above command is a nested list with two sub-lists with the following two sub-lists • {Algebraic relation}: It gives the algebraic relation for the product of propagators, Eq.(3.3).
• {Values}: It is a list of the values obtained for P i and x i .
Consider the example of Bubble integral. To obtain the result for it we can use the following command Due to its length, the second element of the output (i.e., the substitution list) is not shown fully. It contain the values of the x[1],x [2] and P [1] as given in Eq.(2.7). We remark that in the output all the vectors such as P [1], q [1] and q [2] will appear in bold. The scalar product q [1].q [2] appears q[1]q [2]. So care must be taken while doing the numerical checks for the same.
In the next section, we look at a few one-loop and two-loop examples where such a procedure is helpful. For cases where it was feasible to perform numerical checks using the integration, we perform them using FIESTA5 [26] and the corresponding results are given in appendix C.

Results
We now look at results for one loop and higher loop cases that are obtained with the help of the AlgRel.wl package. All the results are also presented in the Mathematica file Example.nb.

One-loop vertex integral
We consider the reduction of the one-loop vertex integral corresponding to Fig.3, which is given by We proceed as described in the previous section. We use the generalized propagators and do the substitutions accordingly so the result reduces to Eq.(4.1). This can be done using following command The result is a relation which is a sum of 4 terms, as follows Integrating Eq.(4.2) over loop momenta k we get vertex integral written as a sum of vertex integrals but with just one massive propagator.

One loop box integral
We now consider one loop box integral corresponding to Fig.4 which can be written as We can get the algebraic relation using the following command Substitute q 1 = 0, q 2 = p 1 , q 3 = p 1 + p 2 , q 4 = p 1 + p 2 + p 3 and M i = 0, i = 1 · · · 7 and simplifying we get where the value of unknowns can be obtained from the Mathematica notebook Examples.nb . Integrating Eq.(4.4) over loop momenta k we get box integral written as a sum of 8 box integrals but with just one massive propagator.

One-loop pentagon integral
We can get the algebraic relation using the following command Doing the substitution as before and simplifying we get where the values of x i and P i can be obtained from the Mathematica notebook Examples.nb . Integrating Eq.(4.6) over loop momenta k we get pentagon integral written as a sum of 16 pentagon integrals but with just one massive propagator. The six-point integral corresponding to the Fig.6, is

One loop six-point integral
As in the previous examples we use the following command to obtain the algebraic relations  To illustrate the method for higher loop integrals let us consider an example of the two-loop box integral, corresponding to the diagram Fig.7. The integral is as follows

Two-loop box integral
The propagators are numbered such that i represents the propagator d i . Firstly we find the algebraic relation for the product of propagators numbered 1 and 2, which has only the loop-momenta k 1 we can use the following command The final relation that we obtain, with M i = 0 is( see Examples.nb ) Multiplying both sides of Eq.(4.9) by 1 ((k 1 −k 2 +p 1 +p 2 ) 2 −m 2 5 ) gives the required algebraic relation for the two-loop box integral.

Two-loop double box integral
The propagators are numbered such that i represents the propagator d i . To find the algebraic relation for the product of propagators numbered 1,5 and 6 we can use the following command Substituting value of q 1 , q 5 and q 6 corresponding to the Feynman integral we get Similarly, for propagators numbered 2,3 and 4, we can use the following command which gives the following result after substituting the value of q 2 , q 3 and q 4 corresponding to the Feynman integral All the values of the parameters P i ,Q i ,x i and y i can be obtained from the Mathematica notebook Examples.nb . To get the algebraic relation for the integrand in Eq.(4.11) we multiply Eq.(4.12) and (4.13) together and then multiply both the sides of the equation We see that, unlike the one-loop case, we now have 3 massive propagators in each term. In fact with the present procedure to find the algebraic relation for any two-loop integral with all non-zero different masses, we have at least 3-massive propagators in each integral. Due to this reason, the present procedure won't be helpful for the case of integrals like the sunset integral where there are only 3-propagators. The three-loop ladder integral corresponding to Fig.9 is

Three-loop ladder integral
We use a similar strategy as before for this case too, to obtain the algebraic relation. The result contains 32 terms and is presented in the Mathematica file Examples.nb .

Limitation at higher loops
Using sunset as an example we now demonstrate the limitation of the method for higher loop integrals. As a demonstrative example, we consider the two-loop sunset integral as a starting point. The corresponding diagram is shown in Fig.10. The integral is given by Using the package we find the algebraic relation between the propagators d 1 and d 3 . For this case, we obtain the following values of the coefficients x 1 and x 2 . We further assume M i = 0 We notice that for the case of sunset q 1 = 0 and q 3 = p − k 2 . Because of the fact that the coefficients x 1 and x 2 are dependent on loop momenta k 2 , the final expression after finding the algebraic relation will include x 1 and x 2 in the integral. However, we can still consider the higher loops cases by taking only the propagators which are dependent on only single loop momenta together and propagators which are dependent on more than one loop momenta are excluded.

Reduction of Hypergeometric functions
In this section, we study the examples when the Feynman integral evaluation gives results in terms of hypergeometric functions. The formalism to find algebraic relation for the product of propagators of Feynman integrals can be employed to find relations between hypergeometric functions [16,[27][28][29]. In this section, we point out some analytic results on N −point function [17,18] and various hypergeometric relations that can be obtained from them with the present analysis. It is well-known that the general one-loop N −point function with zero external momenta and different masses (m i , i = 1, . . . , N ), with unit powers of propagators, can be expressed in terms of Lauricella F D function [17] D represents the Lauricella function of L−variables given by Using the method presented, we can write the N −mass integral as a sum of integrals with just one mass, thus N − 1 masses vanish. Then using Eq.(5.1) and (5.3) we can write each of these integrals as a term dependent only on mass m i . The whole result can be then expressed as a sum of terms each dependent on some m i . To evaluate the result of Eq. (5.1), outside its associated region of convergence, one has to explicitly perform analytic continuation which is difficult to obtain at times for multi-variable hypergeometric functions. Such a reduction of the result is helpful when the analytic continuation of the Eq.(5.1) is required. Such a result should also be viewed as a reduction formula for the Lauricella F L D , obtained using a physical problem [16] which can otherwise be hard to obtain.
For a general one-loop N − point function with non-zero external momenta, the general result can be written as a generalized Lauricella hypergeometric function with (N −1)(N +2) 2 variables [18]. For the case of general vertex integral, the result is a generalized Lauricella function with 5 variables [18]. On the other hand, using Eq.(4.2) the result can be written in terms of a hypergeometric function of 3−variables [18]. Comparing Eq. (2.13) and (2.10) we see that the evaluation of bubble integral has reduced from the evaluation of Appell F 4 which has two variables to that of hypergeometric 2 F 1 with one variable. Such a result can be viewed as a general reduction formula without any explicit relation to the Feynman integrals it has been obtained from. Substituting a = = y, we get the following relation y a− 1 F 4 (a, 1, a, a, x, y) − F 4 (2 − a, 1, a, 2 − a, x, y Here a can take any value except negative integers and positive integers greater than 2. We can further simplify the above relation by using the following relation of F 4 [30] xy . (5.5) For our case α = 1, β = a. Thus we get ( (x + y − 1) 2 − 4xy + x + y − 1) 2 4xy (5.7) As a consequence of this we get F 4 (1, 1; 1, 1; x, y) F 4 (1, 1; 1, 1; x, y) = 1 We can also consider the result for the bubble integral with general masses and unit power of propagators, for which the result is given as follows [18] With the help of the reduction procedure, the result for bubble integral is given by Eq. (2.13). This equality of Eq.(5.9) and (2.13) thus provides a reduction formula for the hypergeometric series in Eq.(5.9), which can be written as follows We can also obtain new hypergeometric relations by deriving other functional equations for the Feynman integrals(see appendix A). Using Eq.(A.3), (A.4) and and we get We further obtain We provide a list of various reduction formulae that can be derived using Eq.(5.4) and (5.12) in the appendix B. The right-hand side of Eq.(5.4) and (5.12) can further be equated to give the relation between the sum of hypergeometric 2 F 1 functions. An interesting consequence of this relation can be obtained with a = 3 As before, such a reduction also helps if the analytic continuation has to be performed to reach a certain kinematical region. We can find the analytic continuations for the series in Eq.(5.9) using automated tools [23], but it still does not guarantee that the parameter space has been covered. In contrast, the complete list of analytic continuations for the hypergeometric 2 F 1 [25] is available and well implemented in software like Mathematica . The complexity of the analytic continuation procedure also increases with the increase in the number of variables of the hypergeometric function due to the increase in difficulty to find the ROC of the resulting series.
We notice that the procedure is sufficiently general and one can obtain a large number of reduction formulae using it by doing the following steps • We take the Eq.(5.1) or any other general result for N −point integral from [17,18].
• For a N − point function we have a product of N − propagators. We take any two propagators and find the algebraic relation. This results in a sum of 2 terms, for which the number of variables in the result, as in Eq.(5.1), is reduced by one. This yields a reduction formula between, say L (which is a function of N ) variable hypergeometric function and (L − 1) variable hypergeometric function.
• We apply the previous step again, thus resulting in a relation between L variable hypergeometric function and (L − 2) variable hypergeometric function. Also using the previous step it gives a relation between (L − 1) variable hypergeometric function and (L − 2) variable hypergeometric function.
• We apply the procedure recursively until we have an algebraic relation for the product of N massive propagators as sum of 2 N −1 terms, such that each term contains product of N −propagators with just one massive propagator.
• The final result of the procedure would be a collection of relations between L, (L − 1), ...1 variable hypergeometric functions.

Summary and Discussion
We have presented an automatized package AlgRel.wl for finding the algebraic relation for the product of propagators. These relations were used by Tarasov [1,4,10,11] to derive the Functional relations for Feynman integrals. The results obtained using the package are also sufficiently general and can be used further to obtain the functional relations for the Feynman integrals by appropriately choosing the arbitrary parameters. In the present work we focused on automatizing the method to derive algebraic relation for the propagators by suitably implementing a recursive algorithm (a slight modification to the Tarasov's algorithm [1]). Furthermore, using a loop-by-loop approach we provided a systematic way so as to use these relation for higher loop integrals too. These relation occur with free parameters which can be chosen suitably.
Using various examples up to three-loops, we focused on how with a simple choice of these free parameters we can reduce integrals with large numbers of massive propagators into integrals with fewer massive propagators [10], which can thus be computed easily. For the one-loop case, we obtained results for up to 6point integral with the procedure and wrote them as a sum of 2 N −1 (for N − point integral) integrals with one massive propagator. We also showed how the procedure can be used for higher-loop integrals too where a loop-by-loop strategy has been applied for finding the relations.
Since the general results for the one-loop N −points integral are explicitly known for various cases in terms of multi-variable hypergeometric functions, we show how the present work can be used to obtain a large list of reduction formulae for these functions. As a demonstrative example of the same, we used the one-loop bubble integral where the reduction of the Appell F 4 to hypergeometric 2 F 1 can be obtained. We also derive another reduction formula for a 2-variable hypergeometric series, Eq.(5.9) in terms of hypergeometric 2 F 1 . The relations thus obtained, can be treated as general reduction formulae for these functions without making reference to the Feynman integral they were derived from. These relation hence provides a way to derive non-trivial reduction formulae for multi-variable hypergeometric function using physical problems. They are also helpful, especially for situations where the analytic continuation of multi-variable hypergeometric functions has to be obtained to evaluate them outside their ROC, which is not easy to derive otherwise.
The present procedure of finding algebraic relation for the product of propagators can be used only if the propagators are dependent on just one loop-momenta. For this reason, the procedure cannot be applied with full generality to multi-loop integrals and a loop-byloop approach has to be adopted. Hence the procedure is not helpful for integrals such as the sunset integral or in the cases where for each loop momenta k i there is just one propagator. To apply such a procedure to sunset-like integrals, a generalization of the procedure for the multi-variable case, when the propagators can depend on more than one loop momenta has to be developed.
As we have seen that the algebraic relation obtained reduces the complexity of the Feynman integral. Specifically for the simple case of one loop bubble (in Section 2), we saw that the result for general bubble integral, which was expressed in terms of double variable hypergeometric function Appell F 4 , was reduced to 2 F 1 which is a single variable hypergeometric function. It would be worth studying such reduction in complexity for other non-trivial cases of Feynman integrals which result in multi-variable hypergeometric functions of even higher variables. Since obtaining analytic expressions might not be feasible for such cases, a detailed numerical study for the same would be an important application of these algebraic relations after the proper function relations have been obtained by the proper choice of arbitrary variables. We would also like to point out to the possibility of using Lemma B.3., given in [31], to find the similar relations as presented here. 3 .

A Functional reduction with M i ̸ = 0
In this appendix, we point out other possibilities for the choice of arbitrary parameters M i [11]. This choice leads to different functional reduction equations than already presented. Also, this gives rise to different reduction formulae as has been done in section 5. Consider the bubble integral considered in section 3. This time we choose a different non-zero value of M 1 . Since the Feynman integrals are relatively easier to compute with equal masses a suitable choice is M 1 = m 1 . With this choice, we get, similar to Eq.(2.9), the following relation I 2 (p 2 , m 1 , m 2 ) = x 1 I 2 ((P 1 + p) 2 , m 1 , m 2 ) + x 2 I 2 (P F 4 (0, 1; 2, 0; x, y) = −2(x + 1)y + (x − 1) 2 + y 2 + x − y + 1 2x (x + y − 1) 2 4xy = −2(x + 1)y + (x − 1) 2 + y 2 − (x + y − 1) (1 − a)(x + y − 1) (x − 1) 2 + y 2 − 2(x + 1)y − (x + y − 1) We can use the following relation as given in [30] and obtain formulae for F 1 γ; x, xy) (B.9) We can further exploit the relation between F 1 and F 2 to derive reduction formulae for F 2 . Wherever possible exploiting such relations amongst various hypergeometric functions we can derive reduction formulas for other hypergeometric functions.

C Numerical results
In this appendix, we present the results of the numerical checks using FIESTA5 [26]. We use the AlgRel.wl package to obtain the algebraic relation and then perform the numerical integration of both the left-hand side, which is the original integral and the right-hand side which is the sum of integrals obtained using the algebraic relation. We give the results with 5 significant digits and O(ϵ) in the Laurent expansion where ϵ = 4−d 2 .
We find that the results obtained using the algebraic relations are numerically consistent.