Four-dimensional differential equations for the leading divergences of dimensionally-regulated loop integrals

We invent an automated method for computing the divergent part of Feynman integrals in dimensional regularization. Our method exploits simplifications from four-dimensional integration-by-parts identities. Leveraging algorithms from the literature, we show how to find simple differential equations for the divergent part of Feynman integrals. We illustrate the method by an application to heavy quark effective theory at three loops.


Introduction
Our ability to evaluate Feynman integrals is crucial for computing perturbative results in quantum field theory. This subject has a long and rich history of various methods being developed [1]. In recent years, the method of differential equations [2][3][4][5] has become the main method for evaluating loop integrals that are needed for collider phenomenology. In particular, insights about writing the differential equations in a canonical form [6], has allowed both for an unprecedented degree of automation of the calculations, and has significantly pushed the boundaries of what is achievable.
Calculations are typically done within dimensional regularization, with D = 4 − 2ϵ, and the results are expanded in a Laurent series in the dimensional regulator ϵ. One advantage of the differential equations method is at the same time a drawback: in principle, the differential equations describe the functions under consideration at any order in ϵ. In practice, however, one often wishes to know results up the the finite part only, so that a lot of unnecessary information is kept at intermediate steps. What is more, the higher-order terms in ϵ are typically scheme-dependent, which means that they are more complicated than the physically relevant terms. Therefore it is desirable to develop methods that directly compute the physically relevant information, and profit from the simplifications. Such a method exists for the case of (infrared-and ultraviolet-)finite Feynman integrals [7]. The authors showed that it is possible to directly compute integration-by-parts-relations (IBP) in four dimensions, taking however special care about possible contact terms that relate integrals of different loop orders. There are several advantages to this. Firstly, the restriction of finite integrals means that significantly fewer integrals and IBP relations between them need to be considered, and the resulting number of master integrals is also less than in the conventional D-dimensional case. The simplifications in terms of handling the IBP relations in [7] are so substantial that even certain three-loop four-point integrals with massive internal lines could be computed that were beyond the reach of conventional techniques. Secondly, when truncating to four dimensions, the canonical differential equations matrix takes a simple, 'block-diagonal' structure, which reflects transcendental weight of the different Feynman integrals. A similar decoupling of the differential equations was also noticed in [8]. This makes it possible to iteratively solve for the master integrals, with increasing transcendental weight.
In many practical situations however, one often has to deal with divergent Feynman integrals. In the context of dimensional regularization, one is really interested in the coefficient functions appearing in the Laurent series in the dimensional regulator ϵ. In practice, only a finite number of terms is needed. Therefore it would be desirable to have a method that directly computes those finite coefficient functions. Ideally, in order to profit from the advances mentioned above, this new method should be compatible with the canonical differential equations approach. Our aim is to develop such a method.
As a first step in this direction, in this paper we focus on Feynman integrals that are free from subdivergences, i.e. that in dimensional regularization only diverge as 1/ϵ, and we develop a four-dimensional integration-by-parts and canonical differential equation method for the coefficient function. We leave the generalization to the general case for future work. One might initially think that most Feynman diagrams have subdivergences, which limits the scope of the present paper. However, the structure of both ultraviolet (UV) and infrared divergences (IR) is very well understood, e.g. in terms of the famous Hopf algebra of renormalization [9], and in terms of factorization properties (see [10][11][12] and references therein). Therefore, at least in principle, one could use our method on Feynman integrals in those cases as well, for some suitably-defined subtracted integrals (see e.g. [13] for recent work on an infrared subtraction for Feynman integrals).
We find it useful to focus on a particular type of divergences to present the new method. In this paper we focus on integrals in heavy quark effective theory (HQET). These integrals arise for example when massive particles emit soft radiation and capture infrared divergences associated to this radiation. Equivalently, they can be used to compute ultraviolet divergences of certain Wilson lines with cusp (related to the angle-dependent cusp anomalous dimension).
Working with these integrals as a case in point, we explain the main ingredients of our method. We first define HQET integrals that are free of subdivergences, and specify their leading divergence. We call the latter I(ϕ). Our goal is to compute I(ϕ). The main advantage of focusing on this coefficient, as opposed to the full Laurent expansion is that we can disregard unnecessary information. We explain how one obtains four-dimensional integration-by-parts (IBP) identities valid for I(ϕ). The master integrals appearing in these relations are all free of subdivergences of themselves, and there are less of them compared to the generic case. In this work, we develop a new syzygy IBP method to forbid integrals with soft divergence. (The original syzygy IBP method was developed to reduce an IBP system's size [14,15].) Moreover, we find, as in [7], that there are integral identities that related different loop orders. As a result, the lower-loop information can be recycled. Using these simplified IBP relations, we obtain differential equations for the master integrals. Finally, we leverage ideas [16] of how to transform these equations into a simple canonical form.
This paper is organized as follows: in the section 2, we briefly review the Feynman integrals from web diagrams in HQET. In the section 3, our method of generating IBPs for the divergent part of Feynman integrals, based on graded IBP operators and syzygy, is introduced. In the section 4, we invent a four-dimensional version of the initial algorithm [16], to generate canonical differential equations without ϵ. In section 5, we provide a threeloop application of our method. Then we provide the summary and the outlook of our new method in the section 6.

From dimensionally-regularized Feynman diagrams to finite functions
Let us define the class of Feynman integrals that we study in this paper. We focus on integrals relevant for describing soft divergences occurring in the scattering of massive particles. These divergences can be described by an eikonal approximation, which leads to Wilson line correlators in position space. Equivalently, in momentum space, one obtains heavy-quark-effective theory (HQET) integrals [17].  Consider a soft exchange in the scattering of two massive particles. The leading divergence for such a process is given by the eikonal diagram in Fig. 1(a). It depends on the scattering angle ϕ, and is given by

Cusped Wilson lines in position space
(2.1) In the following we set v 2 1 = v 2 2 = 1 without loss of generality. The coefficient of the 1/ϵ pole is the angle-dependent cusp anomalous dimension [18][19][20]. Our goal is to compute the coefficient for this and similar higher-loop diagrams directly, using four-dimensional methods, without having to deal with the higher-order terms in ϵ.
When extracting information from divergent Feynman integrals it is paramount to consider carefully whether the computational steps are consistent and well-defined. In principle, there are several sources of divergences in Feynman integrals, in particular ultraviolet (e.g. from renormalization of the Lagrangian or of composite operators) and infrared (soft and collinear divergences associated to massless particles). Here we will focus on Feynman integrals that have one type of divergence only, leaving the case of integrals with divergences in several simultaneous regions for future work.
The knowledge of what region of loop integration produces the divergences allows us to extract its coefficient in a well-defined way. Let us now illustrate this for the case of the one-loop integral shown in Fig. 1, and then explain how we approach the problem in general. In position space, the integral of Fig. 1(a) is given by Here e −τ 1 −τ 2 is a cutoff that makes sure that the integral gives only the divergence at τ 1 , τ 2 ∼ 0 that is of interest to us. 1 We can make the latter manifest by an overall rescaling,

HQET integrals in momentum space
One can equivalently write the integral discussed in subsection 2.1 in HQET language in momentum space. We briefly recall some notations of HQET (see [21] for a detailed review).
HQET is an effective theory describing the dynamics of heavy quarks. In the heavy mass limit, the heavy quark's propagator becomes In fact, we are extracting the UV divergence of this integral. This may seem strange as our starting point was a soft exchange, i.e. an IR effect. The explanation is that these UV and IR divergence are related ΓUV = −ΓIR, via the formal relation where p = mv + k. v is the classical velocity of the heavy quark, normalized as v 2 = 1, while k is its off-shell momentum. In this limit, the mass parameter factorizes out and the quantum field theory computation is simplified significantly. In practice, we may use an infrared regulator δ for the linear propagator in (2.4). Since δ is the only scale for the HQET integrals, it is safe to set δ = 1. 2 In HQET language, diagram Fig. 1(a) is given by (up to an overall factor, D = 4 − 2ϵ) .
In eq. (2.5), the divergences comes from an overall rescaling where k µ ∼ ∞. Now, in the above example, one could directly compute the coefficient of the 1/ϵ pole by setting ϵ = 0 in the second factor of eq. (2.3). Something similar can be done in momentum space. We prefer however to develop a formalism that can be applied directly for covariant, four-dimensional Feynman integrals. The reason for this is that in this way we obtain a more general setup that we expect may be used for larger classes of Feynman integrals. For this reason we define Our goal is then to derive four-dimensional IBP relations that are valid for J.

Higher-loop diagrams free of subdivergences
Let us discuss the generalization to higher loops. Again, we focus on Feynman integrals with the divergence coming from one region of loop integration. In the case of our HQET integrals, we assume that they are free of subdivergences, and that the overall divergence is associated to a region where all loop momenta become large. See Fig. 1(b) for a two-loop example. Fig. 1(b) is an example of a so-called web diagram [22,23]. By definition, webs do not have one-particle irreducible subdiagrams, which ensures the absence of subdivergences. Web diagrams occur naturally in the study of Wilson lines in the context of non-Abelian eikonal exponentiation [22,23]. In a nutshell, they provide a setup for computing the cusp anomalous dimension from Feynman diagrams without subdivergences. So they are ideal objects to study for our purposes. To extract the overall divergence, we define in analogy to eq. (2.6), 3 In this paper we focus on HQET web diagrams built from two eikonal lines (corresponding to the heavy quarks), and an arbitrary number of massless exchanges. We impose the following conditions on them: 1. The integrand has the overall scaling dimension −4L, for the transformation k i → λk i , i = 1, . . . , L. This follows from eikonal Feynman rules. 2. No UV divergence in any subloop subdiagram.

No infrared (IR) divergences. This forbids higher powers of massless propagators.
We call such integrals "admissible". It follows from these properties that in dimensional regularization, admissible integrals have an overall ϵ −1 divergence only. We are interested in the coefficient of the divergence in (2.7). The three conditions formulated in subsection 2.3 can be translated to constraints on the Feynman integrals' indices. While this restricts the space of integrals of interest to us, in general the solution space is infinite dimensional. In practice it is sufficient to work with a finite-dimensional subspace. We do this by placing further constraints on degrees of numerators or denominators. Then we can find all admissible integrals are subject to these additional conditions with the help of computer algebra.

Finding admissible integrals
Let us illustrate this for the diagram shown in Fig. 1(b). We define the crossed-ladder integral family T (2) [a 1 , · · · a 7 ] as follows.
We formulate the following additional criterion (based on experience) for these two-loop integrals: if the number of propagators is p, then we allow up to 7 − p doubled propagators. For example, for six propagators, we allow one doubled propagator; for five propagators, we allow two doubled propagators, and so on. Of course, we could later relax this condition, if necessary. However we will see below that they are sufficient for our purposes. We then find the following admissible integrals (are subject to the additional conditions).
In top sector with six propagators (i.e. indices a 1 , . . . a 6 positive), we have the scalar crossed ladder diagram (see Fig. 2) (2.10) Moreover, there are two additional admissible integrals that involve numerators. The first one is given by the scalar crossed ladder integrand, but with additional factor inserted. Because of eq. (2.11), this integral (and a similarly constructed one) can be written as a difference of integrals in subsectors, as follows In order to make this manifest, we introduce the auxiliary noation . (2.14) where we define two reducible scalar products (RSPs) D 8 = −2k 1 · v 2 , D 9 = −2k 2 · v 1 . In this way, admissible integrals in (2.10) can be neatly written as, A final comment is in order. There are further integrals related to the above by symmetry relations (e.g. the v 1 ↔ v 2 symmetry implies a flip symmetry of graphs), which we do not show. One could generate these additional integrals, and then eliminate them automatically with the help of the integration-by-parts relations discussed in the next section. For sectors with five propagators, we find  For sectors with four propagators, we find These integrals are shown in Fig. 2. Since the integrals above are all admissible, the limit of eq. (2.7), that we are mostly interested in, is well defined. From now on we will use the same notation both for the integral and its value in the limit of eq. (2.7), hoping that this does not lead to confusion.

Four-dimensional integral reduction identities
Our next goal is to find linear relations between the ϵ −1 parts of admissible integrals. We denote V as the linear space spanned by admissible integrals. We then find a basis of this space, taking into account the following additional relations: 1. Apply the zero sector condition [24,25] to V to remove vanishing integrals, and also the Pak algorithm [26] to identify equivalent sectors. This reduces the number of admissible integrals that need to be considered. 2. Find IBP operators that generate relations between admissible integrals. Firstly, the operators need to be graded, which means that all terms have the same scaling dimension w.r.t. the loop momenta. This is necessary in order to be able to apply the UV conditions number 1 and 2 in subsection 2.3, for admissible integrals. Secondly, in order to fullfil the condition number 3 in subsection 2.3, the operators need to satisfy certain syzygy conditions like those in [14,15]. 3. For L-loop admissible integrals with bubble sub-diagrams, use a reduction identity for integrating out the bubble, leading to an (L − 1)-loop admissible integral.

Graded IBP operators and the syzygy method
A generic IBP relation has the schematic form, with vectors u i . However, this generic form of IBP relations contains non-admissible integrals in general. The reason is that the vectors u i may change the power counting of sub-loop diagrams and thus can introduce sub-loop UV divergences. Moreover, when ∂/∂k µ i acts on gluon propagators, it increases the propagator power and may worsen the infrared properties. Our goal is to generate valid IBP relations between admissible integrals. (Below we also use a similar strategy for finding differential equations between admissible integrals.) In order to achieve this, we develop a new way of generating IBP relations with the following constraints: 1. The IBP operator in (3.2) should be graded. This means that under the scaling of each loop momenta, k i → λk i , i = 1, . . . , L, The tuple β i , i = 1, . . . , L are the scaling dimensions for each individual loop momentum. An overall scaling dimension condition for IBP vectors has been used for the linear algebra algorithm of IBP relations without doubled propagators [15]. Here we use individual scaling dimensions to impose a well-defined power counting on the IBP differential operator for each sub loop. When deriving IBP relations we will use this power counting condition to ensure the absence of subdivergences. In practice, it is easy to make an ansatz for graded IBP operators with given a scaling dimensions β i , i = 1, . . . , L. 2. The IBP operator in (3.2) should not produce higher power of gluon propagators. Let D k be a gluon propagator, then, the requirement reads, where g k is a polynomial in the scalar products involving loop and external momenta. This is the syzygy relation for IBP relations [14]. The original purpose of using syzygy relations is to reduce the number of IBP relations for the reduction step. See [15,[27][28][29][30] for the development of the syzygy method for IBP reduction. In this paper, (3.4) is imposed to control the soft divergence and it is only on massless (gluon) propagators. Note that although this condition is not imposed on every propagator, it still significantly reduces the number of IBP relations.
To find suitable IBP differential operators, we use these conditions together. It is natural to make an ansatz for graded differential operators (first condition) and then solving the syzygy relation (3.4) with the linear algebra algorithm in [15]. Given an ansatz for the graded IBP operator, (3.4) becomes a sparse linear algebra problem. In practice, we wrote a proof-of-concept code in Mathematica to find such differential operators, which we make available in the auxiliary file "demo/graded vector demo.wl". There is also an available package NeatIBP [31] to generate small-size IBP relations by syzygy method.
Then we apply such differential operators on Feynman integrals with suitable indices, to get relations between admissible integrals. It is notable that not every term in such a relation is admissible, but suitable combinations of terms are admissible. In practice, to find such combinations, we first set up an integrand basis of admissible integrals, and then expand the relation over this integrand basis. This step is powered by a finite-field linear algebra code, which depends on SpaSM [32]. We provide this code in the auxiliary file "demo/3lHQET graded IBP demo.wl". After this, we have only admissible integrals in the relation. Therefore, if the coefficients of admissible integrals contain linear functions of ϵ, it is safe to set ϵ → 0. The final relation is thus independent of ϵ.
In the subsection 3.3, we illustrate how this strategy is applied for finding relevant IBP relations between the admissible integrals of our two-loop HQET example. Another useful technique for generating linear relations for admissible integrals, is to integrate out bubble sub-diagrams. To be specific, for an L-loop HQET integral with the loop momenta k 1 only appearing in a bubble, with a 1 > 2,

Integrating out bubbles: cross loop order reduction
Note that the new propagator's index depends on ϵ which is not always desirable. However, for the ϵ −1 part only, as we show presently, the following replacement is valid, for a 1 > 2, (3.6) The factor (L − 1)/L requires a comment. It takes into account the fact that the overall divergence we consider is computed from a volume integral, which in turn depends on the dimensionality, and hence the loop order. The factor compensates for the mismatch when comparing overall divergences of L and (L − 1)-loop integrals in dimensional regularization.
Applying this relation to the case of our family of crossed ladder integrals, we find the relations shown in Fig. 3.

Application to crossed ladder integrals
We construct IBP vectors, or equivalently differential operators as in eq. (3.2), which do not double propagators number five and six (see Fig. 2 and admissibility condition number 2 in subsection 2.3). We find that it suffices to consider those with the lowest scaling dimension in {k 1 , k 2 }, In the top sector, the IBP relations, P 2 T (2) [1, 1, 1, 1, 1, 1, 0] = 0 , P 3 T (2) [1, 1, 1, 1, 1, 1, 0] = 0, (3.10) together with symmetry relations yield As a result, we keep only one master integral in the six-propagator sector, namely A 1 , as A 2 and A 3 can be reduced to lower-sector integrals. It is possible to verify that finite integrals with two or more dots can be reduced in similar ways by applying the same set of IBP vectors on finite seeds which have more dots but fixed scaling dimension.
In summary, we find that the admissible crossed ladder integrals can be reduced to the following basis integrals: (3.20) In the next section, we show how to compute these integrals from differential equations.

Algorithm for four-dimensional canonical differential equations
Having identified the IBP relations between admissible integrals, and hence a set of master integrals, we can now derive differential equations for the latter. The HQET integrals we consider depend on a single variable ϕ. It is convenient to trade the latter for another variable x = exp(iϕ) We may differentiate the master integrals with respect to x by means of a differential operator acting on the integrand, e.g.
Note that this operator does not double any gluon propagators and preserves the scaling dimension in each loop momenta. Thus for any admissible integral f ∈ V, df ∈ V. Given a basis of master integrals ⃗ f , one can write down a linear system of differential equations of degree N = dim(V) : Let us illustrate how this works in the case of the crossed double ladder A 1 . The crossed ladder integral family has a basis of five admissible master integrals Let us now derive the differential equation that A 1 satisfies, and then determine it from now. Taking the derivative of A 1 by means of the differential operator given in (4.1), we have Taking into account eq.(3.11), this becomes One may iterate this procedure for B 1 , B 2 , and then for C 1 , C 4 . One finds that the integrals satisfy the differential equations (4.2) with In order to solve this system of differential equations, it is useful to first simplify the equations, following [6,7] (see also [33]), or using the Initial algorithm [16]. The system above is of course very simple, so many different methods could be used. The key idea proposed in [6] is to convert the DE to a canonical form. The latter is obtained by choosing the basis integrals astutely. In fact, in order to find the transformation to a canonical form, under certain assumptions it is sufficient to choose one integral in a good way [16]. Here we can choose This choice can be justified by an integrand analysis [6], as explained in detail in [17]. The factor ensures that F 1 has constant leading singularities. Applying the Initial algorithm, we find that F 1 satisfies a fourth order Picard-Fuchs equation, which defines a rank-four canonical differential system where We see that eq. (4.8) is in canonical form. Moreover, the matrix is block triangular [7]. This means that we can iteratively solve this equation. At each integration, there is one boundary constant to be fixed. The first integral, F 4 , is constant, and we find its value by direct computation. The other boundary constants are fixed by observing that F 1 , F 2 , F 3 vanish at x = 1. In this way, we find This is in agreement with the known answer for A 1 , see e.g. eq. (27) in [34]. This determines four of the five admissible master integrals. In case one wished to solve for the remaining one, one can do so by repeating the above steps with F 5 = 1 2 (x − 1/x) 2 B 2 as the input to the Initial algorithm. In this case, one finds (4.11) Solving this system, we find Hence we have found the complete basis of admissible UT integrals in the crossed ladder family.

Three-loop application
Here we apply our new four-dimensional method to three-loop integrals initially computed in [17,35]. We find that in comparison, much fewer master integrals are needed, and the computation requires only a fraction of the computing time. For example, the standard differential equation approach for our three-loop HQET example via a publicly-available multiple-thread IBP solver took several hours, while our approach only took minutes with one core on the same machine. The propagators of the three-loop integrals in [17] are, For example, say, our goal is to calculate the ϵ −1 part of G[1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1], which is an admissible integral (see Fig. 4). To set up a four-dimensional differential equation system, we first list four-dimensions IBPs for admissible integrals, as described in section 3, with a bound on the propagator indices. Besides the four-dimensional IBP relations, we also use the cross loop order relations introduced in the section 3. We find that some specific integrals in two sectors can be reduced to lower loop integrals. Combining the four-dimensional IBPs and cross loop-order relations, we get a linear system for the ϵ −1 part of the admissible integrals. This linear system is very sparse, free of the parameter ϵ, and can be solved by FiniteFlow [36,37] in two minutes with one CPU core. We find 13 irreducible integrals (cf. the auxiliary file "output/3lHQET 4D MI.txt") that satisfy the differential equation, where A is a 13 × 13 matrix. The explicit expression of A is given in the auxiliary file "output/3lHQET 4D DE.txt".
A straightforward integrand analysis, as explained in detail in [17], suggests that the integral Solving the differential equations and fixing the boundary constants, we find that there is degeneracy, and hence the dimension of the finite system is reduced to six. The solution for the top integral J 1 can be easily found as, We find the following extensions of the work in this paper promising: 1. Apply the method to cutting-edge applications that go beyond the state of the art.
Potential applications include anomalous dimensions of composite operators, and the soft anomalous dimension matrix [38]. 2. Extend the method to Feynman integrals with subdivergences, truncating the Laurent expansion in the dimensional regulator at a given order. This requires a careful identification of the relevant regions of loop integration. This will likely lead to an increase in the number of finite master integrals to be considered, but we still expect a significant simplification compared to the general, D-dimensional case. 3. Combine the method with insights into the structure of the Feynman graph polynomials from tropical geometry. The latter allow to identify the relevant divergent regions and to write their (finite) coefficients in terms of integrals that appear naturally from the geometry [39,40].
also acknowledges the Institute of Theoretical Physics, Chinese Academy of Sciences, for the hospitality through the "Peng Huanwu visiting professor program".