Two-loop diagrams in non-relativistic QCD with elliptics

TL;DR

This paper computes complex two-loop Feynman diagrams with elliptic structures in non-relativistic QCD, providing multiple representations and explicit results relevant for particle decay and production processes.

Contribution

It introduces new methods for evaluating two-loop diagrams with elliptic structures in NRQCD, including series, integral, and hypergeometric function representations.

Findings

Results valid up to ( ext{epsilon}) in 4-2 dimensions

Explicit series in elliptic constants for equal masses

Applicable to parapositronium decay and top pair production

Abstract

In this paper we consider two-loop two-, three- and four-point diagrams with elliptic structure in the case of two different masses $m$ and $M$. The latter diagrams generally arise within NRQCD matching procedures and are relevant for parapositronium decay and top pair production at threshold. We present the obtained results in several different representations: series solution with binomial coefficients, integral representation and representation in terms of generalized hypergeometric functions. The results are valid up to $\mathcal{O}(\varepsilon)$ terms in $d=4-2\varepsilon$ space-time dimensions. In the limit of equal masses $m=M$ the obtained results are written in terms of elliptic constants with explicit series representation.

Equations

$$q_1^2=q_2^2=0,p_1^2=p_2^2\equiv p^2=m^2,s=(q_1+q_2{)}^2=4m^2,t=(p-q_1{)}^2=-m^2.$$

$$\displaystyle e^{2\varepsilon\gamma_{E}}\left(M^{2}\right)^{2\varepsilon+1}J_{01022}^{mM}=_{4}F_{3}\left(\begin{array}[]{c}\frac{1}{2},1,1,1\\ \frac{3}{4},\frac{5}{4},\frac{3}{2}\end{array}\Bigg{|}-\frac{x^{2}}{4}\right)$$

$$e^{2\epsilon {\gamma }_E}{\left(M^2\right)}^{2\epsilon }{J}_{01012}^{mM}=-\frac{1}{2{\epsilon }^2}-\frac{1}{2\epsilon }-\frac{1}{2}(1+{\zeta }_2)+\frac{x}{4}(5-2\ln x)$$

$$e^{2\epsilon {\gamma }_E}{\left(M^2\right)}^{2\epsilon -1}{J}_{01011}^{mM}=-\frac{1}{{\epsilon }^2}(1+\frac{x}{2})+\frac{1}{\epsilon }\left[-3+\frac{1}{4}(-7+4\ln x)\right]-7-{\zeta }_2$$

$$\mbox{{ \begin{picture}(60.0,30.0)(0.0,13.0) \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \raise 27.0pt\hbox to0.0pt{\kern 30.0pt\makebox(0.0,0.0)[b]{$\scriptstyle a,\,M^{2}_{1}$}\hss} \ignorespaces \raise 3.0pt\hbox to0.0pt{\kern 30.0pt\makebox(0.0,0.0)[t]{$\scriptstyle b,\,M^{2}_{2}$}\hss} \ignorespaces \raise 12.0pt\hbox to0.0pt{\kern 0.0pt\makebox(0.0,0.0)[t]{$q$}\hss} \ignorespaces \end{picture} }}\hskip 5.69054pt~{}=~{}i^{1+d}\frac{\Gamma(a+b-d/2)}{\Gamma(a)\,\Gamma(b)}\int\limits_{0}^{1}\frac{{\rm d}s}{s^{a+1-d/2}(1-s)^{b+1-d/2}}\times\hskip 8.53581pt\raisebox{2.84526pt}{{ \begin{picture}(70.0,30.0)(0.0,4.0) \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \put(0.0,0.0){} \raise 7.0pt\hbox to0.0pt{\kern 33.0pt\makebox(0.0,0.0)[b]{$\scriptstyle a+b-d/2$}\hss} \ignorespaces \raise 7.0pt\hbox to0.0pt{\kern 33.0pt\makebox(0.0,0.0)[b]{}\hss} \ignorespaces \raise 1.0pt\hbox to0.0pt{\kern 33.0pt\makebox(0.0,0.0)[t]{$\scriptstyle M^{2}_{1}/s+M^{2}_{2}/(1-s)$}\hss} \ignorespaces \raise-5.0pt\hbox to0.0pt{\kern-3.0pt\makebox(0.0,0.0)[b]{$q$}\hss} \ignorespaces \end{picture} }}\hskip 8.53581pt,$$

