Meromorphic solutions of recurrence relations and DRA method for multicomponent master integrals

We formulate a method to find the meromorphic solutions of higher-order recurrence relations in the form of the sum over poles with coefficients defined recursively. Several explicit examples of the application of this technique are given. The main advantage of the described approach is that the analytical properties of the solutions are very clear (the position of poles is explicit, the behavior at infinity can be easily determined). These are exactly the properties that are required for the application of the multiloop calculation method based on dimensional recurrence relations and analyticity (the DRA method).


Introduction
The ability to perform multiloop calculations is very essential for obtaining high-precision theoretical predictions in particle physics. Such predictions are, in particular, indispensable for the ongoing searches for New Physics. Since 1980s the field of multiloop calculations has experienced an explosive development in terms of the available methods and tools. Two major technical achievements in this region are the integration-by-parts reduction [1,2] and differential equations method [3,4]. Thanks to these two techniques, the present frontiers of the multiloop calculations reside somewhere close to NNLO calculation of the differential cross sections with up to 4 or 5 particles involved (i.e., two loops and 4 parameters). However, some quantities, which depend only on one invariant, deserve and allow for the N 3 LO and N 4 LO accuracy. The differential equations can not help, at least directly, in this case. Probably, the most celebrated example are perturbative contributions to the anomalous magnetic moment of the electron and muon. Due to the efforts of Kinoshita's group (see Ref. [5] and references therein), the QED contributions to electron g − 2 were known numerically up to the four-loop accuracy, and only recently these results have been independently verified by Laporta [6]. In fact the approach of Ref. [6] would lead to completely analytical form of the four-loop QED contribution to electron g −2 if we understood better which transcendental numbers might enter the final expression. This is, probably, another great achievement of contemporary multiloop methods -learning how to guess the analytical results from the high-precision numerical ones.
One of the methods which may be used for the calculation of the one-scale integrals is the method of Ref. [7] based on the difference equations with respect to the power of one of the massive propagators. Another natural idea is to use the recurrence relations with respect to the space-time dimension d, [8]. In Ref. [9] the DRA method was formulated which uses the dimensional recurrence relations and the analyticity of the integrals as functions of d. The key idea of the DRA approach is to use the analytical properties of the integrals in order to fix the form of the homogeneous solutions up to several constants which should be determined by other methods (e.g. from Mellin-Barnes representation). The results of the DRA method have the form of multiple convergent sums with factorized summand. The latter property allows for their fast high-precision evaluation. The corresponding algorithms were implemented in the SummerTime program [10], which allows one to obtain the ǫ expansion of the integrals near any dimensionality d with high-precision numerical coefficients. This method was successfully applied in many physical calculations [11][12][13][14][15][16][17][18][19][20].
In order to apply the DRA method it is necessary to construct the homogeneous solutions with known analytical properties. While for the first-order recurrence relations this is a trivial task, for the higher-order recurrence relations this is not so. In Ref. [17] the homogeneous solutions of the second-order differential equations were obtained from the explicit evaluation of maximally cut integrals. This approach related a concrete example of the integrals and can be hardly generalized. Therefore, the application of the DRA method to the topologies containing the sectors with several master integrals (multicomponent master integrals in terminology of Ref. [17]) is complicated by the necessity to construct the homogeneous solution of the higher-order recurrence relations with known analytical properties.
To explain the main goal of this paper let us consider the following example. Suppose that we want to find the meromorphic function f (ν) (ν ∈ C) which obeys the first-order recurrence relation where P and Q are some polynomials. This is a simple problem: we can take, e.g., the function where p, q, a k , and b l are determined by the factorization of the polynomials P and Q: Suppose now that we have the second-order recurrence relation, e.g.
Then the task of finding the meromorphic solution is not simple anymore. In principle, one can try some integral transformations, e.g. the Mellin transformation, in order to turn Eq. (1.4) into the differential equation. However, the order of this differential equation and the number of its singular points grows rapidly with the degree of the polynomial coefficients P , Q, and R which makes this approach impractical in multiloop calculations. The main goal of this paper is to suggest the alternative approach based on searching the solutions in the form of the sum over poles with coefficient defined recursively. We will show on several examples that the formulated approach can be successfully used for finding the general homogeneous solutions. In the parallel paper [21], we present the package DREAM which allows one to automatically construct the inhomogeneous solutions. In that paper, we use this package and the homogeneous solutions found in the present paper to calculate the corresponding multiloop integrals.

