Iterated Elliptic and Hypergeometric Integrals for Feynman Diagrams

TL;DR

This paper develops advanced mathematical techniques involving elliptic and hypergeometric integrals to compute complex 3-loop Feynman integrals in quantum field theory, introducing new classes of functions and solutions.

Contribution

It introduces a novel class of iterative integrals with elliptic and modular function representations for solving complex Feynman integrals, extending beyond known polylogarithmic functions.

Findings

Solutions expressed in terms of hypergeometric functions and elliptic integrals.

New iterative integrals with non-traditional letters are formulated.

Numerical results demonstrate the effectiveness of the methods.

Abstract

We calculate 3-loop master integrals for heavy quark correlators and the 3-loop QCD corrections to the $\rho$-parameter. They obey non-factorizing differential equations of second order with more than three singularities, which cannot be factorized in Mellin-$N$ space either. The solution of the homogeneous equations is possible in terms of convergent close integer power series as $_2F_1$ Gau\ss{} hypergeometric functions at rational argument. In some cases, integrals of this type can be mapped to complete elliptic integrals at rational argument. This class of functions appears to be the next one arising in the calculation of more complicated Feynman integrals following the harmonic polylogarithms, generalized polylogarithms, cyclotomic harmonic polylogarithms, square-root valued iterated integrals, and combinations thereof, which appear in simpler cases. The inhomogeneous solution of the corresponding differential equations can be given in terms of iterative integrals, where the new innermost letter itself is not an iterative integral. A new class of iterative integrals is introduced containing letters in which (multiple) definite integrals appear as factors. For the elliptic case, we also derive the solution in terms of integrals over modular functions and also modular forms, using $q$-product and series representations implied by Jacobi's $\vartheta_i$ functions and Dedekind's $\eta$-function. The corresponding representations can be traced back to polynomials out of Lambert--Eisenstein series, having representations also as elliptic polylogarithms, a $q$-factorial $1/\eta^k(\tau)$, logarithms and polylogarithms of $q$ and their $q$-integrals. Due to the specific form of the physical variable $x(q)$ for different processes, different representations do usually appear. Numerical results are also presented.

Equations

$$\mathbf{M}[f(x)](N)={\int }_0^{1}dx{x}^{N}f(x),$$

$$q=\exp (i\pi \tau ),$$

$$\tau =i\frac{\mathbf{K}(1-z(x))}{\mathbf{K}(z(x))}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{with}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\tau \in \mathbb{H}=\left\{z\in \mathbb{C},\mathsf{Im}(z)>0\right\}$$

$$\frac{1}{{\eta }^k(\tau )},\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}k>0,\text{leavevmode}\text{leavevmode}\text{leavevmode}k\in \mathbb{N}.$$

$$J=\frac{G_2^3(q)}{216000\Delta (q)}$$

$$\Delta (q)=(2\pi {)}^{12}{q}^2\prod _{k=1}^{\infty }(1-q^{2k}{)}^{24}=(2\pi {)}^{12}{\eta }^{24}(\tau ).$$

$$\left[\frac{d^2}{d{x}^2}+p(x)\frac{d}{dx}+q(x)\right]\psi (x)=N(x)\text{leavevmode},$$

$$r(x)=\sum _{k=1}^{n_r}\frac{b_k^{(r)}}{x-a_k^{(r)}},\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}{a}_k^{(r)},b_k^{(r)}\in \mathbb{Z}\text{leavevmode}.$$

$$W(x)={\psi }_1^{(0)}(x)\frac{d}{dx}{\psi }_2^{(0)}(x)-{\psi }_2^{(0)}(x)\frac{d}{dx}{\psi }_1^{(0)}(x)\ne 0\text{leavevmode}.$$

$$\frac{d}{dx}W(x)=-p(x)W(x)\text{leavevmode},$$

$$W(x)=\prod _{k=1}^{n_1}{\left(\frac{1}{x-a_k^{(1)}}\right)}^{b_k^{(1)}}\text{leavevmode},$$

$$\psi (x)$$

$$n(x)=\frac{N(x)}{W(x)}$$

$$\displaystyle\frac{d}{dx}\left(\begin{array}[]{c}f_{8a}(x)\\ f_{9a}(x)\end{array}\right)=\left(\begin{array}[]{cc}\frac{4}{x}&\frac{6}{x}\\ \frac{4(x^{2}-3)}{x(x^{2}-9)(x^{2}-1)}&\frac{2(x^{4}-9)}{x(x^{2}-9)(x^{2}-1)}\end{array}\right)\otimes\left(\begin{array}[]{c}f_{8a}(x)\\ f_{9a}(x)\end{array}\right)+\left(\begin{array}[]{c}N_{8a}(x)\\ N_{9a}(x)\end{array}\right)$$

$$\displaystyle\frac{d}{dx}\left(\begin{array}[]{c}f_{8b}(x)\\ f_{9b}(x)\end{array}\right)=\left(\begin{array}[]{cc}\frac{4}{x}&\frac{2}{x}\\ \frac{4(3x^{2}-1)}{x(9x^{2}-1)(x^{2}-1)}&\frac{2(9x^{4}-1)}{x(9x^{2}-1)(x^{2}-1)}\end{array}\right)\otimes\left(\begin{array}[]{c}f_{8b}(x)\\ f_{9b}(x)\end{array}\right)+\left(\begin{array}[]{c}N_{8b}(x)\\ N_{9b}(x)\end{array}\right),$$

$$N_{8a}(x)$$

$$N_{9a}(x)$$

$$N_{8b}(x)$$

$$N_{9b}(x)$$

$$0$$

$$f_{9a}(x)$$

$$0$$

