Two-loop diagrams in nonrelativistic QCD with elliptics

We consider two-loop two-, three-, and four-point diagrams with elliptic subgraphs involving two different masses, m and M . Such diagrams generally arise in matching procedures within nonrelativistic QCD and QED and are relevant, e.g., for top-quark pair production at threshold and parapositronium decay. We present the obtained results in several different representations: series solution with binomial coefﬁcients, integral representation, and representation in terms of generalized hypergeometric functions. The results are valid up to terms of O (ε) in d = 4 − 2 ε space-time dimensions. In the limit of equal masses, m = M , the obtained results are written in terms of elliptic constants with explicit series representations.


Introduction
During the last two decades, great progress has been made in the calculation of multiloop Feynman diagrams. This was possible thanks to the appearance of numerous new techniques. In particular, many advancements became possible due to the development of the theory of multiple polylogarithms [1][2][3][4][5][6]. A summary of their algebraic and numerical algorithms may be found in Ref. [7]. At special values of their arguments, polylogarithms may transform to multiple Euler-Zagier sums [8][9][10][11], "sixth-root of unity" constants [12][13][14][15][16][17][18][19], multiple binomial sums [15,16,[20][21][22][23][24][25][26][27][28][29][30], and others. Calculations yielding polylogarithms are often related to massless problems, sometimes also to massive ones where the unitary cuts of a Feynman diagram do not cross more than two massive lines. It should be noted that, in dimensional regularization with d = 4 − 2ε space-time dimensions, all one-loop diagrams with arbitrary masses and kinematics can, at least in principle, be expressed in terms of polylogarithms at any order in ε. 2 Going beyond the class of polylogarithms is a challenging task, which requires additional investigations. If we are aiming at a representation of a result that is more complicated than just a number, then we should answer the question as to what kinds of functions are required for that. Is it possible to classify, enumerate, and build higher functions, so that also the results at higher orders in the ε expansion, and eventually at higher loop orders, can be expressed within that class of functions? Finally, we need a stable and fast numerical library for all such new functions introduced. The knowledge of the whole class of functions and its algebraic structure enables us to apply a very powerful method of calculation, based on differential equations [32][33][34][35][36]. Knowing the structure of the answer, we are able to construct an appropriate ansatz for the solutions and to solve the corresponding differential equations. It is probably fair to say that the above program has been elaborated, to a large extent, only for the class of polylogarithmic functions so far. Nevertheless, there has recently been a lot of progress in understanding the simplest functions beyond multiple polylogarithms, the so-called elliptic polylogarithms . However, it is already clear that, even at two loops, elliptic polylogarithms are not sufficient for all possible applications, as there can either be several elliptic curves [58,59] or completely new functions present [45,[60][61][62][63].
In this paper, we proceed with the study of Feynman diagrams possessing elliptic structure. We consider diagrams with only one or two independent parameters. The former case is represented by the so-called single-scale integrals, which, by dimensional reasons, can be written as a product of a scale and a numerical factor. In the latter case, the considered diagram is expressed in terms of a function of one variable. The analytic structure of such functions can be analyzed with the help of differential equations which these functions obey. The single-scale limit is obtained by setting an existing variable to the appropriate fixed value.
To be specific, we consider diagrams for 2 → 2 processes with two real photons (γ γ ) or gluons (gg) in the initial state and two on-shell massive fermions (ff ) in the final state. In addition, we stick to threshold kinematics. For example, we have γ (q 1 The physical motivation to study such processes with threshold kinematics is the following. Such kinematics appear, for example, in the matching procedures of QCD or QED to the corresponding nonrelativistic effective theories, such as NRQCD and NRQED, where the above processes define the corresponding hard Wilson coefficients. The NRQCD applications are related to heavy-quarkonium production and decays (see Refs. [64][65][66][67][68] and references cited therein) as well as to near-threshold tt production (see Refs. [58,59,[69][70][71][72][73][74][75][76][77][78][79][80][81][82] and references cited therein), which offers the opportunity to considerably improve the accuracy of t -quark mass measurements at the CERN Large Hadron Collider and future linear colliders. As for NRQED applications, the best known one is probably the calculation of the parapositronium decay rate [83][84][85][86][87]. As already mentioned above, there are no new functions beyond polylogarithms at the oneloop level. So, the elliptic structures first appear at next-to-next-to-leading order (NNLO). One of the most complicated diagrams containing elliptic structure is the nonplanar one shown in Fig. 1(a). Applying integration-by-parts (IBP) identities [88,89], this diagram can be reduced to a sum of so-called master integrals with rational coefficients. When both massless lines are contracted, we obtain the diagram shown in Fig. 1(b) upon the replacement M = m. In this paper, we evaluate the latter diagram and all its subtopologies obtained by contracting some of its lines. Moreover, as by a product, we even solve the more general problem with two different masses, as shown in Fig. 1(b).
Here and in the following, we work in dimensional regularization with d = 4 − 2ε space-time dimensions. In our previous paper [90], we presented the results for sunset diagrams obtained by contracting lines 1 and 2 in Fig. 1(b). The results of Ref. [90] were obtained up to and including the O(ε) terms and represented as fast converging series in the mass ratio x = m 2 /M 2 . The O(1) terms may also be found in Section 4.1. It is then easy to rewrite a given such series in terms of generalized hypergeometric functions p F p−1 [91], as d dδ where ζ n = ζ(n) denotes Riemann's zeta function. A similar form was found for the Baxter Q function in Refs. [92,93]. Following our study in Ref. [90], it was shown in Ref. [94] that a similar hypergeometric representation is applicable for the considered sunset diagrams even without ε expansion. In the present paper, we extend the results of Ref. [90] for the case of threeand four-point diagrams (see Fig. 1) with O(1) accuracy. This paper is organized as follows. In Section 2, we recall the central relation between the one-and two-loop diagrams, Eq. (3). Section 3 contains the details of the one-loop calculations. In Section 4, we derive the series, integral, and hypergeometric representations of the considered two-loop diagrams in the general case of non-equal masses, M = m. Subsection 4.4 contains our results for the equal-mass case. Finally, in Section 5, we present our conclusions. In Appendix A, we explain the Frobenius solution for a system of differential equations.

Relation between one-and two-loop diagrams
As explained in Section 1, all two-loop Feynman diagrams considered in this paper have self-energy insertions built from two massive propagators as indicated in Fig. 1(b). These can be reduced to one-loop diagrams by representing the loop with the two massive propagators as an integral whose integrand contains a new propagator with a mass that depends on the variable of integration [20,21,32,33,90,95]. Graphically, this procedure has the following form: where the loop involving the two propagators with mass square M 2 is replaced by one propagator with mass square M 2 1 /s + M 2 2 /(1 − s). The numbers a and b, indicating the powers of the scalar propagators, are called the propagator indices.
Equation (3) is easily derived from the Feynman parameter representation and was introduced in Euclidean and Minkowski spaces in Refs. [32,33] and [20,21], respectively. Here, we work in Minkowski space and thus follow Refs. [20,21].
As in Ref. [90], we adopt the following strategy. First, applying Eq. (3), we write expressions for the considered two-loop diagrams as integrals of one-loop diagrams involving propagators with masses that depend on the variable of integration. Next, in Section 3, we evaluate the oneloop integrals thus obtained and then, in Section 4, use these expressions to reconstruct the results for the two-loop diagrams involving two different masses with the help of Eq. (3). A similar strategy was also adopted for the calculation of certain four-loop tadpole diagrams in Ref. [95].

One-loop integrals
Let us consider the following one-loop box integral with two different masses, m and M, and propagator indices a 1 , a 2 , a 3 , a 4 : where the kinematics in Eq. (1) is implied. Using IBP, the integral in Eq. (4) can always be reduced to a set of scalar master integrals with propagator indices 1 or 0. Some of these master integrals are shown in Fig. 2. , respectively. Notice that, in the course of the reduction, more master integrals appear: two bubbles, I 1000 and I 0001 , two more vertices, I 1110 and I 1101 , and one more self-energy, I 0110 . However, all these additional integrals only contain simple polylogarithmic structures, and their two-loop counterparts have no elliptic content. We omit these diagrams in our discussion, except for those that enter as subgraphs in the integrals I mmmM 0111 and I mmmM 1111 .

Two-point case
In the case of the two-point diagrams, IBP provides us with the relations which can be considered as the differential equations for the initial diagrams I mM 0101 and I mM 0011 , as The integrals I M 000a and I m 00a0 in Eq. (5) are tadpoles with masses M and m, respectively, Similarly, using IBP and substituting m = 1 and d/dM 2 = −x 2 d/dx, we can deduce a differential equation for I (1) 12 = I mM 0102 , which is valid up to O(ε) corrections. Equation (8) is easy to solve and, taking the boundary conditions at x = 0 into account, we have or, in terms of a series solution, Now, IBP relations allow us to write down answers for the I (1) a 1 a 2 integrals with other values of propagator indices. For example, we have For m = M (x = 1), we have In the case of the integral I (2) 12 = I mM 0012 , the application of IBP relations leads to the following differential equation: which is again valid up to O(ε) terms. The solution of Eq. (13) is straightforward and, taking the boundary conditions at x = 0 into account, yields or, in terms of a series solution, The expressions for the I (2) a 1 a 2 integrals with other values of propagator indices can again be obtained with the help of IBP relations. For example, for I (2) 21 and I (2) 11 , we have In the limit m = M (x = 1), we have

Three-point case
In the case of three-point diagrams, IBP yields the relation which, for example, allows us to obtain expressions for the integral I (3) 112 = I mM 0112 , or, in terms of a series solution, To obtain an expression for the I (3) 111 integral, we may consider the following, easy-to-obtain differential equation 3 : The solution of Eq. (21), with account of the corresponding boundary conditions at x = 0, 4 is given by where In terms of a series, the solution takes form In the limit m = M (x = 1), we have

Four-point case
The solution of the four-point case is simple. In this case, we may exploit the fact that the propagators of the integral I mM a 1 a 2 a 3 a 4 in Eq. (4) are linearly dependent in the special kinematics of Eq. (1), so that diagrams with four propagators are reduced to ones with at least one propagator less. Namely, we have the relation which can be used to obtain the following expressions: 3 As before, this follows from IBP relations and is valid up to O(ε) terms. 4 In the limit M 2 → ∞, the integral I 111 goes to zero. where

Two-loop diagrams
We are now ready to evaluate the two-loop diagrams J mM 010a 1 a 2 , J mM 011a 1 a 2 , and J mM 111a 1 a 2 , which can be obtained from J mM b 1 b 2 b 3 a 1 a 2 shown in Fig. 1(b). The latter integral has the expression where again the special kinematics of Eq. (1) where ã i = d/2 − a i and a = a 1 + a 2 − d/2. We note that, strictly speaking, in Eq. (30), we should use the one-loop integrals I mM ...,k+ε , where k = 0, 1, 2. On the other hand, in Section 3, we presented results for the I mM ...k integrals with O(1) accuracy only. However, the latter are sufficient to determine the most complicated contributions from the two-loop diagrams, containing series in n. The less complicated terms can be obtained directly from the two-loop diagrams, either by using asymptotic expansions or differential equations. Moreover, the knowledge of the Frobenius solutions of the differential equations together with the ansätze for the occurring sums is sufficient to determine the required expressions for the two-loop diagrams (see Appendix A).

Series representations
Let us start with the two-point sunset diagrams. It is known, e.g. from Ref. [96], that for the special kinematics in Eq. (1) there exist three master integrals. Choosing the latter to be J mM 01011 , J mM 01012 , and J mM 01022 and following the above general procedure, we obtain after some calculations 122 + J (2) 122 112 + J where Here, for brevity, we have omitted the arguments of the harmonic sums S 1 (n) and introduced the short-hand notations For the practical applications mentioned in Section 1, we also need the O(ε) terms of the sunset diagrams. We do not present them here for a general value of x, since the corresponding expressions are too lengthy. However, they may be found in Ref. [90].
Similarly, for the three-and four-point two-loop integrals, we have 5 where

Integral representations
Although the series expansions obtained in Subsection 4.1 are rapidly converging, it is useful to also find integral representations for the considered diagrams. To derive integral representations, it is helpful to rewrite the above series representations in the following form: Now, we write where J (1) Starting from this point, the corresponding integral representations can be obtained in two steps. First, we rewrite the functions (4n + 2a(n, δ)), with some a(n, δ) depending on the series under consideration, in the denominators as products of two functions, (2n + a(n, δ)) and (2n + a(n, δ) + 1/2). Second, we represent the ratio of the function (2n + b(n, δ)), with some b(n, δ) depending on the series under consideration, in the numerator and the function (2n + a(n, δ) + 1/2) in the denominator as (2n + b(n, δ)) (2n + a(n, δ) The remaining series can now be summed, and we are left with the sought integral representations over t. This way, we obtain the following integral representations for the two-point diagrams: where A = xt/4, Notice that the above integral representations for the sunset diagrams are simpler than those in Ref. [90], which involved some dilogarithms. Similarly, for the functions entering the expressions for the three-and four-point integrals, we have where I 1 (A) is defined in Eq. (23). Such correspondence reveals a deep relation between the subintegral expressions of the two-loop integrals and the corresponding one-loop results, in agreement with the general prescription in Eq. (30).
To check that the above integral representations are simplest, it is convenient to use the following relations between the functions J (i) ... 12 and J (i) where δ i 2 is the Kronecker symbol, which directly follow from the corresponding series representations.
Indeed, in the simplest forms of the integral representations, these relations can easily be reproduced. As a first example, let us consider Eq. (48) Now, notice that the expression in brackets on the r.h.s. of Eq. (51) only depends on A. So, we can replace the derivative with respect to x acting on it with Similarly, we have so that, on the r.h.s. of Eq. (51), we can make the replacement that is Next, using IBP, we have d dx The first term on the r.h.s. of Eq. (56) is equal to zero, and d dt So, finally, we have d dx Then, the evaluation of the terms proportional to A 2 on the r.h.s. of Eq. (59) gives 1 4x so that, finally, we have d dx where we have used Eq. (42) in the second equality. We thus observe that Eq. (48) with i = 1 is also correct.

Representations in terms of generalized hypergeometric functions
For the two-point sunset-type diagrams, the results in terms of hypergeometric functions were already given in Eq. (2). So, here we present such results for the three-and four-point diagrams only. The series representations obtained in Subsection 4.1 can be readily expressed in terms of hypergeometric functions and their derivatives. We thus obtain d dδ

Equal-mass case
In Ref. [90], we showed that, in the equal-mass case, the results for the two-loop sunset diagrams with the special kinematics in Eq. (1) can be written in terms of five "elliptic" sums, where φ = S 1 − 3S 1 + 2S 1 . Specifically, up to O(ε) terms, we have 6 Two of the "elliptic" sums above were shown to be expressible in terms of elliptic integrals of the first and second kinds. Namely, we have where 2 sin 2 ϕ = 1.8829167613 . . . , 6 In the equal-mass case, only two of the sunset diagrams correspond to master integrals, namely J mm 111 and J mm 112 .
with p = √ 5 and ϕ = arctan √ p. Here, F and E are the elliptic integrals of the first and second kinds, respectively. They are defined as However, to write equal-mass expressions for the three-and four-point functions, we need three additional "elliptic" sums, Then, the expressions for the three-and four-point functions are given by We are dealing here with incomplete elliptic integrals; see Section 6 of Ref. [97] for more details.

Conclusions
In this paper, we considered a class of two-loop diagrams with two different masses, m and M, in the special kinematic regime defined by Eq. (1), which play a crucial role in nonrelativistic effective field theories like NRQCD and NRQED. For all these diagrams, we presented analytic results in three different forms. First, we found explicit expressions for the coefficients of the series expansions in mass ratio x = m 2 /M 2 . These allow for the numerical restoration of the results at finite values of x using, for example, Padé approximants. Second, we provided integral representations. These should allow us to express the obtained results also in term of elliptic polylogarithms , which will be the subject of one of our subsequent publications. Third, we also found representations in terms of generalized hypergeometric functions, up to terms of O(ε). Following Ref. [94], we plan to find the exact results for the considered Feynman diagrams without performing ε expansions. For x = 1, our results could be expressed in terms of fast-converging, alternating series. The latter, in general, satisfy a set of relations, which can be found, for example, with the help of the PSLQ algorithm [98]. We were able to express some of these sums in terms of elliptic integrals of the first and second kinds. In general, we expect relations between our sums and elliptic multiple zeta values [41,[99][100][101][102][103][104], which we leave for future work. so that I 00a 3 a 4 a 5 a 6 a 7 = J mM a 3 a 4 a 5 a 6 a 7 .
Now, introducing the vector of master integrals, and using IBP, we obtain the following system of differential equations 7 : where x(x 2 +4) We seek for the solution of Eq. (75) using the Frobenius method as presented in Ref. [107]. 8 To this end, we first need to reduce Eq. (75) to Fuchsian form, i.e. to obtain an equivalent system that only has first-order poles in the variable x as x → 0. The required transformation matrix is given by 9 where 7 We used the program package LiteRed [105,106] for this purpose. 8 For other presentations of the Frobenius method appropriate to linear systems of differential equations, see Refs. [19,[108][109][110][111]. 9 It was found with the program package Fuchsia [112] and then checked with the program package Libra [113].
The algebraic multiplicity of eigenvalue λ is the number of times it shows up in the polynomial det(L 0 (λ)). Another characteristics associated with eigenvalue λ is its geometric multiplicity m g (λ). The latter is defined as the dimension of the eigenspace associated with λ such that L 0 (λ)v = 0 with n-dimensional eigenvector v. In our case, the geometric multiplicities are equal to the algebraic ones, which implies that we have as many Jordan blocks in matrix L 0 (λ) with zero eigenvalue as there are eigenvectors for λ, and each has size 1. Also, this implies that, for non-resonant eigenvalues, we have as many logarithm-free solutions as corresponds to the number of their geometric multiplicities. In general, one should consider root polynomials associated with the eigenvalues λ for each partial multiplicity and account for solutions containing logarithms (see Ref. [107] for more details). From the expressions for the eigenvalues λ in Eq. (84), we see that we have resonant eigenvalues λ = −1 + ε, −1/2 + ε, −2 + 3ε, and the solution presented above is not valid in these cases. In principle, we could use balance transformations [114] to get rid of the resonant eigenvalues in our system from the very beginning. Still, the method of Ref. [107] allows us to find the corresponding series solutions even in the case of resonances. Indeed, following Ref. [107], let us consider the resonant eigenvalue λ 1 such that Re λ 1 < · · · < Re λ r and λ i − λ 1 ∈ N * . Setting α = α(λ 1 ) = r i=2 m a (λ i ), 12 let us consider the modified non-homogeneous system, where g 0 (λ) is again an arbitrary n-dimensional, λ-dependent vector. Substituting the ansatz into Eq. (86), we obtain the following recurrence relation for the coefficients of the logarithm-free series solution: for k ≥ 1. The term proportional to (λ − λ 1 ) α on the right-hand side of Eq. (88) guarantees that we can always solve Eq. (88) for g k (λ) in the vicinity of λ = λ 1 and that the obtained solution is regular at λ = λ 1 . The found solution is, however, a multiple of a regular solution associated with the non-resonant eigenvalue λ r calculated before. To find linearly independent solutions, we first find m g (λ 1 ) independent eigenvectors g 0,i (λ 1 ) of the system L 0 (λ 1 )g 0,i (λ 1 ) = 0. Then, the sought solutions are given by the derivatives d α I(λ, x, g 0,i )/dλ α . 13