Meromorphic solutions of recurrence relations
Let us consider the n-th order homogeneous recurrence relation where p k (ν) are some given polynomials and we shifted the arguments of p k for convenience of further considerations. We will assume that deg p 0 = deg p n = N , deg p i N . We are interested in n independent solutions f 1 (ν), . . . , f n (ν), each being a meromorphic function of the complex variable ν. Then the general solution can be represented as where ω 1 (ν), . . . , ω n (ν) are arbitrary periodic functions with period 1. Note that for n = 1 we can explicitly write such a solution as where p kN is the leading coefficient of p k (ν) = p kN ν N + . . ., and ν ki are the zeros of p k (ν). When n > 1 the problem of finding the meromorphic solutions is much more difficult. Let us search for the solution in the form where the coefficients c(σ) satisfy the recurrence n k=0 p k (σ + k)c(σ + k) = 0 (σ ∈ S + Z) , (2.5) λ is some number, and S = {σ 1 , . . . , σ m } is some finite set with no resonances, i.e. the difference between any two distinct elements of S is non-integer, σ i − σ j ∈ Z. For the moment we assume that the sum over σ converges absolutely for any ν ∈ C\(S + Z). Let us substitute the form (2.4) in the left-hand side of (2.1). We have It is easy to check that Eq. (2.6) defines an entire function of ν, thanks to Eq. (2.5). Indeed, the poles are only possible at ν ∈ S + Z, but taking the residue and using the recurrence relation (2.5) we see that those poles cancel. Now we will determine the conditions at which the entire function (2.6) vanishes. Let us write p k (ν+k) Then we have is the polynomial of the degree N − 1, the last equality in (2.9) defines the polynomials Q l (σ). The second term in the right-hand side of Eq. (2.8) can be easily shown to vanish by shifting σ → σ + k: Therefore, we have to require that the following N conditions hold: The last condition is somewhat special, since Q N −1 does not depend on σ. Let p k (ν) = p kN ν N + . . . Then it is easy to see that Q N −1 = n k=0 p kN λ k . If λ is a root of the characteristic equation n k=0 p kN λ k = 0, (2.14) we have Q N −1 = 0. From now on we will assume assume that λ is one of the roots of Eq. (2.14). Therefore, we have N − 1 equations (2.11)-(2.12). Suppose that we fix somehow the set S. Let us count the number of free parameters which we might tune to secure constraints (2.11)-(2.12). Since c(σ) satisfy the n-th degree recurrence relation (2.5), we might fix arbitrarily the coefficients c(σ l ), c(σ l + 1), . . . , c(σ l + n − 1) for each σ l ∈ S. Therefore, we have n × m parameters (m = |S|). A naive counting would be that when n × m N , we have a nontrivial solution of (2.11)-(2.12). However, in general the convergence requirements should also be treated properly. Suppose that we fix c(σ l ), c(σ l + 1), . . . , c(σ l + n − 1) arbitrarily. Then the asymptotics of the sequence c(σ l + k) is likely to behave as λ k max when k → +∞ and λ k min when k → −∞. Here λ max (λ min ) denote the roots of the characteristic equation (2.14) with maximal (minimal) absolute value, respectively. Therefore, whatever λ we take, the sums in (2.11)-(2.12) are likely to diverge at the upper and/or at the lower limit. Fortunately, there is a way out of these convergence problems. Let us fix S = {σ 1 , . . . σ N }, where σ 1 , . . . σ N are the zeros of the polynomial p n (σ). Then putting c(σ l + n) = 0 for all negative n and c(σ l ) = const = 0 is obviously consistent with the recurrence relations (2.5). In other words, for this specific choice of S the summation limits in (2.11)-(2.12) are restricted from below, σ ∈ S + Z + , where Z + = {0, 1, 2, . . .}. Therefore, we should only care for the convergence at the upper limit. The choice λ = λ max eliminates at least exponential divergence and we shall see on the specific examples that the sums converge. Then, N − 1 equations (2.11)-(2.12) form a linear system for N variables c(σ 1 ), . . . , c(σ N ) from which these variables can be determined up to the arbitrary common factor.
Therefore, we have a receipt to find at least one meromorphic solution of the recurrence relation (2.1). In fact, setting S to be the set of zeros of p 0 (σ) and repeating the similar analysis, we may find the second solution for which the summation goes downwards. For the second-order recurrences these two solutions constitute a complete set, apart from possible degeneracy.
In the next Section we demonstrate that the presented approach can be successfully applied to very nontrivial multiloop integrals.