$$f_{9b}(x)$$

$$\frac{d}{dx}{f}_{10a}(x)$$

$$\frac{d}{dx}{f}_{10b}(x)$$

$$H_{b,\vec{a}}(x)$$

$$H_0,H_1,H_{-1},H_{0,-1},H_{0,1},H_{0,0,-1},H_{0,0,1},H_{0,-1,-1},H_{0,-1,1},H_{0,1,-1},H_{0,1,1},$$

$${\psi }_{1a}^{(0)}(x)$$

$${\psi }_{2a}^{(0)}(x)$$

$$z=\frac{x^2(x^2-9{)}^2}{(x^2+3{)}^3}\text{leavevmode}.$$

$$b\sum _{k=0}^{\infty }{\tau }_k(c\cdot z{)}^k=\sum _{k=0}^{\infty }{m}_k{z}^k,\text{leavevmode}\text{with}\text{leavevmode}\text{leavevmode}{\tau }_k,b\in \mathbb{Q},\text{leavevmode}{m}_k\in \mathbb{Z},$$

$$W(x)=x(9-x^2)(x^2-1).$$

$${\psi }_{1b}^{(0)}(x)$$

$${\psi }_{2b}^{(0)}(x)$$

$$z(x)=-\frac{16x^3}{(x+1)(x-3{)}^3}\text{leavevmode}.$$

Full text

1 Introduction

Many single scale Feynman integrals arising in massless and massive multi-loop calculations in Quantum Chromodynamics (QCD) [QCD] have been found to be expressible in terms of harmonic polylogarithms (HPLs) [Remiddi:1999ew], generalized harmonic polylogarithms [Moch:2001zr, Ablinger:2013cf], cyclotomic harmonic polylogarithms [Ablinger:2011te], square-root valued iterated integrals [Ablinger:2014bra], as well as more general functions, entering the corresponding alphabet in integral iteration. After taking a Mellin transform

[TABLE]

they can be equivalently expressed in terms of harmonic sums [Blumlein:1998if, Vermaseren:1998uu] in the simpler examples and finite sums of different kinds in the other cases [Moch:2001zr, Ablinger:2013cf, Ablinger:2011te, Ablinger:2014bra], supplemented by special numbers like the multiple zeta values [Blumlein:2009cf] and others appearing in the limit $N\rightarrow\infty$ of the nested sums, or the value at $x=1$ of the iterated integrals in Refs. [Blumlein:1998if, Vermaseren:1998uu, Remiddi:1999ew, Moch:2001zr, Ablinger:2013cf, Ablinger:2011te, Ablinger:2014bra].

In many higher order calculations a considerable reduction of the number of integrals to be calculated is obtained using integration by parts identities (IBPs) [IBP], which allow to express all required integrals in terms of a much smaller set of so called master integrals. Differential equations satisfied by these master integrals [DEQ, Caffo:1998du] can then be obtained by taking their derivatives with respect to the parameters of the problem and inserting the IBPs in the result. What remains is to solve these differential equations, given initial or boundary conditions, if possible analytically. One way of doing this is to derive an associated system of difference equations [DIFFEQ, Ablinger:2010ty, Ablinger:2014nga, Ablinger:2015tua] after applying a mapping through a formal Taylor series or a Mellin transform. If these equations factorize to first order equations, we can use the algorithm presented in Ref. [Ablinger:2015tua] for general bases to solve these systems analytically and to find the corresponding alphabets over which the iterated integrals or nested sums are built. The final solution in $N$ and $x$ space is found by using the packages Sigma [SIG1, SIG2], EvaluateMultiSums and SumProduction [EMSSP].

However, there are physical cases where full first order factorizations cannot be obtained for either the differential equations in $x$ or the difference equations in $N$.111There are cases in which factorization fails in either $x$- or $N$-space, but not in both, cf. [Ablinger:2013eba]. This opportunity has to be always checked. The next level of complexity is given by non-factorizable differential or difference equations of second order. Examples of this are the massive sunrise and kite integrals [SABRY, TLSR1, Caffo:1998du, Caffo:2002ch, Laporta:2004rb, TLSR2, TLSR2a, Bailey:2008ib, Broadhurst:2008mx, TLSR3, BLOCH2, TLSR4, Adams:2014vja, Adams:2015gva, Adams:2015ydq, TLSR5, Remiddi:2016gno, Adams:2016xah, TSLR6].

In the present paper, we will address the analytic solution of typical cases of this kind, related to a series of master integrals appearing in the 3-loop corrections of the $\rho$-parameter in [Grigo:2012ji]. It turns out that these integrals are more general than those appearing in the sunrise and kite diagrams, due to the appearance of also the elliptic integral of the second kind, ${\bf E}(z)$, which cannot be transformed away. The corresponding second order differential equations have more than three singularities, as in the case of the Heun equation [HEUN]. For the sake of generality, we will seek solutions of the second order homogeneous differential equations which are given in terms of Gauß’ ${}_{2}F_{1}$ functions [GAUSS] within the class of globally bounded solutions [VANH1], cf. also [BOSTAN]. Here the parameters of the ${}_{2}F_{1}$ function are rational numbers and the argument is a rational function of $x$. The complete elliptic integrals ${\bf K}(z)$ and ${\bf E}(z)$ [TRICOMI, ELLINT, ABEL, JAC1]222As a convention, the modulus $k^{2}=z$ is chosen in this paper, also used within the framework of Mathematica. are special cases of this class.

