Solving differential equations for Feynman integrals by expansions near singular points

We describe a strategy to solve differential equations for Feynman integrals by powers series expansions near singular points and to obtain high precision results for the corresponding master integrals. We consider Feynman integrals with two scales, i.e. non-trivially depending on one variable. The corresponding algorithm is oriented at situations where canonical form of the differential equations is impossible. We provide a computer code constructed with the help of our algorithm for a simple example of four-loop generalized sunset integrals with three equal non-zero masses and two zero masses. Our code gives values of the master integrals at any given point on the real axis with a required accuracy and a given order of expansion in the regularization parameter ϵ.


Introduction
Evaluating Feynman integrals with differential equations (DE) initiated in [1,2] and formulated as a method to evaluate master integrals in [3][4][5][6] became one of the most powerful methods already much time ago. Still this method is under development. In [7] it was suggested to turn from the basis of primary master integrals (i.e. revealed when solving integration by parts relations [8]) to the so-called canonical basis for which the right-hand side of the of system of DE is just proportional to = (4 − D)/2 and the singularities of the matrix on the right-hand side of DE are Fuchsian, i.e. have only simple poles in all the singular points of DE. The first algorithm to arrive at the canonical form was constructed in the case of one variable in ref. [9] (Such form of DE was called -form there). Besides the private implementation of this algorithm by its author and several other private implementations, two public implementations, Fuchsia [10,11] and epsilon [12], of the algorithm of ref. [9] are now available. 1 Once DE for master integrals are converted into an -form, i.e. one finds an appropriate linear transformation to a canonical basis, solving DE becomes straightforward, order-byorder in , at least in the form of iterated integrals, with possible complications connected with letters expressed in terms of square roots. Typically, the corresponding results are expressed naturally in terms of harmonic polylogarithms [15] or multiple polylogarithms [16]. These functions are very well studied. For harmonic polylogarithms, one can apply the package HPL [17] which encodes various analytical properties and provides the possibility of numerical evaluation with a very high precision. For multiple polylogarithms, one can use the computer code [18] based on the GiNaC library [19] to obtain high-precision numerical values, up to several thousand digits and more.
It is well known that the -form of DE for a given set of the master integrals is not always achievable by rational transformations. For massive internal lines it is often required to consider also transformations involving square roots. However, even using JHEP03(2018)008 transformations from this extended class it is not always possible to obtain an -form. 2 The simplest example where an -form is impossible is given by the two-loop propagator sunset diagram with three identical masses. In this example, as well in other known examples without -form, DE can still be reduced to the form where the right-hand side of the differential system is a linear function of .
However, 'integrating out' the constant term in such a form of DE appears to be an essentially more complicated problem. This can be seen in the known examples where results are expressed in terms of elliptic functions. In practice, it can happen that such 'elliptic' master integrals appear only in a small number of sectors. (A sector is specified by a distribution of the set of indices (powers of propagators) into positive and non-positive values.) A first example of a calculation of a full set of the master integrals with 'elliptic sectors' can be found in ref. [21]. A recent example of a calculation relevant to a physical problem is ref. [22], where elliptic functions appear only in two sectors and final results are expressed either in terms of multiple polylogarithms or, for the elliptic sectors and for some sectors involving complicated letters with square roots, in terms of two and three-fold iterated integrals suitable for numerical evaluation. Moreover, in refs. [23,24] a strategy to obtain parametric representations for master integrals applicable also in situations without -form was described and illustrated through one-, two-and three-loop examples.
Other examples of calculations of individual Feynman integrals in situations where -form is impossible can be found in [25,26] (see also references therein), where results are expressed in terms of elliptic generalizations of polylogarithms [25] or iterated integrals of modular forms [26]. One more class of elliptic generalization of multiple polylogarithms was recently introduced in ref. [27]. In particular, it includes functions appearing in the -expansion of the imaginary part of the two-loop massive sunset diagram. However, these new functions do not have the same status as harmonic polylogarithms and multiple polylogarithms, at least in the practical sense, i.e. there are no codes to evaluate them at a given point with a desired precision. Anyway, it looks like we are very far, even in lower loops orders, from answering the following question: 'What is the class of functions which can appear in results for Feynman integrals in situations where -form is impossible' ?
On the other hand, thinking positively, we may say that knowing a differential system and the corresponding boundary conditions gives almost as much information about Feynman integrals as knowing their explicit expressions in terms of some class of functions. In fact, some properties of the integrals are even more accessible via DE. In particular, singularities of DE provide a way to examine the branching properties of integrals. Numerical values of the integrals can be obtained from a numerical solution 3 of the differential system. Many computer algebra systems contain tools to solve this task (e.g. NDSolve procedure in Mathematica system). However, there is one complication that does not allow to use these tools immediately. Namely, we would like to keep as a variable and evaluate solutions of DE as series expansions in .

