Evaluating the last missing ingredient for the three-loop quark static potential by differential equations

We analytically evaluate the three-loop Feynman integral which was the last missing ingredient for the analytical evaluation of the three-loop quark static potential. To evaluate the integral we introduce an auxiliary parameter $y$, which corresponds to the residual energy in some of the HQET propagators. We construct a differential system for 109 master integrals depending on $y$ and fix boundary conditions from the asymptotic behaviour in the limit $y\to \infty$. The original integral is recovered from the limit $y\to 0$. To solve these linear differential equations we try to find an $\epsilon$-form of the differential system. Though this step appears to be, strictly speaking, not possible, we succeed to find an $\epsilon$-form of all irreducible diagonal blocks, which is sufficient for solving the differential system in terms of an $\epsilon$ expansion. We find a solution up to weight six in terms of multiple polylogarithms and obtain an analytical result for the required three-loop Feynman integral by taking the limit $y\to 0$. As a by-product, we obtain analytical results for some Feynman integrals typical for HQET.


Introduction
Three-loop corrections to the static potential between two heavy quarks 1 were evaluated in Refs. [2][3][4]. However, these results involved three constants which were obtained only in a numerical form. In order to find these constants in an analytical form one might try to apply the DRA approach formulated in Ref. [5]. An important step in this direction was made in Ref. [6]. And indeed, the analytical results for the two of the three missing constants were obtained withing this approach and presented in Ref. [1]. However the calculation of the last constant appeared to be much more involved for the DRA method, mostly, due to overwhelming numerical issues.
The goal of the present paper is to attack the problem from a different direction, by using the differential equations approach, along the lines of Ref. [7]. We successfully apply this strategy thus finding the last missing ingredient for a completely analytical result for the three-loop static quark potential. The result for the potential is presented in Ref. [1].
The three-loop Feynman integral which up to now was considered as the most complicated ingredient for the analytical evaluation of the three-loop quark static potential is with the usual −i0 implied in all propagators, v · q = 0, is a non-planar family of dimensionally regularized (with D = 4 − 2ǫ) Feynman integrals with twelve indices. As the dependence on both v 2 and q 2 is trivially factorized, we put v 2 = −q 2 = 1 in what follows. There is only one master integral in the sector with a i > 0, i = 1, 2, . . . , 9, 11, a 10 , a 12 ≤ 0 and it is more convenient [8] to choose it to be I 1 , rather than the corner integral F 1,...,1,0,1,0 of this sector. The integral I 1 = I 1 (ǫ) is finite in four dimensions, D = 4. The numerical value I 1 (0) ≈ −20.9484 was obtained in the calculation of Refs. [2,3] (see also [8]) using a Mellin-Barnes representation. In order to use the method of differential equations [9][10][11][12][13][14][15][16] one has to introduce one auxiliary parameter [7]. Though it might look like unnecessary complication, the idea behind this step is the following. An auxiliary parameter, when chosen judiciously, allows one to "deform" the original family of integrals to the point which is more accessible for other methods, in particular, to DRA method. Then, putting boundary conditions at this point and using powerful machinery of the method of differential equations, one can obtain results for the original integrals.
In the next section we explain how we introduce an auxiliary parameter and solve differential equations for the corresponding set of the master integrals. In Section 3 we describe how boundary conditions are fixed. As a by-product, we obtain analytical results for a family of Feynman integrals typical for the Heavy Quark Effective theory (HQET).