The hypergeometric function obeys different relations like the Euler- and Pfaff-transformations [HYP, SLATER], the 24 Kummer solutions [KUMMER, RIEMANN] and the 15 Gauß’ contiguous relations [HYP, SLATER]. There are more special transformations for higher than first order in the argument [KUMMER, GOURSAT, VANHVID, VID2]. In the present case, equivalent ${}_{2}F_{1}$ representations are obtained by applying arithmetic triangle groups [TAKEUCHI]. The corresponding algorithm has been described in Ref. [IVH] in its present most far reaching form. The relations of this type may be useful to transform a found solution into another one, which might be particularly convenient. In the case a function space of more solutions is considered, these relations have to be exploited to check the independence of the basis elements.

The main idea of the approach presented here is to obtain the factorization of a high order scalar difference or differential equation, after uncoupling [Zuercher:94, ORESYS, NewUncouplingMethod] the corresponding linear systems, to all first order parts and its second order contributions. While the first order parts have been algorithmically solved in Ref. [Ablinger:2015tua], the treatment of second order differential equations shall be automated.333For more involved physical problems also irreducible higher order differential equations may occur. The class of ${}_{2}F_{1}$ solutions has an algorithmic automation to a wide extent [IVH] and it seems that this class constitutes the next one following the iterated square-root valued letters in massive single-scale 3-loop integrals. Applying this method, we obtain the corresponding ${{}_{2}F_{1}}$ functions with (partly) fractional parameters and rational argument, and ir(rational) pre-factors, forming the new letters of the otherwise iterated integrals. These letters contain a definite integral by virtue of the integral representation of the ${{}_{2}F_{1}}$ function, which cannot be fully transformed into an integral depending on the follow-up integration variable only through its integration boundaries. In general, we have therefore to iterate new letters of this kind. Through this we obtain a complete algorithmic automation of the solution also when second order differential operators contribute, having ${}_{2}F_{1}$ solutions.

As it will be shown, in a series of cases the reduction of the ${}_{2}F_{1}$ functions to complete elliptic integrals ${\bf E}(r(z))$ and ${\bf K}(r(z))$ is possible. Therefore we also study special representations in terms of $q$-series, which have been obtained in the case of the sunrise graph, cf. [BLOCH2, Adams:2014vja, Adams:2015gva, Adams:2015ydq, Adams:2016xah], before. More general representations are needed for the integrals considered in the present paper and we describe the necessary extension.

In performing a higher loop calculation, in intermediary steps usually more complicated nested integrals and sums occur than in the final result444For a simple earlier case, see e.g. [Vermaseren:2005qc, Ablinger:2010ty].. The various necessary decompositions of the problem that have to be performed, such as the integration by parts reduction and others, account for this in part. It appears therefore necessary to have full control on these occurring structures first, which finally may simplify in the result. Moreover, experience tells that in more general situations, more and more of these structures survive, cf. [Ablinger:2014nga] in comparison to [Vermaseren:2005qc]. If the mathematical properties of the quantities occurring are known in detail, various future calculations in the field will be more easily performed.

The paper is organized as follows. In Section 2 we present the linear systems of first order differential equations for master integrals in Ref. [Grigo:2012ji] which cannot be solved in terms of iterated integrals. We first perform a decoupling into a scalar second order equation and an associated equation for each system. Using the algorithm of Ref. [Ablinger:2015tua], the non-iterative solution both in $x$- and Mellin-$N$-space is uniquely established. In Section 3 we first determine the homogeneous solutions of the second order equations, which turn out to be ${}_{2}F_{1}$ solutions [VANH1] and obey representations in terms of weighted complete elliptic integrals of first and second kind at rational argument. In Section 4 we derive the solutions in the inhomogeneous case, which are given by iterated integrals in which some letters are given by a higher transcendental function defined by a non-iterative, i.e. definite, integral in part. We present numerical representations for $x\in[0,1]$ deriving overlapping expansions around $x=0$ and $x=1$. The methods presented apply to a much wider class of functions than the ones being discussed here specifically. These need neither to have a representation in terms of elliptic integrals, nor of just a ${}_{2}F_{1}$ function. The respective letter can be given by any multiple definite integral.

Owing to the fact that we have elliptic solutions in the present cases we may also try to cast the solution in terms of series in the nome

[TABLE]

where

[TABLE]

denotes the ratio of two complete elliptic integrals of first kind and $z(x)$ is a rational function associated to the elliptic curve of the problem. It is now interesting to see which closed form solutions the corresponding series in $q$ obey. All contributing quantities can be expressed in terms ratios of the Dedekind $\eta(\tau)$ function [DEDEKINDeta], cf. Eq. (6.14). However, various building blocks are only modular forms [KF, KOECHER, MILNE, RADEMACHER, DIAMOND, SCHOENENBERG, APOSTOL, KOEHLER, ONO, MIYAKE, SERRE, MART1] up to an additional factor of

[TABLE]

We seek in particular modular forms which have a representation in terms of Lambert–Eisenstein series [LAMBERT, EISENSTEIN] and can thus be represented by elliptic polylogarithms [ELLPOL]. However, the $\eta$-factor (1.4) in general remains. Thus the occurring $q$-integrands are modulated by a $q$-factorial [SLATER, ERDELYI1] denominator.

Structures of the kind for $k>0$ are frequent even in the early literature. A prominent case is given by the invariant $J$, see e.g. [HURWITZ],

[TABLE]

with $G_{2}(q)$ an Eisenstein series, cf. Eq.(LABEL:eq:EISEN1), and the discriminant $\Delta$

[TABLE]

In the more special case considered in [BLOCH2, Adams:2014vja, Adams:2015gva, Adams:2015ydq, Adams:2016xah] terms of this kind are not present.