JHEP03(2018)008
The goal of the present paper is to describe an algorithm which enables one to find a solution of a given differential system in the form of an -expansion series with numerical coefficients. We describe such an algorithm in the case of Feynman integrals depending on one variable, i.e. with two scales where the variable is introduced as the ratio of these scales. As a proof of concept, we provide a computer code where this algorithm is implemented for a simple example of a family of Feynman integrals where the -form is impossible. The general idea behind our approach is to use generalized power series expansions near the singular points of the differential system and solve difference equations for the corresponding coefficients in these expansions. This idea is very well known in mathematics: already in the classical textbook [34] the algorithm to obtain first terms of expansion near a singular point is described. 4 In high-energy physics, its application to Feynman integrals can be found, for example, in refs. [21,[30][31][32]35] to evaluate expansions of solutions of DE at a given singular point by difference equations.
Though the approach of the present paper follows approximately the same lines, there are some peculiarities connected with our goals. Namely, we are aimed at the automatic calculation of multiloop integrals at an arbitrary kinematic point up to an arbitrary order in with an arbitrary precision. It is the latter goal that requires effective automatic generation of many terms of power series expansion 'on the fly'. In the next section, we present an algorithm to solve difference equations for coefficients of the series expansions at a given singular point. The advantage of the presented algorithm is that its computational complexity grows only linearly with the required number of terms. In section 3, we describe a matching procedure which enables one to connect series expansions at two neighboring points. In section 4, we describe a Mathematica code based on our algorithm and the matching procedure to evaluate master integrals in a simple four-loop example. This code as well as the corresponding auxiliary files can be downloaded from https://bitbucket. org/feynmanintegrals/dess. Then we conclude with a discussion of perspectives.

Generalized series expansion near a singular point
Let us have a differential system where J is a column-vector of N functions, and M is an N × N matrix with entries being rational functions 5 of x and . Below we will suppress in the arguments for brevity. We assume that all the singular points of the differential system are regular. This implies that we can reduce the differential system to a local Fuchsian form in any singular point. The general solution of this linear system has the form

JHEP03(2018)008
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 is well known that the expansion has the form where S is a finite set of powers, K λ 0 is an integer number corresponding to the maximal power of the logarithm. We have introduced the factor 1/k! for convenience. Our goal is to determine S, K λ , and the matrix coefficients C (n + λ, k). As to the latter, we are going to determine them via recurrence relations equipped with initial conditions. Since we assume that the differential system has only regular singular points, we can reduce it at x = 0 to normalized Fuchsian form [20] by means of rational transformations. For the sake of presentation, we will assume that the system is in 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. Note that the -form is only one example of normalized Fuchsian form, so we allow for a much wider class of differential systems which seems to be sufficient for any applications in multiloop calculations. In particular, the 'elliptic' cases, as a rule, can easily be reduced to a global normalized Fuchsian form. Besides, it is easy to generalize our algorithm properly if needed.
In order to obtain a recurrence relation of a finite order, we will first multiply both sides of eq. (2.1) by the common denominator xQ(x), where (2.6) By construction we have q 0 = 0. We will also define the polynomial matrix B (x, α) and its coefficients B m (α) by Then the recurrence relations read

JHEP03(2018)008
Here    denotes a (K + 1)N × N matrix built from blocks C (α, k) and the three-letter notation BJF stands for 'Block Jordan Form' defined as . Now note that the operator U , eq. (2.3), is determined up to a multiplication by a constant matrix from the right. We fix it by the condition This condition is, strictly speaking, mathematically incorrect when the distance between some eigenvalues of A 0 is larger than one, but it should be understood as the constraint on the leading terms of the expansion for each distinct eigenvalue. This condition gives us a way to determine S, i.e. the set of distinct eigenvalues of A 0 , and K λ , i.e. the highest power of the logarithm in front of x λ in x A 0 for each λ ∈ S, and the leading coefficients C(λ, k). We simply determine these parameters by representing Now note that the matrix −BJF(B 0 (λ + n), −q 0 , K λ ) on the left-hand side of eq. (2.8) is invertible for λ ∈ S and n > 0. Indeed, (2.11) and both q 0 = 0 and det(A 0 − λ − n) = 0, the latter is due to the absence of resonances in A 0 (since if det = 0 both λ and λ + n would be the eigenvalues of A 0 ). Therefore, we can rewrite recurrence relations (2.8) as where and use (2.12) together with the initial conditions determined 6 by eq. (2.10) in order to construct the generalized power series expansion (2.4). Note that the finite-order recurrence relation results in a linear growth of the computational complexity with the number of expansion terms.

