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 JHEP04(2018)061 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][21].
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 JHEP04(2018)061 on several examples that the formulated approach can be successfully used for finding the general homogeneous solutions. In the parallel paper [22], 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 the following conditions are satisfied where N > 0 is some integer number. Basically, this condition corresponds to the moderate growth of the general solution of eq. (2.1), restricted by some power of ν when ν → ±∞, see ref. [23]. 1 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 1 As our main concern will be the application to the multiloop integrals, it is instructive to check whether this condition holds for the dimensional recurrence relations. We note that the multiloop integrals grow factorially when ν → −∞ due to the presence of the usual Γ(I −Lν) prefactor in the Feynman representation (I is the number of internal lines, L is the number of loops, ν = d/2), therefore, the condition (2.2) is violated. However, it is sufficient to pass from the master integral J(ν) to a new function f (ν) via J(ν) = Γ(I − Lν)f (ν) to eliminate this factorial growth.
In the examples presented in the next section we will apply a more general function change J(ν) = L i=1 Γ(ai − ν) f (ν), which also eliminates the factorial growth at µ → −∞, and use the freedom of choice of ai to reduce the degree N of the coefficients.

JHEP04(2018)061
where the coefficients c(σ) satisfy the recurrence n k=0 p k (σ + k)c(σ + k) = 0 (σ ∈ S + Z) , (2.6) λ 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. The convergence of the sum in eq. (2.5) depends on the asymptotics of the coefficients c(σ). This asymptotics can be determined from the recurrence relation (2.5) by the method described ref. [23]. In what follows, we assume that the convergence has been established and holds for any ν ∈ C\(S + Z). We admit though that the convergence issues here and below should be examined on a case-by-case basis, as we do in the subsequent section. Let us substitute the form (2.5) in the left-hand side of (2.1). We have It is easy to check that eq. (2.7) defines an entire function of ν, thanks to eq. (2.6). Indeed, the poles are only possible at ν ∈ S + Z, but taking the residue and using the recurrence relation (2.6) we see that those poles cancel. Now we will determine the conditions at which the entire function (2.7) vanishes. Let Then we have is the polynomial of the degree N − 1, the last equality in (2.10) defines the polynomials Q l (σ). The second term in the right-hand side of eq. (2.9) can be easily shown to vanish by shifting σ → σ + k:

JHEP04(2018)061
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 From now on we will put λ to be one of the roots of eq. (2.15). Therefore, we have N −1 equations (2.12)-(2.13). Suppose that we fix somehow the set S. Let us count the number of free parameters which we might tune to secure constraints (2.12)-(2.13). Since c(σ) satisfy the n-th degree recurrence relation (2.6), 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.12)-(2.13). 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.15) with maximal (minimal) absolute value, respectively. Therefore, whatever λ we take, the sums in (2.12)-(2.13) 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.6). In other words, for this specific choice of S the summation limits in (2.12)-(2.13) 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.12)-(2.13) 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 figure 1. The first master integral is trivial, S 1 = −Γ 3 (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. (3.33) is not of the form of eq. (2.1) because the degrees of the polynomial coefficients do not satisfy the conditions deg However, we can easily fix this by passing to a new function s 2 , e.g. with where

JHEP04(2018)061
Substituting this in eq. (3.3), we have 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.10). 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 JHEP04(2018)061 care of. Namely, we tacitly assumed that the sum σ∈S+Z p k (σ)λ k−σ c(σ) ν+k−σ in the first expression of eq. (2.11) 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. [23]. 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.11). 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.

JHEP04(2018)061
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 (ν) 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     and, therefore, where ω(ν) is some periodic function. Again, using educated guess and PSLQ we find
The first three properties secure that  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

Four-loop watermelon integral
Let us consider the four-loop watermelon tadpole topology. There are three master integrals depicted in figure 2. The first master integral is trivial, T 1 = Γ 4 (1 − ν), and the last two satisfy a coupled system of dimensional recurrence relations: These equations lead to the second-order recurrence relation for T 4 : We pass to a new function t 4 with Substituting this in eq. (3.33), we have