$$I_{a_1{a}_2{a}_3{a}_4}^{mmmM}$$

$$I_{0101}^{mM}(d-3)\left(M^2-2m^2\right)$$

$$I_{0011}^{mM}(d-3)M^2$$

$$I_{0102}^{mM}=\frac{d}{d{M}^2}{I}_{0101}^{mM},I_{0012}^{mM}=\frac{d}{d{M}^2}{I}_{0011}^{mM}.$$

$$I_{000a}^M=\frac{i(-1{)}^a}{(M^2{)}^{a-d/2}}\frac{\Gamma (a-d/2)}{\Gamma (a)},I_{00a0}^m=\frac{i(-1{)}^a}{(m^2{)}^{a-d/2}}\frac{\Gamma (a-d/2)}{\Gamma (a)}.$$

$$(1-2x)(1-4x)\frac{d}{dx}{I}_{12}^{(1)}-4x{I}_{12}^{(1)}+i(1-2x+2x\ln x)=0,$$

$$\frac{M^2}{i}{I}_{12}^{(1)}=\frac{\ln x}{2x}+\frac{(1-2x)}{2x\sqrt{1-4x}}\ln \frac{1+\sqrt{1-4x}}{1-\sqrt{1-4x}}+\mathcal{O}(\epsilon ),$$

$$\frac{M^{2}}{i}I_{12}^{(1)}=-1-\sum_{n=1}^{\infty}\left(\begin{array}[]{c}2n\\ n\end{array}\right)\frac{n}{n+1}x^{n}\left[\ln x-2S_{1}(n-1)+2S_{1}(2n-1)-\frac{1}{n+1}\right]+\mathcal{O}(\varepsilon)\,.$$

$$\frac{M^2}{i}{I}_{21}^{(1)}$$

$$e^{\epsilon {\gamma }_E}\frac{(M^2{)}^{\epsilon }}{i}{I}_{11}^{(1)}$$

$$\frac{M^2}{i}{I}_{12}^{(1)}$$

$$e^{\epsilon {\gamma }_E}\frac{(M^2{)}^{\epsilon }}{i}{I}_{11}^{(1)}$$

$$\left(1+4x^2\right)\frac{d}{dx}{I}_{12}^{(2)}+4x{I}_{12}^{(2)}+i(1+2x+2x\ln x)=0,$$

$$\frac{M^2}{i}{I}_{12}^{(2)}=-\frac{1}{2x}\ln x+\frac{1}{4x\sqrt{1+4x^2}}\left(\ln \frac{\sqrt{1+4x^2}-1}{\sqrt{1+4x^2}+1}+\ln \frac{\sqrt{1+4x^2}-2x}{\sqrt{1+4x^2}+2x}\right)+\mathcal{O}(\epsilon ),$$

$$\frac{M^2}{i}{I}_{12}^{(2)}$$

$$\frac{M^2}{i}{I}_{21}^{(2)}$$

$$e^{\epsilon {\gamma }_E}\frac{(M^2{)}^{\epsilon }}{i}{I}_{11}^{(2)}$$

$$\frac{M^2}{i}{I}_{12}^{(2)}$$

$$e^{\epsilon {\gamma }_E}\frac{(M^2{)}^{\epsilon }}{i}{I}_{11}^{(2)}$$

$$I_{011a}^{mM}=\frac{1}{2m^2(a-1)}\left[(a-1)I_{001a}+I_{002a-1}-(a-1)I_{010a}-I_{020a-1}\right],$$

$$\frac{(M^2{)}^2}{i}{I}_{112}^{(3)}$$

$$\frac{d}{dx}{I}_{111}^{(3)}+\frac{1}{x(1-2x)}{I}_{12}^{(1)}+\frac{1}{x}{I}_{12}^{(2)}-\frac{i\ln x}{1-2x}=0.$$