Master integrals and differential equations
In order to introduce an auxiliary parameter, one might think of relaxing the condition v · q = 0. However, this choice complicates the situation considerably because one effectively obtains three scales, or two variables in differential equations. Instead, we introduce parameter y by the following replacement and similarly for propagator number 11. This prescription corresponds to the introduction of the residual energy in HQET propagators. The factor 1/2 is introduced for convenience. We arrive at the following family of Feynman integrals Using FIRE [17][18][19] combined with LiteRed [20,21] we reveal 109 master integrals. The derivation of differential equations for a family of master integrals is a straightforward procedure. We take derivatives of the master integrals in y with the help of LiteRed [20,21] and then apply FIRE to reduce the resulting integrals to master integrals. As a result we obtain a system of linear differential equations where F is the column-vector of primary master integrals and A is a 109 × 109-matrix.
According to the strategy suggested in [15] and first applied in [7,22,23] and then in many other papers, in particular at the four-loop level [24], it is reasonable to try to pass to a new basis (canonical basis in terms of Ref. [15]) where the differential equations take the form with a factorized ǫ-dependence on the right-hand side (or ǫ-form for brevity).
To do this we use the algorithm [16] introduced by one of the author of the present paper. We successfully obtain a globally Fuchsian form and an ǫ-form for each diagonal block (which basically correspond to sectors). However, we fail to perform the last step (namely, the factorization of the ǫ-dependence) for the matrix as a whole. Therefore, we arrive at the form where A is a block-diagonal matrix independent of ǫ and B is the strictly lower-triangular matrix with zeros in each diagonal block ofÃ. Besides, they both have a global Fuchsian formÃ with y k running over the set {−1, −1/2, 0, 1, 1/2}. Since we want to obtain the boundary conditions from the limit y → ∞, we reduce the matrix residue at infinity, − k (ǫa k +b k (ǫ)), to a Jordan normal form, which, in particular, means that we secure that Although the resulting system (2.4) is not in an ǫ-form, this form is completely sufficient for our purposes. Indeed, for integrals of a certain sector, we have the system where the inhomogeneous term R(y, ǫ) is a linear combination of simpler master integrals, the coefficients being rational functions of y and ǫ. We assume that these simpler integrals have been already calculated at this point in a sense that we know their expansion in ǫ at fixed y up to a sufficiently high order as well as their asymptotic behavior at fixed ǫ and y → y k . Then we also know the expansion and the asymptotic behavior of R(y, ǫ), irrespectively on whether the coefficients of the linear combination are proportional to ǫ. We only note that, the condition (2.6) essentially simplifies finding the asymptotics of R(y, ǫ) at y → ∞.
To solve the linear system (2.4) in a power expansion in ǫ we pass to the variable x = 1/y and obtain results in terms of multiple polylogarithms of the argument x, up to expansion coefficients of 109 unknown functions C i (ǫ). The multiple polylogarithms [25] are defined recursively by G(a 1 , . . . , a n ; z) = z 0 dt t − a 1 G(a 2 , . . . , a n ; t) (2.8) with a i , z ∈ C and G(; z) = 1. In the special case where a i = 0 for all i one has by definition The letters a i for our result belong to the alphabet {−2, −1, 0, 1, 2}.

Fixing boundary conditions
To fix the unknown constants in the result, we consider the leading order asymptotic behaviour of the solution of Eqs. (2.4) in the limit y → ∞. Terms of the corresponding expansion can be described in the language of expansion by regions [26,27]. In our case, all the contributions of regions are obtained by considering each loop momentum to be hard (i.e. ∼ y) or soft (i.e. ∼ y 0 ). The crucial point in the matching is that the leading-order contributions are classified according to the power dependence on y. This parameter enters with powers of the form −kǫ where k ∈ {0, 2, 4, 6} and this is seen both from the point of view of differential equations and expansion by regions.
To fix the constants in our solutions of differential equations we evaluate, first, the leading asymptotic behaviour associated with the region where all three loop momenta are hard, i.e. ∼ y, All these contributions are described as integrals of the following family . (3.1) We imply that a 8 ≤ 0. These integrals are typical to the HQET. Using FIRE we reveal 9 master integrals in this family, The first 5 master integrals can easily be evaluated by a consecutive integration over loop momenta and expressed in terms of gamma functions at general ǫ. To evaluate the four more complicated integrals we apply the DRA method [5] and use the dedicated package SummerTime [28] for the calculation of multiple sums appearing underway.
In addition to the information about this region, we use the leading order behaviour for the region where all three loop momenta are soft, i.e. ∼ y 0 . The corresponding contributions are either zero or expressed in terms of well-known massless propagator integrals. It turns out that after this it is sufficient to add information about 'intermediate' contributions for a small number of master integrals to fix all the constants.
Upon fixing all the constants in our solution we obtain analytic results for all the 109 elements of the canonical basis from which we derive analytical results for the primary basis. We checked our results numerically with FIESTA [29].
The integral I 1 = F 1,...,1,−2,1,0 at ǫ = 0 can be obtained as the value at y = 0 of the corresponding integralĪ 1 (y) depending on y and belonging to the family (2.2). Since we obtained our results for the 109 master integrals in terms of multiple polylogarithms with the variable x = 1/y we obtain also an expression forĪ 1 (1/x) in a similar form. We used two ways to analytically evaluate the limiting values at x → ∞. First, we evaluateĪ 1 (1/x) at the point x = 10 200 with the computer implementation [30] of GiNaC [31] with the precision of 200 digits and this happens to be enough to arrive at an analytical result using the PSLQ algorithm [32]. Second, we rewrite the result forĪ 1 (y) in terms of multiple polylogarithms of the variable y using our implementation of an algorithm presented in Ref. [31], then set y = 0 and again apply PSLQ. This time, increasing precision can be done much easier. We successfully checked the result obtained in the first way with 500 digits.
As a result we obtain the following value of the missing constant I 1 (0) for the three-loop quark static potential which corresponds to I 18 in the notation of Ref. [1].

Conclusions
Using the method of differential equations we analytically evaluated the three-loop Feynman integral which was the last missing ingredient for the analytical evaluation of the threeloop quark static potential. In the accompanying paper [1] this result is used to present completely analytical results for the three-loop quark static potential. As a by-product, we obtained analytical results for some Feynman integrals typical for the Heavy Quark Effective theory. Our results for the master integrals involved are available 2 in a computer-readable format at http://theory.sinp.msu.ru/~smirnov/tlp.