JHEP04(2018)061
The characteristic equation 3600 + 1098λ − 720λ 2 = 0 has the solutions The zeros of p 2 (ν) are So, we search for the solution in the form where the sequences {a 0 , a 1 , . . .} and {b 0 , b 1 , . . .} are determined recursively: In order to determine the only non-fixed coefficient b 0 we again apply the convergence constraint (similar to the previous example, the polynomial Q(ν, σ) vanishes) lim n→∞ n(a n + b n ) = 0 . 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 3 , 2 3 . However, again, it is simpler to use the symmetry of the homogeneous part of eq. (3.36) under the replacement followed by ν → 2 − ν. Therefore, we may write the second solution as (3.46) It is easy to check that t 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. (3.21) and (3.22) with p k and r now defined by eq. (3.37). The periodic factors can be fixed in a similar way as in the previous example. We finally find

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 figure 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 Passing to a new function t 10 , p k (ν + k)t 10 (ν + k) = 0 , The zeros of p 2 (ν) are So, according to the prescriptions of the previous section, we search for the solution in the form a −1 = 0, a n≥1 = 455n 3 −1050n 2 +691n−120 9n(3n−2)(13n−17) a n−1 + (n−1)(3n−4)(13n−4) 3n(3n−2)(13n−17) a n−2 , (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.10). We have Note that Q 2 = 0 as it should be. Thus, we have two conditions, ∞ n=0 a n λ −ν 1 −n Let us remark here that for practical reasons it might be simpler to use another approach to discover numbers in eq. (3.62). Namely, we calculate with high precision the right-hand side of the recurrence relation (3.53) substituting t 10 (ν) with the contribution of each term in square brackets in eq. (3.57) for two random values of ν. Then we have a linear system of the form 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 It is easy to check numerically that t 10h (ν). Examining the analytical properties of the integrals T 9,10 in a similar way as for the previous example, we obtain the final result

Three-loop box in special kinematics
All previous examples related the integrals with massive internal lines. In order to demonstrate the applicability of our method also to the integrals with massless lines, we present here the result for the integral corresponding to the diagram in figure 4. This integral is relevant for the three-loop hard contributions to the energy levels of the weakly bound QED systems, like positronium. There are two master integrals in the highest sector, so the integral without squared propagators satisfies a second-order recurrence relation. Passing to j 5 (ν) connected with the original integral via we have the following recurrence relations 3(ν − 1)(4ν − 5)j 5 (ν) − 28ν 2 − 28ν + 9 j 5 (ν + 1) + 4ν(4ν + 1)j 5 (ν + 2) = . . . , (3.69) where dots in the right-hand side denote simpler master integrals. Using the described approach, we find the following homogeneous solutions for j 5 : where The inhomogeneous solution can be found in a straightforward way. In particular, in a parallel paper [22], we consider the unitarity cut of the diagram 4, where all massless lines are cut. It appears that for the cut diagram the specific homogeneous solution is zero in this case and the whole result comes from the inhomogeneous terms.

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,25] 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 figure 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 [22] 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:

JHEP04(2018)061
The high-precision numerical computation of S 2 near d = 2 can be obtained equally easily, Unfortunately, the basis of transcendental numbers involved in the above expansion appears to be not quite clear to us 3 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.5), (2.6). 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 f (ν) →f (ν) = n−1 k=0 c k (ν)f (ν + k) is very desirable. These transformations may lead to essential simplification of the coefficients of the recurrence relations (2.1).
Finally, we note that the current method is also applicable to the integrals depending on several parameters. There are some complications though. Namely, when solving the characteristic equation (2.15), one obtains λ 1 , . . . λ n as functions of the parameters. The whole space of parameters is, therefore, split into regions where in the neighboring regions the definitions of λ max and/or λ min changes. Also, the rate of convergence of the sums depends on the values of parameters. Therefore, we would suggest to try the differential equations approach to such integrals, while using the present method for fixing the boundary conditions.