Three-loop massive sunrise integral on the pseudo-threshold
Let us consider the three-loop massive sunrise topology on the threshold. There are three master integrals depicted in Fig. 1. The first master integral is trivial, Figure 1. Master integrals. Solid lines correspond to the propagators 1/(k 2 − 1). and the last two satisfy a coupled system of dimensional recurrence relations: Here and below ν = d/2, where d is space-time dimension, and α n = α(α + 1) · · · (α + n − 1) is Pochhammer symbol. These equations lead to the second-order recurrence relation for S 2 : The homogeneous part of Eq. However, we can easily fix this by passing to a new function s 2 , e.g. with where We now apply the method described in the previous section to find the homogeneous solution of Eq. (3.6), i.e., solution of The characteristic equation 144λ 2 + 108λ − 36 = 0 has the solutions The zeros of p 2 (ν) are Here we note that λ 1 = −1 which determines the behavior of the homogeneous solution at ν → +∞ is negative. Therefore it is convenient to pass to the functions 2 (ν) = a(ν)s 2 (ν), where a(ν) = −a(ν + 1) is some anti-periodic function which we choose as a(ν) = sin(πν). Such a substitution flips the signs of p 1 (ν) and λ 1,2 . Then, according to the prescriptions of the previous section, we search for the solution in the form where we have taken into account that −λ 1 = 1 and the sequences {a 0 , a 1 , a 2 , . . .} and {a 0 , b 1 , b 2 , . . .} are defined recursively: In order to determine the unfixed coefficient b 0 in Eq. (3.40), we calculate the polynomial Q(ν, σ), Eq. (2.9). It appears that Q = 0 and one might think that Eq. (3.11) is a solution for arbitrary b 0 . There is, however, a convergence issue which we have to take care of. Namely, we tacitly assumed that the sum σ∈S+Z p k (σ)λ k−σ c(σ) ν+k−σ in the first expression of Eq. (2.10) converges. Specialization to the case under consideration is that should converge. The asymptotics of the sequences a n and b n can be found along the lines of Ref. [22]. We have a n , b n ∼ 1 n , (3.14) therefore, a necessary condition of the convergence of (3.13) reads lim n→∞ n(a n + b n ) = 0 . Note that the sum in Eq. (3.13) diverges logarithmically even after we put b 0 to the value (3.17), however this is sufficient for the consistency of summation variable shifts and changing the summation order used to prove Eq. (2.10). Therefore, explicitly, we have the following homogeneous solution where a n and b n are defined by Eqs. (3.12) and (3.17).
The second homogeneous solution can be found in a similar way by taking S to be the set of zeros of the trailing coefficient p 0 (x), i.e., S = 1 6 , 5 6 . However it is simpler to use the hidden symmetry of Eq. (3.8). Namely, Eq. (3.8) is invariant under the replacement followed by ν → 2 − ν. Therefore, we may write the second solution as It is easy to check numerically that s (2) 2h (ν) is independent of s 2h (ν), i.e., that their ratio is not a periodic function.

Fixing periodic factors
Let us explain now how the homogeneous solutions found can be used within the DRA method and allow one to fix the form of the result for s 2 (and for S 2 ). We write the general solution of (3.6) in the form 2h (ν + 1) s 2h (ν + 1) s 2h (ν) where ω 1,2 (ν) are arbitrary periodic functions and t 4ih (ν) is the special solution of the inhomogeneous equation. The latter can be written as where and p 0 , p 1 , p 2 , r are defined in Eq. (3.7). Note that s 2ih (ν) does not have poles in the region Re ν < 1. Now, according to the standard prescription of the DRA method, we rewrite Eq. (3.20) as where T (ν) = s (1) 2h (ν + 1) s (2) 2h (ν + 1) s 2h (ν)s (2) 2h (ν + 1) . (3.25) Let us remind that the Casoratian W (ν) satisfies the equation and, therefore, where ω(ν) is some periodic function. Again, using educated guess and PSLQ we find 2h (ν) · T (ν) = (0, 1), in particular, at ν = ± arccos 3 2 2π . The first three properties secure that = 0 for b 1 and b 2 . Using PSLQ for the last constraint, we obtain At this point we have only one constant to be fixed. We fix it by explicitly calculating S 2 at d = 1 to find S 2 (1/2) = π 3/2 3 . We obtain Thus, our result is where s 2ih , s 2h , and s (2) 2h are defined in Eqs. (3.21), (3.18), and (3.19), respectively. In the next Section we present some numerical results obtained with this representation.
It is easy to check that t (2) 4h (ν) is independent of t (1) 4h (ν). Moreover, using PSLQ, we obtain We write the general solution of (3.36) in the form where ω 1,2 (ν) are arbitrary periodic functions. The inhomogeneous solution t 4ih (ν) is defined in the same way as in Eqs.