For the present solutions, we develop the formalism in Section 5. We discuss possible extensions of integral classes to the present case in Section 6 and of elliptic polylogarithms [ELLPOL], as has been done previously in the calculation of the two-loop sunrise and kite-diagrams [BLOCH2, Adams:2014vja, Adams:2015gva, Adams:2015ydq, Adams:2016xah]. Here the usual variable $x$ is mapped to the nome $q$, expressing all contributing functions in the new variable. This can be done for all the individual building blocks, the product of which forms the desired solution. Section 7 contains the conclusions.

In Appendix LABEL:sec:8 we briefly describe the algorithm finding for second order ordinary differential equations ${}_{2}F_{1}$ solutions with a rational function argument. In Appendix LABEL:sec:10 we present for convenience details for the necessary steps to arrive at the elliptic polylogarithmic representation in the examples of the sunrise and kite integrals [BLOCH2, Adams:2014vja, Adams:2015gva, Adams:2015ydq, Adams:2016xah]. Here we compare some results given in Refs. [BLOCH2] and [Adams:2014vja]. In Appendix LABEL:sec:11 we list a series of new sums, which simplify the recent results on the sunrise diagram of Ref. [Ananthanarayan:2016pos].

In the present paper we present the results together with all necessary technical details and we try to refer to the related mathematical literature as widely as possible, to allow a wide community of readers to apply the methods presented here to other problems.

2 The Differential Equations

The master integrals considered in this paper satisfy linear differential equations of second order

[TABLE]

with rational functions $r(x)=p(x),q(x)$, which may be decomposed into

[TABLE]

The homogeneous equation is solved by the functions $\psi_{1,2}^{(0)}(x)$, which are linearly independent, i.e. their Wronskian $W$ obeys

[TABLE]

The homogeneous Eq. (2.1) determines the well-known differential equation for $W(x)$

[TABLE]

which, by virtue of (2.2), has the solution

[TABLE]

normalizing the functions $\psi_{1,2}^{(0)}$ accordingly. A particular solution of the inhomogeneous equation (2.1) is then obtained by Euler-Lagrange variation of constants [EULLAG]

[TABLE]

with

[TABLE]

and two constants $C_{1,2}$ to be determined by special physical requirements. We will consider indefinite integrals for the solution (2.6), which allows for more singular integrands. For the class of differential equations under consideration, $N(x)$ can be expressed by harmonic polylogarithms and rational functions, $W(x)$ is a polynomial, and the functions $\psi_{1,2}^{(0)}(x)$ turn out to be higher transcendental functions, which are even expressible by complete elliptic integrals in the cases considered here. Therefore Eq. (2.6) constitutes a nested integral of known functions [Remiddi:1999ew, Moch:2001zr, Ablinger:2013cf, Ablinger:2011te, Ablinger:2014bra] and elliptic integrals at rational argument.

We consider the systems of differential equations [Grigo:2012ji] for the $O(\varepsilon^{0})$ terms of the master integrals

[TABLE]

and

[TABLE]

with

[TABLE]

By applying decoupling algorithms [Zuercher:94, ORESYS, NewUncouplingMethod] one obtains the following scalar differential equation

[TABLE]

The above differential equations of second order contain more than three singularities. We seek solutions in terms of Gauß’ hypergeometric functions with rational arguments, following the algorithm described in Appendix LABEL:sec:8. It turns out that these differential equations have ${}_{2}F_{1}$ solutions.

Two more master integrals are obtained as integrals over the previous solutions. They obey the differential equations

[TABLE]

with $\zeta_{k},k\in\mathbb{N},\leavevmode\nobreak\ k\geq 2$ the values of the Riemann $\zeta$-function at integer argument and the harmonic polylogarithms $H_{\vec{a}}(x)$ are defined by [Remiddi:1999ew]

[TABLE]

Subsequently, we will use the shorthand notation $H_{\vec{a}}(x)\equiv H_{\vec{a}}$. The harmonic polylogarithms occurring in the inhomogeneities of Eqs. (2.34, 2.35) can be rewritten as polynomials of

[TABLE]

cf. [Blumlein:2003gb].

3 Solution of the homogeneous equation

In the following we will derive the solution of the homogeneous part of Eqs. (2, 2.32) as examples in detail, using the algorithm outlined in Ref. [IVH], see also Appendix LABEL:sec:8.

The homogeneous solutions of Eq. (2) read

[TABLE]

with

[TABLE]

The ${}_{2}F_{1}$ solutions (3.1, 3.2) are close integer series [VANH1] obeying

[TABLE]

with $c=27$. The Wronskian for this system is

[TABLE]

The solutions are shown in Figure 1.

Equivalent solutions are found by applying relations due to triangle groups [TAKEUCHI], see Appendix LABEL:sec:8,

[TABLE]

where

[TABLE]

These solutions have the Wronskian (3.5) up to a sign555The sign can be adjusted by $\psi_{1b}^{(0)}\leftrightarrow\psi_{2b}^{(0)}$. but differ from those in (3.1, 3.2). The ratios of the homogeneous solutions are given by

[TABLE]

The hypergeometric functions appearing in (3.6, 3.7) are given in terms of complete elliptic integrals [TRICOMI]

[TABLE]

Both solutions obey (3.4) with $c=16$.

We also used the relation [HYP1]

[TABLE]

noting that it is always possible to map a ${}_{2}F_{1}(a,b;c;x)$ function with $2a,2b,c\in\mathbb{Z},c>0$ into complete elliptic integrals. Their integral representations in Legendre’s normal form [LEGENDRE] read

[TABLE]

In going from (3.1, 3.2) to (3.6, 3.7) also a contiguous relation had to be applied, leading to a linear combination of two hypergeometric functions. The solutions are shown in Figure 2.

