Analytic results for planar three-loop integrals for massive form factors

We use the method of differential equations to analytically evaluate all planar three-loop Feynman integrals relevant for form factor calculations involving massive particles. Our results for ninety master integrals at general $q^2$ are expressed in terms of multiple polylogarithms, and results for fiftyone master integrals at the threshold $q^2=4m^2$ are expressed in terms of multiple polylogarithms of argument one, with indices equal to zero or to a sixth root of unity.


Introduction
Form factors of heavy quarks are important in top quark physics, see e.g. [1][2][3]. While the form factors are interesting at general momentum transfer, a particularly important kinematical regime is the threshold expansion. One of the main technical obstacles in studying this regime is the lack of higher-order analytic results for the threshold expansion of Feynman integrals contributing to it.
We consider the following family of planar vertex Feynman integrals F a 1 ,...,a 12 = 1 (iπ D/2 ) 3 . The parametrization of eq. (1.1) can be used to describe any planar Feynman integral of this type. It is based on dual or region coordinates. When considering Feynman integrals with less than 9 propagators, it typically happens that their graph can be represented as a subdiagram of more than one of the families shown in Fig. 1. Different representations can most easily seen to be equivalent by a permutation of the (dual) integration variables. The integrals also have a flip symmetry, since the integrated result only depends on p 1 , p 2 through q 2 . In this way, a given diagram can be represented in many equivalent ways. One may also use computer programs such as [42] to find such equivalences.
The outline of the paper is as follows. In the next section, we explain how we evaluate MI for integrals (1.1) at general q 2 and, in Section 3, we obtain analytical results for MI for integrals (1.1) at q 2 = 4m 2 using these general results and matching at threshold. For convenience, we provide the main results, as well as further key information, in terms of ancillary files. Appendix A contains a pedagogical one-loop example of all steps of the calculation that can be followed in detail.
2 Integrals at general q 2 : solving DE To solve integration by parts (IBP) relations [38] using FIRE [39][40][41] combined with LiteRed [42] we reveal 90 MI at general q 2 while the family of threshold MI has 51 MI.
Suppose that we are evaluating MI for a given family of Feynman integrals. Let us denote the kinematical variables by x = (x 1 , . . . , x n ), the set of MI by f = (f 1 , . . . , f N ), and let us work in D = 4 − 2ǫ dimensions. The general set of DE takes the form where ∂ i = ∂ ∂x i , and each A i is an N × N matrix. In [4], it was suggested to turn to a new basis of the master integrals having constant leading singularities, for which the DE should take the following form One essential difference with respect to (2.1) is that the matrix in this equation is just proportional to ǫ.
In the differential form, we have The matricesÃ α k are constant matrices and the arguments of the logarithms α i (letters) are functions of x. Let us deal with the case of two scales, i.e. n = 1 so that x is just one variable. Then the desirable form of the DE is where x (k) is a set of singular points of the DE, and the N × N matrices a k are independent of x and ǫ. For integrals (1.1) considered at general q 2 , let us introduce the variable Note the x ↔ 1/x symmetry of this definition. The values x = 0, x = −1 and x = 1 correspond to the high energy limit q 2 = ∞ (or m 2 = 0), the threshold limit q 2 = 4m 2 , and to the soft limit s = 0, respectively. Below, the latter limit is used as a boundary point when solving differential equations. To convert DE for this family of integrals into the form (2.3) we follow the strategy of [4,29]. The main point of the method is to choose integrals having constant leading singularities. (Sometimes we also used small additional basis transformations to 'integrate out' unwanted terms.) In this way, we obtain where A i , i = 1, 2, 3, 4 are constant (x-and ǫ-independent) matrices. Our basis choice f , as well as the corresponding differential equation matrix on the r.h.s. of eq. (2.7), is given in an ancillary file, for convenience of the reader. We see that eq. (2.7) takes the form of eq. (2.3), with the letters are x, 1 + x, 1 − x, 1 − x + x 2 . This is very interesting, since for analogous integrals up to two loops, only the letters x, 1 + x, 1 − x appeared [1,27].
Those letters, and the corresponding singularities have a clear physical interpretation, in terms of threshold, soft, and high-energy limits. The presence of the letter 1 − x + x 2 is a new feature at three loops. In terms of the original variables, it corresponds to a pseudo-threshold at q 2 = m 2 . We stress that we did not have to guess the presence of this letter, but rather found it systematically by setting up the differential equations. It would be interesting to derive the presence of this (spurious) singularity from an analysis of the Landau equations. See ref. [25] for recent work in this direction.
The differential equations are most straightforwardly solved directly from (2.7), in terms of iterated integrals. This also gives the shortest and most flexible representation of the answer. As an example, we present the first term in the ǫ expansion for one of the basis integrals, The graph is shown in Fig. 2. This is also one of the simplest examples where the new letter 1 − x + x 2 makes an appearance. The most natural way of writing the answer is as Chen iterated integrals [24] with boundary point x = 1. We use brackets to denote the latter, e.g.
These integrals have many nice properties. See e.g. [26] for more examples and applications. Using this notation, we have (2.11) To be more explicit, let us carry out the first three integrations. This gives In order to make contact with more commonly used classed of functions, we also give the solution in another form. In order to do this, we first rewrite the differential equation (2.7) in the form (2.5), at the cost of introducing complex roots of the polynomial Interestingly, the latter are 6th roots of unity. In this form, the DE admits a natural solution in an ǫ expansion with coefficients written in terms of Goncharov (multiple) polylogarithms (GPL) [37]. The latter are defined recursively by with a i , z ∈ C and G(z) = 1. In the special case where a i = 0 for all i one has by definition (2.14) Given the alphabet, the a i can take the values 0, ±1, r 3,4 . We find where The main differences to eq. (2.11) are the following: the letter 1 − x + x 2 was factored into linear pieces, the Goncharov polylogarithms have x = 0 as boundary point, and (by convention), the indices are read in the opposite order, compared to the [. . .] notation. We derived analytic results for all integrals up to transcendental weight six. This is expected to be sufficient for three-loop computations. If needed, higher-order terms in the ǫ expansion can be obtained by further expanding. Our analytical results are presented in ancillary files.
3 Integrals at q 2 = 4m 2 : matching at threshold In this section, we explain the results at general q 2 can be used to extract the integrals at threshold q 2 = 4m 2 . A subtle point is that near threshold, different scaling regions contribute. Crucially, the exact knowledge of the differential equations near threshold allows us to properly disentangle those contributions, as explained in this section.
Let us turn to the integrals considered at q 2 = 4m 2 . We reveal a set of 51 master integrals.
The idea is to obtain analytical results for the threshold MI using our results at general q 2 , via the threshold expansion. For a given three-loop Feynman integral at general q 2 , the latter has the form where the summation over n is over integer or half-integer numbers. According to the strategy of expansion by regions [43,44], the threshold expansion is given by a sum over so-called regions where every loop momentum can be of the following four types: hard, potential, soft and ultrasoft. At each loop, a potential and a soft loop momentum gives −ǫ to the exponent of the expansion parameter in (3.1) and an ultrasoft loop momentum gives −2ǫ.
Our goal is to compute the MI of the family of the 'naive' (hard) values at threshold. They correspond to the one-scale integrals F 0,0 (a 1 , . . . , a 12 ; 4m 2 ) defined with q 2 set to 4m 2 , i.e. under integral sign, either in integrals over loop momenta or in Feynman parametric integrals. Such integrals correspond to the contribution of the region where all the loop momenta are hard.
Unfortunately, we cannot just set q 2 = 4m 2 , i.e. x = −1 in our basis because some integrals enter with the coefficients 1/(x+1) and 1/(x+1) 2 . These are spurious singularities which eventually cancel between different terms in the definition of the basis integrals. It might be possible to solve this practical problem by choosing a basis without such spurious singularities. Here, instead, we chose to 'naively' expand in x + 1 some of the Feynman integrals involved in the basis at least up to the second order. In order to do this, we have to deal not only with threshold integrals but also with their ('naive') derivatives. We introduce one more (13th) index for the order of this derivative in q 2 , i.e. we want to deal with the family F ′ (a 1 , . . . , a 12 , a 13 ; q 2 , m 2 ) = ∂ ∂q 2 −a 13 F (a 1 , . . . , a 12 ; q 2 , m 2 ) where the derivative is understood in the naive sense.
Taking naive derivatives at threshold can be illustrated by the simple example of the triangle diagram of Fig. 3, with p 2 1 = p 2 2 = m 2 i.e. the one-loop prototype of our three-loop diagrams. It is given by The corresponding Feynman parametric integral is Its naive derivative is obtained by differentiating in q 2 and setting q 2 = 4m 2 under integral sign. In this simple case, the general term of the corresponding naive series in q 2 − 4m 2 can be evaluated in terms of gamma functions at general ǫ, with the result Of course, the values of naive derivatives at threshold are not derivatives of the full integral.
In particular, the latter would include also the singularity (4m 2 − q 2 ) −1/2−ǫ corresponding to the potential region. Writing down IBP relations for integrals at general q and expanding all the terms naively in q 2 at q 2 = 4m 2 gives fifteen IBP relations for integrals (3.2) with thirteen indices. Following [50] we introduce one more relation which is obtained from (3.1) by a naive differentiation in s. As a result we obtain the possibility to express, by solving these IBP relations, any F ′ (a 1 , . . . , a 12 , a 13 ) in terms of master integrals. To do this, we use FIRE and observe that the corresponding master integrals are all with a 13 = 0, i.e they directly correspond to the 51 MI of the family of the threshold integrals which are the goal of our calculation in this section.
To match our analytic results for the 90 MI at general q 2 and arrive at analytic results for the 51 threshold MI, we analyze our DE (2.5) at the singular point x = −1. Let us change the variable x = y − 1 so that now we are interested in the behaviour of our DE at y = 0. Near this point the DE (2.5) has the form whereÃ ′ (y) = A 0 + yA 1 + y 2 A 2 + . . .. It turns out that the language of DE provides an alternative description of the threshold expansion (3.1): the eigenvalues of the matrix A 0 correspond to contributions of various regions within expansion by regions [43]. In the language of DE, the naive part of the expansion near y = 0 corresponds to the zero eigenvalues of the matrix A 0 , while eigenvalues of the form −kǫ with positive integer k correspond to the other contributions. Sometimes, we also need power suppressed terms in this expansion. For this, we use a trick from the theory of DE (see, e.g., [36]). One is looking for a polynomial P = 1 + r=1 P r y r such that the DE for the function g defined by f = P g takes the form yg ′ (y) = ǫA 0 g(y) where A 0 is independent of y. Then the solution of this equation is just g = y ǫA 0 g 0 with a boundary value g 0 . The full solution is then given by f (ǫ, y) = (1 + r=1 P r y r )y ǫA 0 g 0 . (3.7) We implemented the algorithm presented in [36] and easily obtained polynomials P r at least up to r = 10.
We perform matching at threshold in the following way: having determined both (3.7), and the full solution f (ǫ, y), we compared the ǫ → 0 expansion of the former to the threshold expansion of the latter. In this way, we identified g 0 . The knowledge of y ǫA 0 g 0 then allowed us to match with F 0,0 (a 1 , . . . , a 12 ; m 2 ) in eq. (3.1).
Solving these equations we obtain coefficients of the ǫ expansion of the MI up to some order written in terms of GPL G(a 1 , . . . , a n ; 1) with a 1 = 1 and a i taken from the sevenletters alphabet {0, r 1 , r 3 , −1, r 4 , r 2 , 1}.
The numbers G(a 1 , . . . , a n ; 1) form an important set of constants which appear in many calculations. They were discussed, in particular, in [45], where a linear basis in this set of constants up to weight 3 was explicitly described in terms of known transcendental numbers. Constants present in results for Feynman integrals up to weight 5 were also discussed in [46][47][48]. For example, one has where G R (a 1 , . . . , a n ) = Re G(a 1 , . . . , a n ; 1) , G I (a 1 , . . . , a n ) = Im G(a 1 , . . . , a n ; 1) , (3.8) refer to real and imaginary parts of the Goncharov polylogarithms. These constants satisfy various relations -see, e.g. [49]. We solved them in [51] up to weight six and presented a table of results for all these constants in terms of elements of some bases. Using these tables we obtained analytical results for the 51 MI presented in an ancillary file. Let us give two examples of these results. For the leading order term of the naive expansion at threshold of (2.8), we have πG I (0, 0, 0, r 2 ) − 751 320 π 4 G R (r 4 ) + 81 4 G I (0, r 2 ) 2 G R (r 4 ) + π 2 G R (r 2 , −1) The second example is the leading order term of the naive expansion at threshold of the The results at threshold typically do not have uniform weight. This is related to the fact that sometimes, we need to include power suppressed terms in order to identify the information about master integrals at threshold. In principle, one could search for a new integrals basis at threshold that has uniform weight.
The results we provide for the threshold integrals are up to certain orders in the ǫ expansion. Given that they originated from the information about integrals up to weight six away from threshold, one would expect these expansions to be sufficient for three-loop computations. At least, this was indeed the case in ref. [52]. In case further terms in the expansions are required, they can be obtained with the methods of this paper.