$$\frac{M^2}{i}{I}_{111}^{(3)}=\frac{1}{4x}\left[{\ln }^2\frac{1+\sqrt{1-4x}}{1-\sqrt{1-4x}}-\frac{1}{4}{\ln }^2\frac{\sqrt{1+4x^2}-1}{\sqrt{1+4x^2}+1}-2I_1(x)\right]+\mathcal{O}(\epsilon ),$$

$$I_1(x)={\text{Li}}_2(-R(x))-{\mathrm{Li}}_2(R(x)),R(x)=\frac{1-\sqrt{S(x)}}{1+\sqrt{S(x)}},S(x)=\frac{\sqrt{1+4x^2}-2x}{\sqrt{1+4x^2}+2x}.$$

$$\frac{M^2}{i}{I}_{111}^{(3)}$$

$$\frac{(M^2{)}^2}{i}{I}_{112}^{(3)}$$

Full text

PACS: 12.38.Bx, 13.20.Gd, 13.25.Gv, 14.40.Pq

1 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\varepsilon$ 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 $\varepsilon$.222At order $\mathcal{O}(\varepsilon)$, this was already stated in Ref. [31].

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 $\varepsilon$ 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 [37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57]. 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\to 2$ processes with two real photons ($\gamma\gamma$) or gluons ($gg$) in the initial state and two on-shell massive fermions ($f\bar{f}$ ) in the final state. In addition, we stick to threshold kinematics. For example, we have $\gamma(q_{1})+\gamma(q_{2})\to f(p_{1})+\bar{f}(p_{2})$, where

[TABLE]