The ratio $\psi_{1b}^{(0)}/\psi_{2b}^{(0)}$ exhibits the interesting form

[TABLE]

with

[TABLE]

Whether ${}_{2}F_{1}$ solutions emerging in single scale Feynman integral calculations as solutions of differential equations for master integrals are always of the class to be reducible to complete elliptic integrals a priori is not known. However, one may use the algorithm given in Appendix LABEL:sec:8 to map a solution to one represented by elliptic integrals, if the parameters of the respective ${}_{2}F_{1}$ solution match the required pattern.

The homogeneous solutions of (2.32) read

[TABLE]

with

[TABLE]

The argument $1-z$ appeared already in complete elliptic integrals by A. Sabry in Ref. [SABRY], Eq. (68), with $x=-\lambda$, calculating the so-called kite-integral at 2 loops, 55 years ago; see also Ref. [Laporta:2004rb], Eq. (A.11), for the sunrise-diagram and [Remiddi:2016gno], Eq. (D.18), with $x=1/\sqrt{u}$ for the kite-diagram. The latter aspect also shows the close relation between the elliptic structures appearing for both topologies, which has been mentioned in Ref. [Adams:2016xah].

Using the Legendre identity [LEGENDRE]

[TABLE]

one obtains the Wronskian of the system (3.18, 3.19)

[TABLE]

cf. (2.5).

The homogeneous solutions (3.18, 3.19), which are complex for $x\in[0,1]$, are shown in Figure 3. The real part of $\psi_{3}^{(0)}(x)$ has a discontinuity at $x=1/3$ moving from $-(4/9)\sqrt{2/(3\pi)}$ to $(4/9)\sqrt{2/(3\pi)}$, while Re$(\psi_{4}^{(0)}(x))$ vanishes for $x>1/3$. Likewise, Im$(\psi_{4}^{(0)}(x))$ vanishes for $0\leq x\leq 1/3$.

We finally consider the ratio $\psi_{3}^{(0)}/\psi_{4}^{(0)}$

[TABLE]

where

[TABLE]

This structure is the same as in (3.16, 3.17) up to the pre-factor.

The solution of the inhomogeneous equations (2, 2.32) are obtained form (2.6) specifying the constants $C_{1,2}$ by physical requirements. The previous calculation of the corresponding master integrals in [Grigo:2012ji] used expansions of the propagators [Q2EXP, Steinhauser:2000ry], obtaining series representations around $x=0$ and $x=1$. The first expansion coefficients of these will be used to determine $C_{1}$ and $C_{2}$. The inhomogeneous solutions are given by

[TABLE]

with

[TABLE]

Eq. (3.25) is an integral which cannot be represented within the class of iterative integrals. It therefore requires a generalization. We present this in Section 4. Efficient numerical representations using series expansions are given in Section 5.

4 Iterated Integrals over Definite Integrals

The elliptic integrals (3.14, 3.15) cannot be rewritten as integrals in which their argument $x$ only appears in one of their integral boundaries.666Iterative non-iterative integrals have been introduced by the 2nd author in a talk on the 5th International Congress on Mathematical Software, held at FU Berlin, July 11-14, 2016, with a series of colleagues present, cf. [ICMS16]. Therefore, the integrals of the type of Eq. (3.25) do not belong to the iterative integrals of the type given in Refs. [Remiddi:1999ew, Ablinger:2011te, Ablinger:2013cf, Ablinger:2014bra] and generalizations thereof to general alphabets, which have the form

[TABLE]

For a given difference equation, associated to a corresponding differential equation, the algorithms of [SIG1, SIG2] based on [SIGBASE] allow to decide whether or not the recurrence is first order factorizable. In the first case the corresponding nested sum-product structure is returned. In the case the problem is not first order factorizable, integrals will be introduced whose integrands depend on variables that cannot be moved to the integration boundaries and over which one will integrate by later integrals. This is the case if the corresponding quantity obeys a differential equation of order $m\geq 2$, not being reducible to lower orders. Examples of this kind are irreducible Gauß’ ${}_{2}F_{1}$ functions, to which also the complete elliptic integrals ${\bf E}(z)$ and ${\bf K}(z)$ belong.

The new iterative integrals are given by

[TABLE]

and cases in which more than one definite integral $F_{m}$ appears. Here the $f_{a_{i}}(y)$ are the usual letters of the different classes considered in [Remiddi:1999ew, Ablinger:2011te, Ablinger:2013cf, Ablinger:2014bra] multiplied by hyperexponential pre-factors

[TABLE]

and $F[r(y)]$ is given by

[TABLE]

such that the $y$-dependence cannot be transformed into one of the integration boundaries completely. We have chosen here $r(y)$ as a rational function because of concrete examples in this paper, which, however, is not necessary. Specifically we have

[TABLE]

The new iterated integral (4.2) is not limited to the emergence of the functions (4.5). Multiple definite integrals are allowed as well. They emerge e.g. in the case of Appell-functions [HYP, SLATER] and even more involved higher functions. These integrals also obey relations of the shuffle type w.r.t. their letters $f_{a_{m}}(y_{m})(F_{m}[r(y_{m})])$, cf. e.g. [SHUF, Blumlein:2003gb].

Within the analyticity region of the problem one may derive series expansions of the corresponding solutions around special values, e.g. $x=0,x=1$ and other values to map out the function for its whole argument range. In many cases, one will even find convergent, widely overlapping representations, which are highly accurate and provide a numerical solution in terms of a finite number of analytic expansion coefficients. We apply this method to the solution of the differential equations in Section 2 in the following section and return to the construction of a closed form analytic representation using $q$-series and Dedekind $\eta$ functions in Section 6.

