Evaluating ‘elliptic’ master integrals at special kinematic values: using differential equations and their solutions via expansions near singular points

This is a sequel of our previous paper where we described an algorithm to find a solution of differential equations for master integrals in the form of an ϵ-expansion series with numerical coefficients. The algorithm is based on using generalized power series expansions near singular points of the differential system, solving difference equations for the corresponding coefficients in these expansions and using matching to connect series expansions at two neighboring points. Here we use our algorithm and the corresponding code for our example of four-loop generalized sunset diagrams with three massive and tw massless propagators, in order to obtain new analytical results. We analytically evaluate the master integrals at threshold, p2 = 9m2, in an expansion in ϵ up to ϵ1. With the help of our code, we obtain numerical results for the threshold master integrals in an ϵ-expansion with the accuracy of 6000 digits and then use the PSLQ algorithm to arrive at analytical values. Our basis of constants is build from bases of multiple polylogarithm values at sixth roots of unity.


Introduction
Analytical results for Feynman integrals can, typically, be expressed in terms of harmonic polylogarithms [1] or multiple polylogarithms [2] which are very well mathematically studied special functions.For harmonic polylogarithms, one can apply the package HPL [3] which encodes various analytical properties and provides the possibility of numerical evaluation with a desirable precision.For multiple polylogarithms, one can either use the computer code [4] based on the GiNaC library [5] to obtain high-precision numerical values or construct a code based on their algebraic properties, for example, to reveal their behavior near singular points.Anyway, with a result in terms of these functions at hand, one can evaluate a Feynman integral at regular points and obtain expansions at singular points.
The possibility to arrive at a result written in terms of these functions exists if, within the method of differential equations [6][7][8][9][10][11], one succeeds to turn to a so-called canonical basis [12] of the master integrals 1 .It is well known that the ǫ-form of DE for a given set of the master integrals is not always possible 2 The simplest counter example is given by the two-loop propagator sunset diagram with three identical masses.In many irreducible cases, the lowest nontrivial term of ǫ-expansion is expressed in terms of elliptic integrals, and in what follows we will also refer to these cases as 'elliptic' in no relation with the functional form of the coefficients of ǫ expansion.
In situations without canonical bases, one can hope that the number of elliptic master integrals is small and try to obtain, in these cases, two and three-fold parametric representations suitable for numerical evaluation, and in all other cases to proceed with canonical subbases -see examples of such an approach in Refs.[20][21][22][23].On the other hand, it is quite natural to try to introduce new functions which would enable us to present results, in elliptic cases, in an analytical form.Multiple suggestions to introduce elliptic generalizations of multiple polylogarithms can be found in Refs.[24][25][26][27][28][29][30][31][32].However, these new functions do not have the same status as harmonic and multiple polylogarithms, as far as a detailed description of their properties and the possibility to evaluate them numerically are concerned.Moreover, the examples of successful treatment of ǫ-expansion in elliptic cases are, at most at the two-loop level.Anyway, we are very far, even in lower loops orders, from obtaining a complete description of a class of functions which can appear in results for Feynman integrals.
We advocated [33] an alternative way to solve differential equations which can be used also in 'elliptic' situations and illustrated it through a four-loop example.We considered multiloop Feynman integrals depending on one variable, i.e. with two scales where the variable is introduced as the ratio of these scales.We described an algorithm to find a solution of a given differential system in the form of an ǫ-expansion series with numerical coefficients.It is based on using generalized power series expansions near singular points of the differential system, solving difference equations for the corresponding coefficients in these expansions and using matching to connect series expansions at two neighbouring points.We provided a computer code where this algorithm is implemented for a simple example of a family of four-loop Feynman integrals where the ǫ-form is impossible.Using this code it is possible to evaluate master integrals at a given point as well as expansions at singular points with a required precision in an ǫ-expansion with a required number of terms.
Our present paper is a sequel of Ref. [33].The goal here is to apply our algorithm which is numerical in its character and the corresponding code in our example, i.e. fourloop generalized sunset diagrams with three massive and two massless propagators, in order to obtain new analytical results.We analytically evaluate the master integrals at threshold, p 2 = 9m 2 , in an expansion in ǫ up to ǫ 1 .Remember that the ǫ-form of the corresponding equations is impossible so that we cannot use the well-known procedure of using a solution at a general point in terms of well-established special functions and then turn to results at this singular point.However, although solutions to the differential equations look too complicated, the values of the master integrals at some special points can be conventional polylogarithmic constants.We will see that this is indeed true for the ǫ-expansion of the integrals of our family at the threshold.We use our code in order to construct a linear operator (a matrix) which renders the boundary conditions in one, suitable chosen, singular point to the coefficients of asymptotic expansion at the other point, p 2 = 9m 2 in our case.From these coefficients we extract the values of the integrals at p 2 = 9m 2 .
After reminding the main points of our setup in Section 2, we explain in Section 3 how we obtain high-precision values, up to 6000 digits, for the threshold integrals and then succeed in finding a relevant basis of constants in order to use the PSLQ algorithm [34].It turns out that the relevant bases of constants can be constructed starting from the bases of multiple polylogarithm values at sixth roots of unity, i.e. of the form G(a 1 , . . ., a w ; 1) where the indices a i are equal to zero or a sixth root of unity, with a 1 = 1, which were constructed in Ref. [35] up to weight six.We discuss various perspectives in Conclusion.