For each λ ∈ S:
(a) the maximal power of the logarithm K λ and the leading coefficients C (λ, 0 . . . K λ ) defined by (2.10). To use the recurrence formula one has to take into account that C (λ + n, k) = 0 for n < 0.
Expansion in . So far we paid no attention to dependence of both the eigenvalues λ i and the coefficients C (λ + n, k). Fortunately, the question of how to expand obtained formulas in is very simple. Taking into account that the coefficients C (λ + n, k) are always rational functions of and the convergence of power series is determined by the factor x n , it is easy to conclude that the -expansion commutes with the summation. Therefore one can simply expand the coefficients C (λ + n, k) under the sum. As we should treat as the 'smallest' parameter, we have | log(x)| 1 and therefore, we also expand the factors x λ i .

Matching
The above considerations enable one to evaluate the evolution operator (2.3) within the convergence region of the power series (2.4). In order to perform an analytical continuation to the whole complex plane, one may use the same approach for the expansion around other singular points. Suppose that the next singular point closest to the origin is x = 1. We can construct the evolution operator (2.3) also in an expansion near this point, In general, due to the above mentioned freedom in the definition of the evolution operator we have where L is some constant matrix. If the convergence regions of the power series in U and U overlap, we may fix L by picking some point in the intersection of these regions. E.g. at x = 1/2 we have 7 L =Ũ −1 (1/2) U (1/2), i.e., finally, in the whole convergence region ofŨ we have Acting in the same way, we may, in principle, extend the definition of U onto the whole complex plane of x. In fact, this is a general approach to the analytical continuation Figure 1.
The generalized sunset graph with two massless and three massive lines with the same mass. of a function defined by a converging power series. In order to reach an arbitrary finite point of the complex plane, we are likely to need also expansions near the regular points (reducible to the considered case by putting A 0 = 0) and/or Möbius transformations of the variable. In the case where the singularities lie on the real axis and if we are interested in the evaluation of Feynman integrals for real x, we can avoid expansions near regular points and rely only on the Möbius transformations. Suppose, e.g., that we have the following sequence of the singular points Then for each 0 k s we make the variable change which maps the points x k−1 , x k , x k+1 to ∓1, 0, ±1, respectively. 8 It is convenient to choose the sign in such a way that the cuts of the non-integer powers and logarithms appearing in the series expansions coincide with the cuts of the integral.

Implementation
The four master integrals we evaluate form a basis of the following family of integrals: where p is the external momentum and m is the mass of three lines. They correspond to the generalized sunset graph shown in figure 1. We introduce x = p 2 /m 2 .

JHEP03(2018)008
We derive DE for J 1 in a straightforward way. When taking derivatives with respect to x one can apply LiteRed [39,40] to do this automatically. The derivatives are then expressed in terms of integrals of the given family. Solving integration by parts relations with an IBP-reduction code, 9 one expresses these derivatives as linear combinations of the primary master integrals and obtains a system of linear DE which has the form (2.1).
The matrix in the corresponding DE, as well as other entries connected with the described Mathematica code, can be can be downloaded from https://bitbucket.org/ feynmanintegrals/dess. One uses {M, T, Ti, Mf} = << "Data/TransformationData"; Here M = M is the matrix on the right-hand side of the original differential system, Mf = M f is the matrix of the transformed system, and T = T and Ti = T −1 are the transformation matrix and its inverse. 10 We turn to the basis J = T −1 · J 1 for which we have the matrix M f = T −1 (M · T − ∂ x T ) with normalized Fuchsian singularities at any singular point in the corresponding DE (2.1). The transformation T to the new basis was found with the help of the algorithm of ref. [9].
To fix boundary conditions we choose the point x = 0 where the integrals of the given family become vacuum integrals. To evaluate the four master integrals at x = 0 we derive onefold Mellin-Barnes representations for them and obtain the possibility to achieve a high precision for any given coefficient in the -expansion. We restricted ourselves to the accuracy of 500 digits but one can increase it to 1000 digits and more.
The singular points are x 0 = 0, x 1 = 1, x 2 = 9 and x 3 = x −1 = ∞. We solve difference equations for coefficients in series expansions near singular points according to the algorithm described in section 2. The corresponding results are encoded in a file present in the package:

{L, cis, cisrule} = Get["Data/BoundaryConditions"];
Here L = L is a constant matrix (see section 2) and the list cis = c defines the required information about the primary master integrals. These two quantities determine the column of constants J 0 in eq. The matching procedure described in the previous section is performed in our example as follows. The variable changes corresponding to the singular points are f 0 = x/(2 − x), f 1 = (x − 1)/(1 + 7x/9), f 2 = (9 − x)/(7 + x), f 3 = −9/(2x − 9). For example, the first function maps 0 to 0, 1 to 1 and infinity to −1. In new coordinates the radius of convergence is equal to 1, however, the convergence is very slow when approaching the boarder of the convergence domain.