5 The Solution of the Inhomogeneous Equation by Series Expansion

The inhomogeneous solutions of type (3.25) can be expanded into series around $x=0$ and $x=1$ analytically using computer algebra packages like Mathematica or maple. One either obtains Taylor series or superpositions of Taylor series times a factor $\ln^{k}(x),k\in\mathbb{N}$. For all solutions both expansions have a wide overlap777This technique has also been used in Ref. [TLSR2a]. and one may obtain in this way a highly accurate representation of all solutions in the complete region $x\in[0,1]$.

In the following we present the first terms of the series expansion for the functions $f_{8(9,10),a(b)}(x)$ around $x=0$ and $x=1$.

For $f_{8a}$ we obtain

[TABLE]

around $x=0$. Here we also applied a series of relations for $\psi^{(k)}$-functions at rational argument, cf. Ref. [Ablinger:2011te].

Likewise, one may expand around $y=1-x=0$. In this case, we can rewrite the inhomogeneous solution given in (3.25) as

[TABLE]

with

[TABLE]

One obtains

[TABLE]

The solution of Eq. (2) around $x=0$ reads

[TABLE]

The corresponding expansion around $x=1$ is given by

[TABLE]

Here the integration constants $C_{1,2}$ and $\overline{C}_{1,2}$ are888We thank P. Marquard for having provided all the necessary constants and a series of expansion parameters for the solutions given in Ref. [Grigo:2012ji] in computer readable form.

[TABLE]

The solution of (2.34) is an integral containing the functions $f_{8a}(x)$ and $f_{9a}(x)$. Its series around $x=0$ reads

[TABLE]

Likewise, one obtains the expansion around $x=1$, which is given by

[TABLE]

with

[TABLE]

with [Broadhurst:1991fi]

[TABLE]

as integration constants in this case.

The series expansion of the solution of Eq. (2.32) is given by

[TABLE]

The solution of Eq. (2.33) reads

[TABLE]

Here the constants in (3.25) have been fixed by comparing to the first expansion coefficients in [Grigo:2012ji]

[TABLE]

The solution of (2.35) is given as an integral containing $f_{8(9)b}(x)$ with the constant

[TABLE]

and reads

[TABLE]

The corresponding solutions around $x=1$ have the expansions

[TABLE]

[TABLE]

[TABLE]

with the integration constants

[TABLE]

obtained by comparing again to the first expansion coefficients in [Grigo:2012ji]. The constants are complex here.

The solutions are illustrated in Figures 4–9. The expansions around $x=0$ and $x=1$ have wide overlapping regions in all cases. We use expansions up to $O(x^{50})$ and $O(y^{50})$, respectively. Due to the constants $C_{i}(\overline{C}_{i}),i=1,2,3,$ which are imposed by the physical case studied, all solutions are real in the region $x\in[0,1]$. The fact, that the homogeneous solutions in the cases $b$, have a branch point, has, however, consequences for the solutions around $x=0$, as will be shown below.

The function $f_{8a}(x)$ is shown in Figure 4. Its boundary values at $x=0,1$ read

[TABLE]

At very small $x$, the expansion around $x=1$ delivers too small values, while at large $x$ the small $x$ expansion evaluates to somewhat larger values, however, well below double precision. $f_{9a}(x)$ is shown in Figure 5 with the values

[TABLE]

at $x=0,1$ and a very similar behaviour for the approximation around $x=0$ and 1 as in the case of $f_{8a}$. Figure 6 shows the function $f_{10a}$, for which the boundaries are

[TABLE]

Here somewhat larger deviations of the series solutions around $x=0$ at 1 and $x=1$ at 0 are visible.

In Figure 7 the behaviour of $f_{8b}(x)$ is illustrated. The series expansion around $x=0$ starts to diverge at $x\sim 0.4$, while the expansion around $x=1$ still holds at $x\sim 0.1$. The boundary values of $f_{8b}$ at $x=0,1$ are

[TABLE]

There is a numerical artefact in Figure 7b at $x\sim 0.14$ implied by the zero-transition of $f_{8b}$ in this region.

A similar behaviour to that of $f_{8b}$ is exhibited by $f_{9b}(x)$, shown in Figure 8. Again the series-solution around $x=0$ starts to diverge for $x\sim 0.4$. However, the one around $x=1$ holds even below $x\sim 0.1$. The boundary values of $f_{9b}$ at $x=0,1$ are

[TABLE]

$f_{10b}(x)$ is shown in Figure 9. The validity of the serial expansions around $x=0$ and 1 are very similar to the cases of $f_{8(9)b}(x)$, discussed above.

The boundary values at $x=0,1$ are

[TABLE]

Notice that the representations (2.32, 5.2) allow for the analytic determination of the $N$th expansion coefficient of the corresponding series around $x=0\leavevmode\nobreak\ (y=0)$ using the techniques of the package HarmonicSums.m [HARMONICSUMS, Ablinger:PhDThesis, Ablinger:2011te, Ablinger:2013cf, Ablinger:2014bra].

The series expansions agree with those obtained by solving the differential equations through series Ansätze in [Grigo:2012ji]. In an attachment to this paper, we present the expansion of the solutions around $x=0$ and $x=1$ up to terms of $O(x^{50})$ for further use. The solutions are well overlapping in wider ranges in $x$. In the case of the functions $f_{8(9,10)a}(x)$ the power series expansion around $x=1$ reflects the branch point at $x=1/3$ in the homogeneous solution. Our general expressions easily allow expansions around other fixed values of $x$, which may be useful in special numerical applications.