The physical motivation to study such processes with threshold kinematics is the following. Such kinematics appears, 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 $t\bar{t}$ 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 one-loop 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 a by 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\varepsilon$ 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 $\mathcal{O}(\varepsilon)$ terms and represented as fast converging series in the mass ratio $x=m^{2}/M^{2}$. The $\mathcal{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

[TABLE]

where $\zeta_{n}=\zeta(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 $\varepsilon$ expansion. In the present paper, we extend the results of Ref. [90] for the case of three- and four-point diagrams (see Fig. 1) with $\mathcal{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. (19). 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\neq 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.

2 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:

[TABLE]

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 (19) 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. (19), 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 one-loop 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. (19). A similar strategy was also adopted for the calculation of certain four-loop tadpole diagrams in Ref. [95].

3 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}$:

[TABLE]

where the kinematics in Eq. (1) is implied. Using IBP, the integral in Eq. (20) 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.

Following the line numbering in Fig. 2(c), the diagrams shown in Figs. 2(a)–(c) correspond to the integrals $I^{mmmM}_{0011}$, $I^{mmmM}_{0111}$, and $I^{mmmM}_{1111}$, 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}$.

3.1 Two-point case

In the case of the two-point diagrams, IBP provides us with the relations

[TABLE]

which can be considered as the differential equations for the initial diagrams $I_{0101}^{mM}$ and $I_{0011}^{mM}$, as

[TABLE]

The integrals $I_{000a}^{M}$ and $I_{00a0}^{m}$ in Eq. (21) are tadpoles with masses $M$ and $m$, respectively,

[TABLE]

Similarly, using IBP and substituting $m=1$ and ${d}/{dM^{2}}=-x^{2}{d}/{dx}$, we can deduce a differential equation for $I_{12}^{(1)}=I_{0102}^{mM}$,

[TABLE]

which is valid up to $\mathcal{O}(\varepsilon)$ corrections. Equation (24) is easy to solve and, taking the boundary conditions at $x=0$ into account, we have

[TABLE]

or, in terms of a series solution,

[TABLE]

Now, IBP relations allow us to write down answers for the $I_{a_{1}a_{2}}^{(1)}$ integrals with other values of propagator indices. For example, we have

[TABLE]

For $m=M$ ($x=1$), we have

[TABLE]

In the case of the integral $I_{12}^{(2)}=I_{0012}^{mM}$, the application of IBP relations leads to the following differential equation:

[TABLE]

which is again valid up to $\mathcal{O}(\varepsilon)$ terms. The solution of Eq. (29) is straightforward and, taking the boundary conditions at $x=0$ into account, yields

[TABLE]

or, in terms of a series solution,

[TABLE]

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_{21}^{(2)}$ and $I_{11}^{(2)}$, we have

[TABLE]

In the limit $m=M$ ($x=1$), we have

[TABLE]

3.2 Three-point case

In the case of three-point diagrams, IBP yields the relation

[TABLE]

which, for example, allows us to obtain expressions for the integral $I_{112}^{(3)}=I_{0112}^{mM}$,

[TABLE]

or, in terms of a series solution,

[TABLE]

To obtain an expression for the $I_{111}^{(3)}$ integral, we may consider the following, easy-to-obtain differential equation:333As before, this follows from IBP relations and is valid up to $\mathcal{O}(\varepsilon)$ terms.

[TABLE]

The solution of Eq. (45), with account of the corresponding boundary conditions at $x=0$,444In the limit $M^{2}\to\infty$, the integral $I_{111}^{(3)}$ goes to zero. is given by

[TABLE]

where

[TABLE]

In terms of a series, the solution takes form

[TABLE]

In the limit $m=M$ ($x=1$), we have

[TABLE]

3.3 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_{a_{1}a_{2}a_{3}a_{4}}^{mM}$ in Eq. (20) are linearly dependent in the special kinematics of Eq. (1) and so that diagrams with four propagators are reduced to ones with at least one propagator less. Namely, we have the relation

[TABLE]

which can be used to obtain the following expressions:

[TABLE]

where

[TABLE]

4 Two-loop diagrams

We are now ready to evaluate the two-loop diagrams $J_{010a_{1}a_{2}}^{mM}$, $J_{011a_{1}a_{2}}^{mM}$, and $J_{111a_{1}a_{2}}^{mM}$, which can be obtained from $J_{b_{1}b_{2}b_{3}a_{1}a_{2}}^{mM}$ shown in Fig. 1(b). The latter integral has the expression

[TABLE]

where again the special kinematics of Eq. (1) is implied. Exploiting the rule of Eq. (19), we replace the $k_{1}$ integration in Eq. (56) with the integration over the variable $s$ from the corresponding one-loop integral $I_{...a}^{mM}$ with the replacement $M^{2}\to\mu^{2}=M^{2}/[s(1-s)]$. Indeed, calling the above two-loop integral $J_{\dots a_{1}a_{2}}^{mM}$ and the corresponding one-loop diagram $I_{\dots a}^{mM}$, we have

[TABLE]

where $\tilde{a}_{i}=d/2-a_{i}$ and $a=a_{1}+a_{2}-d/2$. We note that, strictly speaking, in Eq. (57), we should use the one-loop integrals $I_{\dots,k+\varepsilon}^{mM}$, where $k=0,1,2$. On the other hand, in Section 3, we presented results for the $I_{\dots k}^{mM}$ integrals with $\mathcal{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).

4.1 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_{01011}^{mM}$, $J_{01012}^{mM}$, and $J_{01022}^{mM}$ and following the above general procedure, we obtain after some calculations

[TABLE]

where

[TABLE]

Here, for brevity, we have omitted the arguments of the harmonic sums $S_{1}(n)$ and introduced the short-hand notations

[TABLE]

For the practical applications mentioned in Section 1, we also need the $\mathcal{O}(\varepsilon)$ 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 have555Here, we have summed series corresponding to multiple polylogarithms.

[TABLE]

where

[TABLE]

4.2 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:

[TABLE]

Now, we write

[TABLE]

where

[TABLE]

Starting from this point, the corresponding integral representations can be obtained in two steps. First, we rewrite the functions $\Gamma(4n+2a(n,\delta))$, with some $a(n,\delta)$ depending on the series under consideration, in the denominators as products of two $\Gamma$ functions, $\Gamma(2n+a(n,\delta))$ and $\Gamma(2n+a(n,\delta)+1/2)$. Second, we represent the ratio of the function $\Gamma(2n+b(n,\delta))$, with some $b(n,\delta)$ depending on the series under consideration, in the numerator and the function $\Gamma(2n+a(n,\delta)+1/2)$ in the denominator as

[TABLE]

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:

[TABLE]

where $A={xt}/{4}$,

[TABLE]

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

[TABLE]

where $I_{1}(A)$ is defined in Eq. (47). 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. (57).

To check that the above integral representations are simplest, it is convenient to use the following relations between the functions $J_{...12}^{(i)}$ and $J_{...11}^{(i)}$ $(i=1,2)$:

[TABLE]

where $\delta^{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. (143) with $i=2$. From Eq. (140), we have

[TABLE]

Now, notice that the expression in brackets on the r.h.s. of Eq. (146) only depends on $A$. So, we can replace the derivative with respect to $x$ acting on it with

[TABLE]

Similarly, we have

[TABLE]

so that, on the r.h.s. of Eq. (146), we can make the replacement

[TABLE]

that is

[TABLE]

Next, using IBP, we have

[TABLE]

The first term on the r.h.s. of Eq. (151) is equal to zero, and

[TABLE]

So, finally, we have

[TABLE]

where we have used Eq. (138) in the third equality. We thus observe that Eq. (143) with $i=2$ is correct.

As a second example, let us consider Eq. (143) with $i=1$. Repeating steps similar to those in Eqs. (146)–(153), we have

[TABLE]

Then, the evaluation of the terms proportional to $A^{2}$ on the r.h.s. of Eq. (154) gives

[TABLE]

so that, finally, we have

[TABLE]

where we have used Eq. (137) in the second equality. We thus observe that Eq. (143) with $i=1$ is also correct.

4.3 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. (18). 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

[TABLE]

4.4 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,

[TABLE]

where $\phi=S_{1}-3\overline{S_{1}}+2\overline{\overline{S_{1}}}$. Specifically, up to $\mathcal{O}(\varepsilon)$ terms, we have666In the equal-mass case, only two of the sunset diagrams correspond to master integrals, namely $J_{111}^{mm}$ and $J_{112}^{mm}$.

[TABLE]

where

[TABLE]

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

[TABLE]

where

[TABLE]

with $p=\sqrt{5}$ and $\varphi=\arctan\sqrt{p}$. Here, $F$ and $E$ are the elliptic integrals of the first and second kinds, respectively. They are defined as

[TABLE]

However, to write equal-mass expressions for the three- and four-point functions, we need three additional “elliptic” sums,

[TABLE]

Then, the expressions for the three- and four-point functions are given by

[TABLE]

where

[TABLE]

We are dealing here with incomplete elliptic integrals; see Section 6 of Ref. [97] for more details.

5 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 [37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57], 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 $\mathcal{O}(\varepsilon)$. Following Ref. [94], we plan to find the exact results for the considered Feynman diagrams without performing $\varepsilon$ 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.

Acknowledgments

The authors are grateful to M. Yu. Kalmykov and R. N. Lee for fruitful discussions and valuable comments. The work of B.A.K. and O.L.V. was supported in part by the German Federal Ministry for Education and Research (BMBF) through Grant No. 05H18GUCC1 and by the German Research Foundation (DFG) through Grant No. KN 365/12-1. The work of A.V.K. was supported in part by the Russian Foundation for Basic Research (RFBR) through Grant No. 17-02-00872, by the Russian Ministry of Science and Higher Education under Contract No. 02.A03.21.0003 as of 27.08.2013, and by the Heisenberg-Landau program. The work of A.I.O. was supported in part by the RFBR through Grant No. 17-02-00872 and by the Foundation for the Advancement of Theoretical Physics and Mathematics “BASIS”.

Appendix A Frobenius solution for a system of differential equations

To use IBP and write down a system of differential equations for the integrals under consideration, it is convenient to introduce somewhat more general two-loop integrals,

[TABLE]

so that

[TABLE]

Now, introducing the vector of master integrals,

[TABLE]

and using IBP, we obtain the following system of differential equations:777We used the program package LiteRed [105, 106] for this purpose.

[TABLE]

where

[TABLE]

We seek for the solution of Eq. (193) using the Frobenius method as presented in Ref. [107].888For other presentations of the Frobenius method appropriate to linear systems of differential equations, see Refs. [19, 108, 109, 110, 111]. To this end, we first need to reduce Eq. (193) to Fuchsian form, i.e. to obtain an equivalent system that only has first-order poles in the variable $x$ as $x\to 0$. The required transformation matrix is given by999It was found with the program package Fuchsia [112] and then checked with the program package Libra [113].

[TABLE]

where

\displaystyle\mathbf{T}=\left(\begin{array}[]{cccccccccc}1&0&0&0&0&0&0&0&0&0\\ 0&1&0&0&0&0&0&0&0&0\\ 0&0&-\frac{-2i-(1-2i)x+x^{2}}{x}&0&0&0&0&0&0&0\\ 0&0&0&1-x&0&0&0&0&0&0\\ 0&0&0&0&1-x&0&0&0&0&0\\ 0&0&0&0&0&1-x&0&0&0&0\\ 0&0&0&0&0&0&-\frac{-2i-(1-2i)x+x^{2}}{x}&0&0&0\\ \frac{3(1-2\varepsilon+\varepsilon^{2})x}{4(2\varepsilon-1)}&-\frac{(1-2\varepsilon+\varepsilon^{2})x}{2(2\varepsilon-1)}&0&0&0&0&0&1-x&0&0\\ 0&0&0&0&0&0&-\frac{-2i+2ix+x^{2}}{x}&0&1&0\\ 0&0&0&0&0&0&0&0&0&1\end{array}\right)\,.

(216)

For the reader’s convenience, we use the notation of Ref. [107] in what follows. Then, the system we wish to solve is rewritten as101010When passing to the representation of the system in terms of the $\mathbf{A}_{i}$ matrices, we multiply our original system by the common denominator of the $\mathbf{\widetilde{M}}$ matrix, so that the $\mathbf{A}_{i}$ matrices become polynomials in $x$.

[TABLE]

where $\mathbf{A}_{i}(x)=\sum_{j=0}^{\infty}\mathbf{A}_{i,j}x^{j}\in\mathbb{C}[x]^{n\times n}$ with $i=0,1$ and $n=10$, and $\mathbf{A}_{i,j}$ are equal to zero for sufficiently large values of $j$.

Now, following Barkatou and Gluzeau [107] and starting from the system in Eq. (218), we consider a more general inhomogeneous system,

[TABLE]

where $\mathbf{g}_{0}(\lambda)$ is an arbitrary $n$-dimensional, $\lambda$-dependent vector and we have introduced the following matrix polynomials:

[TABLE]

The reduction to Fuchsian form performed before was necessary in order to have a regular matrix polynomial $\mathbf{L}_{0}(\lambda)$, such that $\det(\mathbf{L}_{0}(\lambda))\neq 0$. To find the logarithm-free series solution, we make the ansatz

[TABLE]

where $\mathbf{g}_{i}(\lambda)$ are arbitrary $n$-dimensional, $\lambda$-dependent vectors, and substitute it into Eq. (219). Extracting coefficients at equal powers of $x$ on both sides of the resulting equation, we obtain

[TABLE]

for $i\geq 1$. If $\lambda=\lambda_{1}$ and $\mathbf{g}_{0}(\lambda_{1})$ is such that $\det(\mathbf{L}_{0}(\lambda_{1}))=0$, $\mathbf{L}_{0}(\lambda_{1})\mathbf{g}_{0}(\lambda_{1})=\mathbf{0}$, and $\det(\mathbf{L}_{0}(\lambda_{1}+i))\neq 0$ for $i\geq 1$, then the found solution is also the solution of our initial homogeneous equation in Eq. (218). The last condition, $\det(\mathbf{L}_{0}(\lambda_{1}+i))\neq 0$, implies that $\lambda_{1}$ is a non-resonant eigenvalue of matrix $\mathbf{L}_{0}(\lambda)$, and we can solve Eq. (222) for $\mathbf{g}_{i}(\lambda_{1})$. For our system, we have seven different eigenvalues111111By eigenvalues, we mean solutions of the equation $\det(\mathbf{L}_{0}(\lambda))=0$. of matrix $\mathbf{L}_{0}(\lambda)$,

[TABLE]

with algebraic multiplicities $m_{a}(\lambda)$,

[TABLE]

The algebraic multiplicity of eigenvalue $\lambda$ is the number of times it shows up in the polynomial $\det(\mathbf{L}_{0}(\lambda))$. Another characteristics associated with eigenvalue $\lambda$ is its geometric multiplicity $m_{g}(\lambda)$. The latter is defined as the dimension of the eigenspace associated with $\lambda$ such that $\mathbf{L}_{0}(\lambda)\mathbf{v}=\mathbf{0}$ with $n$-dimensional eigenvector $\mathbf{v}$. In our case, the geometric multiplicities are equal to the algebraic ones, which implies that we have as many Jordan blocks in matrix $\mathbf{L}_{0}(\lambda)$ with zero eigenvalue as there are eigenvectors for $\lambda$, 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 $\lambda$ for each partial multiplicity and account for solutions containing logarithms (see Ref. [107] for more details).

From the expressions for the eigenvalues $\lambda$ in Eq. (223), we see that we have resonant eigenvalues $\lambda=-1+\varepsilon,-1/2+\varepsilon,-2+3\varepsilon$, 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 $\lambda_{1}$ such that $\operatorname{Re}\lambda_{1}<\cdots<\operatorname{Re}\lambda_{r}$ and $\lambda_{i}-\lambda_{1}\in\mathbb{N}^{*}$. Setting $\alpha=\alpha(\lambda_{1})=\sum_{i=2}^{r}m_{a}(\lambda_{i})$,121212In our case, all values of $\alpha$ for resonant eigenvalues are equal to unity. let us consider the modified non-homogeneous system,

[TABLE]

where $\mathbf{g}_{0}(\lambda)$ is again an arbitrary $n$-dimensional, $\lambda$-dependent vector. Substituting the ansatz

[TABLE]

into Eq. (225), we obtain the following recurrence relation for the coefficients of the logarithm-free series solution:

[TABLE]

for $k\geq 1$. The term proportional to $(\lambda-\lambda_{1})^{\alpha}$ on the right-hand side of Eq. (227) guarantees that we can always solve Eq. (227) for $\mathbf{g}_{k}(\lambda)$ in the vicinity of $\lambda=\lambda_{1}$ and that the obtained solution is regular at $\lambda=\lambda_{1}$. The found solution is, however, a multiple of a regular solution associated with the non-resonant eigenvalue $\lambda_{r}$ calculated before. To find linearly independent solutions, we first find $m_{g}(\lambda_{1})$ independent eigenvectors $\mathbf{g}_{0,i}(\lambda_{1})$ of the system $\mathbf{L}_{0}(\lambda_{1})\mathbf{g}_{0,i}(\lambda_{1})=\mathbf{0}$. Then, the sought solutions are given by the derivatives $d^{\alpha}\mathbf{\widetilde{I}}(\lambda,x,\mathbf{g}_{0,i})/d\lambda^{\alpha}$.131313This is only valid in the case $m_{a}(\lambda_{1})=m_{g}(\lambda_{1})$. See Ref. [107] for a discussion of the general situation.

In this way, we obtain $n=10$ linearly independent solutions of the system in Eq. (218) altogether. To fix ten constants in the general solution, we account for boundary conditions at $x=0$. The latter are obtained using asymptotic expansions of the Feynman diagrams at $x=0$, properly taking subgraphs into account, and read

[TABLE]

Having found the Frobenius solution for the vector of master integrals in Eq. (192), we may use the known basis of occurring sums to reconstruct the series solutions presented in Subsection 4.1.