Conclusion
We have presented two more applications of the strategy to solve DE for Feynman integrals initiated in [4]. Historically, this project started from the evaluation of the threshold integrals which are single scale integrals and for which one cannot immediately apply the method of DE. However, as we explained in [16], one can introduce an extra scale, solve DE for the corresponding integrals and them then to turn back to the single-scale integrals. Of course, the integrals at general q 2 are also interesting in themselves. An application of them is described in the accompanying paper [52].
For convenience, our main results are available in terms of ancillary files.

A Pedagogical example: one-loop case
Here we give a pedagogical example of all steps of the method, at one loop. This has the advantage that all steps can be followed in detail. We consider the family of one-loop integrals of the type shown in Fig. 3, namely (A.1) We recall that p 2 i = m 2 and q 2 = (p 1 − p 2 ) 2 . Integral reduction shows that there are two master integrals. We choose them to be uniform weight integrals, following the procedure described in [4,29]. Our choice of basis is with c = (m 2 ) ǫ /Γ(1 + ǫ). Of course, the tadpole integral can be trivially computed. With our choice of normalization, it is given by f 1 = 1. For future reference, we note that the scalar triangle integral is related to the chosen basis via F 2,1,0 = −ǫ F 1,1,1 .
The next step is to derive differential equations for the master integrals in the kinematic invariants m 2 and q 2 . This is done via the chain rule. E.g., we can write with α = (q 2 − 2m 2 )/q 2 /(q 2 − 4m 2 ) and β = 2m 2 /q 2 /(q 2 − 4m 2 ). We can write the differential equations in a compact way Changing variables according to −q 2 /m 2 = (1 − x) 2 /x, this can be written more simply as or, equivalently, The differential equation has singularities at x = 0, −1, ∞. (The latter singularity can be seen by changing variables according to x → 1/x.) More generally, one finds that all integrals of this type, up to two loops, have singularities only at x = 0, 1, −1, ∞, see e.g. [27]. We can use the soft limit q 2 = 0, i.e. x = 1, to obtain a simple boundary condition, namely Equations like eq. (A.7) are easily solved in a series expansion in ǫ. The class of functions required are iterated integrals, with certain integration kernels d log α. One calls the set of allowed α letters (forming an alphabet specifying the class of functions). This above singularities correspond to the letters {x, 1 + x, 1 − x} (with only x, 1 + x required at one loop).
Solving eq. (A.7) with the boundary condition (A.8), we have, up to ǫ 3 , where H refer to harmonic polylogarithms [28]. Note that, by construction, all terms in the ǫ expansion have uniform weight. The three-loop calculation is very similar, except that we find that for some integrals, a new letter 1−x+x 2 is required, corresponding to a term d log(1−x+x 2 ) = dx (−1+2x)/(1− x + x 2 ) in the differential equations. As discussed above, this is related to sixth roots of unity. While the other singularities have a clear physical interpretation, the appearance of this new letter is somewhat surprising.
Since we chose x = 1 as boundary point, the above formula is valid near that point and can be analytically continued to other regions. We can consider e.g. the non-physical region for real x. There are several physical regions. Below threshold, x is complex and lies on the unit circle. The threshold is at x = −1. Above threshold, x is real and negative, and has a small positive imaginary part (originating from the Feynman i0 prescription).
Next, we consider the threshold limit. We parametrize x = e i(π−z) . In this way, we can analytically continue from our boundary point to the threshold. Taking the limit z → 0 of eq. (A.10), we find We now want to use this information to identify contributions from different scaling regions to this limit. In order to do this, we analyze the z → 0 limit of the differential equation (for fixed ǫ).
The equation takes the form where h is the boundary information at threshold (to be determined). Here .13) and the matrices P k are determined iteratively from the following equations [36] (A 0 − r)P r − P r A 0 = − r−1 s=0 A r−s P s . (A.14) For example, we have 1 + P 1 z + P 2 z 2 = 1 0 iǫz 1+2ǫ 1 + ǫz 2

12
(A. 15) We can see that the only terms in eq. (A.12) for which the ǫ → 0 and the z → 0 limit do not commute are contained in the matrix exponential z A 0 . The z −2ǫ terms correspond to a potential region. We are interested in the hard region, i.e. the limit without the z −2ǫ terms. So, we only missing piece of information is the boundary vector h(ǫ). We determine the latter, perturbatively in ǫ, by matching the small ǫ expansion of eq. (A.12) to eq. (A.10). In this way, we obtain, for the first few orders in ǫ, 16) This information allows us to compute the threshold integral. Throwing away all terms z jǫ with j = 0, and taking into account the factor relating f 2 and F 1,1,1 , we readily reproduce eq. (3.5), expanded in ǫ.