JHEP03(2018)008
For adjacent regions i and i + 1 we search the best possible matching point which is such x that it lies between x i and x i+1 and that |f i (x)| = |f i+1 (x)|. In our case we result in matching points {−3, 3(3 − 2 √ 2), 3, 3(3 + 2 √ 2)}. The matching points are separating the singular points. We have To test our code we ran our procedure with oe = 15 and np = 75 at the sample points −10, −3/2, 1/3, 2/3, 2, 4, 12, 25 which lie between the singular and matching points and confirmed our results with the code FIESTA [41]. For example, at x 0 = 25, we obtain the following result (shown with a truncation to 10 digits) for the first primary integral:

JHEP03(2018)008
The -expansions here are implied with pulling out the standard factor e −γ E per loop where γ E is the Euler's constant.
In fact, the maximal order of expansion in and the maximal accuracy is determined by the boundary conditions where expansion of boundary vacuum integrals is included up to 3 with the accuracy of 500 digits. This results in an -expansion up to 3 of our primary master integrals. However, we recommend to set oe = 15 because high negative powers of appear in calculations. Moreover, we recommend to add the value 25 to the desired precision np, for a similar reason.
One more command of our code is denoted in the same way but has one more argument: The performance of the code depends essentially on the desired accuracy. Of course, the time needed for the evaluation with 500 digits is much more than with 50 digits because the matching procedure in the case of a high accuracy requires much more terms of expansions at singular points. For example, the evaluation at x = 9 with the accuracy of 500 digits takes around 40 minutes on a laptop. However, we did not try to optimize the procedure very much. In particular, one can think of parallelization and/or putting it on c++.
The results for the evolution operator given by DESS[rdatas, x, f(x), oe, np, nt] are linear combinations of powers of (±(x−x 0 )) n+j so that it is possible to select contributions for specific j at this level. For example, one can arrive at results for the naive part of the expansion of the primary master integrals near a given finite singular point by selecting only integer powers. In fact, near x 0 = 0 and x 0 = 1, we have only Taylor expansions of the master integrals in our example. We have exponents x j , with j=1,2,3,4, in the expansion at infinity but this does not mean that there is no naive expansion. The point is that the limit x → ∞ corresponds to the limit, where m 2 |p 2 |, so that the naive expansion in this limit reduces to the expansion of integrands in Taylor series in m 2 . If one is oriented at this very limit it is reasonable to introduce a dimensionless variable in another way, as x = m 2 /p 2 , and then the naive expansion will be in integer powers of this variable.

Conclusion
We have presented an algorithm for the numerical evaluation of a set of master integrals depending nontrivially on one variable at a given real point with a required accuracy. The algorithm is oriented at situations where canonical form of the DE is impossible. We have provided a computer implementation of the algorithm in a simple example. This code is similar in spirit to the well-known existing codes to evaluate harmonic polylogarithms [15] and multiple polylogarithms [16], where the problem of evaluation reduces to summing up appropriate series.
We hope that one can use our algorithm and implement it to evaluate master integrals in situations where an analytic evaluation is problematic. In fact, we have provided more than the code for the evaluation of the four master integrals we considered because our package includes tools for a decomposition of the real axis into domains, a subsequent mapping and an introduction of appropriate new variables. We are thinking of a more general package which would include an automation of as many steps of the presented algorithm as possible. Input data of this package would be a matrix in DE in the normalized Fuchsian form (defined near eq. (2.5)). Output data would be the evolution operator in an epsilon expansion up to a required order with a required accuracy. In addition to the existing tools, the future package needs at least an implementation of the algorithm of section 2 to solve difference equations for series expansions at the singular points.
Of course, one can hardy construct a general algorithm to fix boundary conditions because, usually, the choice of the corresponding point and the way to obtain data for the boundary conditions is done in every situation in a special way. Still we can suggest a format for including information about the boundary conditions for using it in our future package. Anyway, our future package would check if a given system of DE is already in a global normalized Fuchsian form, with singularities on the real axis, and, if this is true, the package would automatically construct the evolution operator in an expansion up to a required order.
We discussed the problem of evaluation of Feynman integrals with two scales, i.e. dependent on one variable, x. However, one can apply DE even in the case of one-scale integrals by introducing an extra scale, solving DE with the respect to the ratio of the two scales, x, and then picking a contribution to the expansion at the point where x tends to its original value [42]. The second form of the call of DESS allows one to find the coefficients of the expansion of the primary master integrals near a given singular point x 0 . Then it is easy to separate the 'naive' part of the expansion, i.e. the contribution of the non-negative integer powers of x − x 0 and to find the 'naive' values of the primary integrals at x = x 0 . For example, for the integrals considered in the previous section, this procedure can provide naive values at x = 1, i.e. integrals considered from the scratch with p 2 set to m 2 which are nothing but typical integrals appearing in the evaluation of the g − 2 factor.

JHEP03(2018)008
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.