Four-loop cat-eye graph
Let us now consider the four-loop cat-eye graph. In the highest sector there are two master integrals depicted in Fig. 3. The dimensional recurrence relations for these two master integrals have the form where the dots in the right-hand side denote the contribution of the simpler masters. The homogeneous part of the second-order recurrence relation for the dotted integral T 10 has the form − 9(13ν − 4)T 10h (ν) − 2(2ν − 1) 455ν 3 − 1050ν 2 + 691ν − 120 νT 10h (ν + 1) Passing to a new function t 10 , where The zeros of p 2 (ν) are So, according to the prescriptions of the previous section, we search for the solution in the form  (3.58) The above recurrence relations determine a n 0 , b n 0 , and c n 0 via the starting coefficients a 0 , b 0 , and c 0 , respectively. In order to determine (up to an overall constant) the starting coefficients, we calculate the polynomial Q l (σ), Eq. (2.9). We have   where c ij are some high-precision numbers. Then the ratios (3.62) can be obtained from this system. Therefore, we have the following homogeneous solution where λ 1 = 1 54 35 + 13 √ 13 , a 0 = 1, b 0 = 7/18, and a n>0 and b n>0 are defined recursively by Eq. (3.58).

Computational issues
First, we would like to note that the representations like (3.18) are ideally fitted to obtaining the ǫ expansion. One should simply expand under the summation sign which amounts to the replacement One might question the possibility to obtain high-precision results for the representations like (3.18). Indeed, the sum converges harmonically, as n 1/n 2 , and the direct summation of N terms would give only log 10 N decimal digits. Fortunately, the convergence acceleration technique described in Refs. [10,24] can be successfully applied. In the ancillary file Meromorphic.nb, we present the Mathematica procedures for the calculation of all homogeneous solutions obtained in this paper, namely, s (1,2) 2h , t (1,2) 4h , t (1,2) 10h , and j (1,2) 5h . In Fig. 5 we present some timing results for the calculation of the homogeneous solutions. As it can be seen, the time of calculation scales as O(p α ), where p is the requested precision and α ≈ 2.3. The almost quadratic dependence on p could have been anticipated if one takes into account that the computational complexity of the convergence acceleration is quadratic in the number of terms, and the latter is approximately proportional to p.
In a parallel paper [21] we present the DREAM program suitable for the required highprecision numerical computation using the results for the homogeneous solutions obtained in the present paper. Here we will only demonstrate the effectiveness of our method by presenting the expansions of S 2 near d = 3: Unfortunately, the basis of transcendental numbers involved in the above expansion appears to be not quite clear to us 2 which prohibited the use of PSLQ algorithm.

Conclusion
In this paper we have formulated an approach to the construction of the solutions of the higher-order recurrence relations in the form of a sum over poles, with coefficients defined recursively, (2.4), (2.5). The main advantage of the described approach is that the analytical properties of the solutions are very clear (the position of poles is explicit, the behavior at infinity can be easily determined). These are exactly the properties that are required for the application of the DRA method. It is quite remarkable that three out of four examined recurrence relations (and also some more not presented here) exhibited a hidden symmetry ν → a − ν with a equal to 2 or 3. The only exclusion is the homogeneous recurrence relation for the cat-eye topology, Eq. (3.53). However for this topology one could speculate that there might exist a linear combinationt 10 (ν) = c 0 (ν)t 10 (ν) + c 1 (ν)t 10 (ν + 1) with coefficients c 0,1 (ν) being the rational functions of ν, such that the dimensional recurrence fort 10 again has the symmetry ν → a − ν. In general, a better understanding of the transformations is very desirable. These transformations may lead to essential simplification of the coefficients of the recurrence relations (2.1).