The above representations constitute a practical analytic solution in the case of iterative non-iterative integrals. Indeed it applies to the whole class of these functions within their analyticity regions. Thus the method is not limited to cases in which elliptic integrals contribute. Since, however, the case in which ${}_{2}F_{1}$ solutions may be related in a non-trivial manner, see ii) and iii) in Section 6, to solutions through elliptic integrals with rational argument is very frequent, we turn now to a more detailed discussion of this case.

6 Elliptic Solutions

As we have seen, in special cases the solutions of a second order differential equation having a ${}_{2}F_{1}$ solution may be expressed in terms of the complete elliptic integrals ${\bf E}(r(z))$ and ${\bf K}(r(z))$. Our general goal is to represent the emerging structures in terms of $q$-series with explicit predicted expansion coefficients in closed form as far as possible, if not even simpler representations can be found.

Different levels of complexity can be distinguished, depending on the structure of $r(z)$ and whether only elliptic integrals of the first kind or also of the second kind necessarily contribute. Furthermore, there are requirements to other building blocks emerging in the solutions, which we will discuss below.

1. (i)

If the complete elliptic integrals are given by ${\bf K}(z)$ or ${\bf K}(1-z)$, choosing the case $z\in[0,1]$, and similarly for ${\bf E}$, one may solve the difference equation, obtained from the differential equation by a Mellin transform. It turns out that this difference equation factorizes to first order, unlike the differential equation in $x$-space; see [vonManteuffel:2017hms] for an example. The Mellin transforms (1.1) are given by

[TABLE]

since

[TABLE]

Here the Mellin convolution is defined by

[TABLE]

Eqs. (6.1) and (6.2) are hypergeometric terms in $N$, which has been shown already in Ref. [Ablinger:2013eba] for ${\bf K}(1-z)$, see also [Ablinger:2014bra]. As we outlined in Ref. [Ablinger:2015tua] the solution of systems of differential equations or difference equations can always be obtained algorithmically in the case either of those factorizes to first order. The transition to $z$-space is then straightforward. In $z$-space also the analytic continuation to the other kinematic regions is performed. 2. (ii)

In a second set of cases, only the elliptic integrals ${\bf K}(r(z))$ and ${\bf K}^{\prime}(r(z))$ contribute, with $r(z)$ a rational function. In transforming from $z$- to $q$-space, furthermore, no terms in the solution emerge which cannot be expressed in terms of modular forms [KF, KOECHER, MILNE, RADEMACHER, DIAMOND, SCHOENENBERG, APOSTOL, KOEHLER, ONO, MIYAKE, SERRE, MART1], except terms $\propto\ln^{k}(q),k\leavevmode\nobreak\ \in\leavevmode\nobreak\ \mathbb{N}$. This is the situation e.g. in Refs. [BLOCH2, Adams:2014vja, Adams:2016xah]. We will show below that here both the homogenous solution and the integrand of the inhomogeneous solution can be expressed by Lambert–Eisenstein series [LAMBERT, EISENSTEIN], also known as elliptic polylogarithms, modulo eventual terms $\ln^{k}(q)$. The remaining $q$-integral in the inhomogeneous term can be carried out in the class of elliptic polylogarithms [ELLPOL], see [Adams:2016xah]. 3. (iii)

In the cases presented in Section 3, the solutions depend both on the elliptic integrals ${\bf K}(r(z))$, ${\bf E}(r(z))$ and ${\bf K}^{\prime}(r(z))$, ${\bf E}^{\prime}(r(z))$, see also Section LABEL:sec:6.2. Both ${\bf E}(r(z))$ and ${\bf E}^{\prime}(r(z))$ can be mapped to modular forms representing them by the nome $q$ according to Eqs. (1.2, 1.3), powers of $\ln(q)$, and polylogarithms, like ${\rm Li}_{0}(q)$ [GH], and the $\eta$-factor given in Eq. (1.4). These aspects lead to a generalization w.r.t. the cases covered by ii), since in a series of building blocks the factor $1/\eta^{k}(\tau)$ has to be split off to obtain a suitable modular form. This factor is a $q$-Pochhammer symbol and also emerges in the $q$-integral in the inhomogeneous solution.

Since the topic of analytic $q$-series representations is a very recent one and it is only on the way to be algorithmized and automated for the application to a larger number of cases appearing in Feynman parameter integrals, we are going to summarize the necessary definitions and central properties for a wider audience in Section 6.1. Then we will show in Section LABEL:sec:6.2 that in the case of the differential equations (2, 2.32) both the elliptic integrals K and E are contributing, which implies the appearance of the additional $\eta$-factor (1.4). In Section LABEL:sec-52 we will then construct the building blocks for the homogeneous and inhomogeneous solutions of all terms through polynomials of $\eta$-weighted Lambert–Eisenstein series, referring to the examples (3.18, 3.19). Here we use methods of the theory of modular functions and modular forms.

6.1 From Elliptic Integrals to Lambert–Eisenstein Series

There are various sets of functions which can be used to express the complete elliptic integrals and their inverse, the elliptic functions, which have been worked out starting with Euler [EUL3], Legendre [LEGENDRE] and Abel [ABEL], followed by Jacobi’s seminal work [JAC1, JAC2] and the final generalization by Weierstraß [WEIER]999For $q$-expansions starting with the Weierstraß’ $\wp$ and $\sigma$ functions see e.g. [LANG2].. We first present a collection of relations out of the theory of elliptic integrals, their related functions and modular forms [KF, KOECHER, MILNE] for the convenience of the reader. They are essential to derive integrals over complete elliptic integrals at rational arguments, which can be represented in terms of elliptic polylogarithms. Later, we will consider the different steps for a representation of the inhomogeneous solution based on the homogeneous solutions $\psi_{3}$ and $\psi_{4}$ given before.