Setup
Differential equations for master integrals have the form where x is a dimensionless ratio of two scales for a family of dimensionally regularized Feynman integrals depending on two scales, J is a column-vector of N functions, and M is an N × N matrix with elements which are rational functions of x and ǫ.
The general solution of this linear system has the form where J 0 is a column of constants, and U is an evolution operator represented in terms of a path-ordered exponential We want to expand this operator in the vicinity of each singular point.Without loss of generality, let us consider the expansion near x = 0.It has the form where S is a finite set of powers, K λ 0 is an integer number corresponding to the the maximal power of the logarithm.
We assume that all the singular points of the differential system are regular so that we can reduce the differential system to a local Fuchsian form in any singular point.Therefore, we can reduce it at x = 0 to normalized Fuchsian form [19] by means of rational transformations.Let us assume that the system is in a global normalized Fuchsian form, i.e., and for any k = 0, . . ., s the matrix A k is free of resonances, i.e. the difference of any two of its distinct eigenvalues is not integer.In particular, the 'elliptic' cases, as a rule, can easily be reduced to a global normalized Fuchsian form.
As it is shown in Ref. [33], S and K λ can be determined and the difference equations for the coefficients C (n + λ, k) in (2.4) can be solved algorithmically.In fact, the idea to use series expansions at singular points and difference equations for the corresponding coefficients is very well known in mathematics.In high-energy physics, this strategy when evaluating Feynman integrals can be found, for example, in Refs.[20,[38][39][40][41]. Let us emphasize that our algorithm, with the current assumptions, provides solutions with no more than a linear growth of computational complexity with respect to a required number of terms.This is very important for a subsequent matching procedure which enables one to connect series expansions at two neighbouring points and thereby to obtain the possibility to evaluate Feynman integrals at any given point -see details in Ref. [33].

Four-loop generalized sunset diagram at threshold
As in Ref. [33] let us consider the example of the following family of four-loop Feynman integrals: where p is the external momentum and m is the mass of three lines.They correspond to the generalized sunset graph shown in Fig. 1.We introduce x = p 2 /m 2 .
Figure 1.The generalized sunset graph with two massless and three massive lines with the same mass.
The singular points of the differential equations are x 0 = 0, x 1 = 1, x 2 = 9 and x 3 ≡ x −1 = ∞.Using our algorithm we presented in Ref. [33] the code DESS to evaluate master integrals at a given point as well as expansions at singular points with a required precision keeping a required number of terms in an ǫ-expansion.
Let us now formulate our current goal: to evaluate master integrals of family (3.1) considered at threshold, p 2 = 9m 2 , i.e. exactly at the singular point x 2 = 9.In fact, for such integrals defined with the same general formula (3.1) we have now three master integrals which can be chosen as This can be done with the code DESS, where boundary conditions at the point x 0 = 0 were implemented.This point is, however, not a neighbour of x 2 = 9 so that matching is used twice when transporting information between x 0 and x 2 .This results in the necessity to evaluate much more terms of the series expansions at the three points {x 0 , x 1 , x 2 } in order to achieve a required precision which, as we will see later, should be very high because of a big number of constants relevant to results at x 2 .
In fact, the choice of the point x 0 to impose boundary conditions encoded in DESS can be explained by the fact that, generally speaking, setting p 2 = 0 for propagator integrals is equivalent to p = 0 and resulting vacuum Feynman integrals turn out to have just less indices.However, vacuum integrals can involve 'more complicated' constants.To solve our current goal, we can make a better choice to impose boundary conditions at the point x 1 for two reasons: this is now a neighbour of x 2 = 9 and the corresponding constants are multiple zeta values, logarithm of two and polylogarithms of one half.Indeed, master integrals at x = x 1 appeared in the calculations presented in Refs.[48,49] where they were evaluated using a onefold Mellin-Barnes representation.In particular, for F 1,1,1,1,1,0,...,0 at p 2 = m 2 = 1, the result is where ζ(. ..) are multiple zeta values.Here the result is restricted to contributions of weight seven and some efforts are needed to go further.
It turns out that the best way to impose boundary conditions is to choose x 3 = ∞ because the corresponding expansion is nothing but the large-momentum expansion [50][51][52][53] where, for our integrals (3.1), any term is a product of one-loop tadpoles and massless propagator integrals and can be evaluated via gamma functions at general ǫ.This provides any required accuracy and any required number of terms in ǫ-expansions in the boundary conditions.In the updated version of our code DESS, we introduce the possibility to impose boundary conditions at an arbitrary singular point.We added one more argument ns to the function DESS[rdatas, x, f(x), oe, np, nt, ns] which means the number of a singular point and this number is 1 for x 0 , 2 for x 1 , and 4 for x 3 .There is no sense to choose x 2 since this point is most complicated from the calculational point of view.We attach also two more auxiliary files: BoundaryConditions1 and BoundaryConditionsInf where analytic results for the boundary integrals are encoded.As before, the code and the auxiliary data can be downloaded from https://bitbucket.org/feynmanintegrals/dess.With the current version of DESS, we have obtained numerical results for the threshold master integrals in an ǫ-expansion up to ǫ 2 with the accuracy of 6000 digits for the corresponding coefficients.This took less than four hours on a desktop.As we will see shortly, such a big accuracy is needed for an application of the PSLQ algorithm.
The crucial point is a choice of a relevant basis of constants.A first hint comes from the known results for the two-loop sunset diagram at threshold [54,55] where one can observe multiple polylogarithm values at sixth roots of unity and π √ 3 .Let us also take into account that, at least according to Refs.[56][57][58], it might be reasonable to include into the basis the constant √ 3 separately.Therefore, we tried to use the bases connected with multiple polylogarithm values at sixth roots of unity and constructed in Ref. [35] up to weight six4 and √ 3. We consider bases of constants by including multiple polylogarithm values at sixth roots of unity up to weight six, i.e. of the form G(a 1 , . . ., a w ; 1) where the indices a i are equal to zero or a sixth root of unity, i.e. taken from the seven-letters alphabet {0, r 1 , r 3 , −1, r 4 , r 2 , 1} with and a 1 = 1.
The multiple polylogarithms are defined as with a i , z ∈ C and G(z) = 1.In the special case where a i = 0 for all i, the corresponding integral is divergent and instead one defines If a w = 0 and ρ = 0, then G(ρa 1 , . . ., ρa w ; ρz) = G(a 1 , . . ., a w ; z) so that one can express such MPL in terms of G(. . .; 1).The length w of the index vector is called the weight.One can consider separately the real and imaginary parts of the MPL G(a 1 , . . ., a w ; 1) = G R (a 1 , . . ., a w ) + i G I (a 1 , . . ., a w ) (3.6) For example, the elements of weight one are chosen in Ref. [35] as Let us denote by B R (w) (B I (w)) the bases generated by G R (a 1 , . . ., a w ) (G I (a 1 , . . ., a w )).
They include not only elements of the form G R/I (a 1 , . . ., a w ) but also products of constants of lower weights.The definitions of the bases can also be found in auxiliary files supplied with Ref. [35].They can be downloaded from http://theory.sinp.msu.ru/~smirnov/mpl6.As we shall see in our case in practice, when using the PSLQ algorithm, it is sufficient to use the bases B(w) = {B R (w), √ 3B I (w)} of weights w = 1, 2, . ... The element √ 3 does not contribute to the weight and it is 'imaginary' in its character, so that elements from √ 3B I (w) are 'real'.To get rid of √ 3 in our results, we can turn to rescaled imaginary elements via The numbers of elements in these bases B(w) are 3, 8, 21, 55, 144 for weights w = 1, 2, 3, 4, 5, correspondingly.If a constant is expected to be uniformly transcendental one can use these bases.Otherwise, one uses The numbers of elements in these bases are 4, 12, 33, 88, 232 for weights w = 1, 2, 3, 4, 5, correspondingly.
In simple situations, the number of available digits per constant in a basis can be as small as 7.In more complicated situations, with cumbersome coefficients in results, it can be more than 15.In our case, the accuracy of 2000 digits was quite enough to obtain results with PSLQ in an ǫ-expansion up to the finite part in ǫ, or, in other words, up to weight 4, in a straightforward way.Still at weight 5, it looks like the coefficients in results are more cumbersome and it is better to simplify our approach.
which differs from the uniformly transcendental basis of weight 5 given by Eq. (3.8) by adding three elements that are proportional to the leading terms of J 1 , J 5 , J 4 in their ǫexpansions.
Running PSLQ on a high-precision numerical value of the ǫ-term of J 6 , with this basis we obtain an analytical result from which we derive the ǫ-term of J 1 and thereby arrive at the following result for the first master integral at threshold − 960 GI (r 2 ) GI (0, 1, 1, r 4 ) + 568 GI (0, r 2 ) GI (0, 1, r 4 ) − 72172 81 GI (r 2 ) 3 GI (0, r 2 ) + 320 27 GI (r 2 ) GI (0, r 2 ) − 1152 GI (r 2 ) GI (0, r 2 , −1) + 14816 9 GI (r 2 ) GI (0, 0, 0, r 2 ) + 288 GI (r 2 ) GI (0, 1, r 2 , −1) + 1600 3 GI (0, r 2 ) 2 + 1680 GI (0, r 2 ) + 1136G R (0, 0, 1, r 2 , r 4 ) + 288G R (r 4 )G R (0, 0, r 4 , 1) − 420G R (0, 0, r 4 , 1) We apply this procedure based on uniformly transcendental bases also to J 2 and J 3 .Here we use, in a similar way, the following two linear combinations Then for the second and the third master integrals at threshold, we obtain the following results up to ǫ 1 : of threshold expansion but also leading terms of the form (9 − x) n−6ǫ (n is integer) and observed that the same bases also work and lead to analytical results via the PSLQ algorithm.Besides, proceeding in a similar way we arrived at the (Taylor) expansions at the singular point x = 0, with coefficients in terms of elements of our bases, in agreement with results for vacuum integrals [61,62].Therefore, we have demonstrated that although a canonical form of differential equations in our example is impossible and we don't know analytical results for the integrals, we can obtain analytical results for these integrals at some special kinematic points where the integrals are expressed in terms of usual polylogarithmic constants.
We have obtained results up to ǫ 1 but we believe that multiple polylogarithms values at sixth roots of unity form bases also at higher weights.The only possible complication when going to higher orders is connected with the fact that the size of the bases rapidly grows with the transcendental weight.We would like to emphasize that the bottleneck here is not connected with obtaining high-precision results using DESS, but rather with subsequent applications of the PSLQ algorithm.In fact, after the current calculation was done, we realized that we might use smaller (by 20-25 percents) bases defined in Ref. [39] via values of harmonic polylogarithms at sixth roots of unity.At least, the results presented in this paper can be expressed also in terms of these constants.Of course, it should be simpler to try to extend these results to higher weights using these bases.