We first summarize a series of properties of Jacobi $\vartheta_{i}$ and the Dedekind $\eta$ functions in Section 6.1.1, followed by the representation of the complete elliptic integrals of the first and second kind by the parameters of the elliptic curve and by the Jacobi $\vartheta_{i}$ and the Dedekind $\eta$ functions in Section 6.1.2. Basic facts about modular functions and modular forms are summarized in Section 6.1.3 for the later representation of the building blocks of the homogeneous and inhomogeneous solutions of the second order differential equations of Section 2. In Section LABEL:sec:6.1.5 we collect some relations on elliptic polylogarithms and give representations of $\eta$-ratios in terms of modular forms weighted by a factor $1/\eta^{k}(\tau)$ in Section LABEL:sec:6.1.4. The modular forms are expressed over bases formed by Lambert–Eisenstein series and products thereof.

6.1.1 The Jacobi $\vartheta_{i}$ and Dedekind $\eta$ Functions

As entrance point we use Jacobi’s $\vartheta_{i}$ functions [JAC2]. The $\vartheta$ functions possess $q$-series and product representations101010In the literature different definitions of the Jacobi $\vartheta$-functions are given, cf. [TRICOMI], p. 305. We follow the one used by Mathematica.

[TABLE]

The elliptic polylogarithms, introduced in (LABEL:eq:ELP1, LABEL:eq:ELP2) below are also $q$-series, containing a specific parameter pattern which allows to accommodate certain classes of $q$-series emerging in Feynman integral calculations. The product representations associated to (6.7–6.9) read

[TABLE]

They are closely related to Euler’s totient function [EUL4]

[TABLE]

the first emergence of $q$-products, and to Dedekind’s $\eta$ function [DEDEKINDeta].111111The $\vartheta$ and $\eta$ functions, as well as their $q$-series, play also an important role in other branches of physics, as e.g. in lattice models in statistical physics in form of Rogers-Ramanujan identities, see e.g. [BAXTER, MCCOY, BOSTAN], percolation theory [KZ], and other applications, e.g. in attempting to describe properties of deep-inelastic structure functions [Scott:1999nm]. In the latter case, the asymptotic behavior of Dedekind’s $\eta$ function at $x\sim 1$ seems to resemble the structure function for a wide range down to $x\leavevmode\nobreak\ \sim\leavevmode\nobreak\ 10^{-5}$. It has a surprisingly similar form as the small-$x$ asymptotic wave equation solution [WAVE], however, with a rising power of the soft pomeron [Donnachie:1996gq].

[TABLE]

One has121212It is usually desirable to work with $\eta$-functions depending on integer multiples of $\tau$ only, cf. [KOEHLER], which can be achieved by rescaling the power of $q$.

[TABLE]

In the following we will make use of series representations of both Jacobi $\vartheta$- and Dedekind $\eta$-functions. We list a series of important relations for convenience :

[TABLE]

Many other identities hold and can be found e.g. in Refs. [KOEHLER, MACDONALD, KAC, LEPOWSKY, ZUCKER, MARTIN, WILLIAMS, KENDRIL, ONO].

6.1.2 Representations of the Modulus and the Elliptic Integrals

For later use we consider also the structure of the differential equation of the Weierstraß’ function $\wp(z)$ [WEIER],

[TABLE]

The functions $g_{2},g_{3},e_{1},e_{2}$ and $e_{3}$ are given by

[TABLE]

and the following representation in terms of Jacobi $\vartheta$ functions holds:

[TABLE]

Here Jacobi’s identity is implied by (6.31) with

[TABLE]

The r.h.s. of (6.28) parameterizes the elliptic curve

[TABLE]

of the corresponding problem. Setting $e_{1}-e_{2}=1$ for the purpose of illustration, the elliptic curves corresponding to the module $k$ is shown in Figure 10, choosing specific values.

The modulus $k$ can be represented in terms of the functions $e_{i}$ by

[TABLE]

cf. (3.8, 3.20). $k$ and $k^{\prime}=\sqrt{1-k^{2}}$ are given by

[TABLE]

cf. (6.28), which implies the following relation for $\eta$ functions

[TABLE]

Further, one may express the elliptic integral of the first kind K by

[TABLE]

Sometimes one also introduces the Jacobi functions $\omega,\omega^{\prime},\eta$ and $\eta^{\prime}$, which are defined by

[TABLE]

with

[TABLE]

The function $\eta^{\prime}$ can be obtained using Legendre’s identity (3.21) in the form

[TABLE]

One obtains the following representations of the elliptic integrals of the second kind by

[TABLE]

Later on we will use the relation [ERDELYI2, ABRSTE] for E

[TABLE]

and the Legendre identity (3.21) to express E′**,

[TABLE]

6.1.3 Modular Forms and Modular Functions

All building blocks forming the homogeneous solutions and the integrand of the inhomogeneous solutions of the second order differential equations considered above can be expressed in terms of $\eta$-ratios. They are defined as follows.

Definition 6.1**.**

Let $r=\left(r_{\delta}\right)_{\delta|N}$ be a finite sequence of integers indexed by the divisors $\delta$ of $N\in\mathbb{N}\backslash\{0\}$. The function $f_{r}(\tau)$

[TABLE]

is called $\eta$-ratio.

These are modular functions or modular forms; the former ones can be obtained as the ratio of two modular forms. In the following we summarize a series of basic facts on these quantities in a series of definitions and theorems needed in the calculation of the present paper, cf. also Refs. [KF, KOECHER, MILNE, RADEMACHER, DIAMOND, SCHOENENBERG, APOSTOL, KOEHLER, ONO, MIYAKE, SERRE, MART1].